Boundary conditions for an FEM approximation of the Laplace operator

Using FEM, I want to approximate the Laplacian

$$u = nabla cdot nabla h , ,$$

where $h(x,y)$ is an FEM approximated scalar field on the same mesh, i.e. piecewise differentiable.

I am using MOOSE to solve the following (hopefully correct!) weak formulation of the above equation:

$$int_Omega uw_i , text{d}Omega = -int_Omega nabla h cdot nabla w_i , text{d}Omega + int_Gamma nabla h cdot vec{n}w_i , text{d}Gamma ,,$$

where $w_i$ is the i-th test function and $vec{n}$ is the normal vector on the boundary $Gamma$ of the computational domain $Omega$. I am not using any boundary conditions on $u$.

The FEM solution that I obtain for this weak formulation is inaccurate at the boundaries (see below). I have tested this by prescribing the field $h=text{sin}(x)$ on a 2D rectangular domain from $-2pi le x le 2pi$
and solving with quadratic Lagrange shape functions for both $h$ and $u$. (it is also inaccurate using linear shape functions, but it is easier to spot with quadratic shape functions).

Is it possible to obtain an accurate $u$ at the boundaries with this approach? Am I missing a boundary condition/contribution?

---

Background:

I am dealing with a thin film approximation of the Navier-Stokes equations in 2D FEM, where the film height is modeled by a dependent variable $h(x,y)$. In order to include surface tension as a contribution to the pressure gradient, I want to find the curvature of the film surface $u(x,y)$ that I approximate by the second derivative

$u approx nabla cdot nabla h$ .

My idea is to introduce curvature $u$ as a variable, solve the above equation by FEM, and couple $u$ into the momentum equation.