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 $-2\pi \le x \le 2\pi$ 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.
