Mean curvature computation in unfitted finite element methods

Computing the mean curvature of an implicitly described surface is a classical challenge in numerical geometry processing and a key ingredient in surface-PDE and geometric-flow solvers. The goal of this project is to combine two modern, complementary variational ideas for curvature — and to make them work in the unfitted (CutFEM / level-set) setting, where the surface is not aligned with the computational mesh.

Why curvature is hard to approximate

Two variational approaches (body-fitted case)

We take a setting where the surface $\Gamma$ is resolved by the mesh (body-fitted) as the reference point, and review two ways of turning the “kinks carry the curvature” idea into a robust finite element computation.

1. Weak mean-curvature vector (Hansbo–Larson–Zahedi, Lou)

The mean curvature vector is the surface Laplacian (Laplace–Beltrami operator) of the identity (position) map and is, equivalently, expressed through the surface divergence of the normal:

$\Delta_\Gamma\,\id_\Gamma \;=\; \mathbf{H} \;=\; -(\divG\,\nn)\,\nn .$

Multiplying by a (vector) test function $\mathbf{v}$ and integrating by parts over $\Gamma$ moves one derivative onto the test function (Dziuk's identity),

$\displaystyle\int_\Gamma \mathbf{H}\cdot\mathbf{v}\,\mathrm{d}s \;=\; \int_\Gamma \nabla_\Gamma\id_\Gamma : \nabla_\Gamma\mathbf{v}\,\mathrm{d}s ,$

so that only first surface derivatives of the geometry appear — which are well-defined for the piecewise-planar $\Gamma_h$, with the curvature entering through the conormal contributions at the kinks. This stabilized cut-FEM formulation was introduced by Hansbo, Larson and Zahedi [HLZ15] for the mean curvature vector. A subtlety is a regularity/stability mismatch: the right-hand side is only $H^1$-stable in the geometry while the left-hand side is only $L^2$-stable. They close this gap with a ghost-penalty term that acts as a regularization. The original analysis is for first-order elements; in Yimin Lou's work in our group [YL23] the construction is carried to higher order in the isoparametric unfitted setting. Combining and comparing these techniques is one core of this project.

2. Distributional curvature with FE lifting (Gopalakrishnan, Neunteufel et al.)

A second, very elegant route interprets the curvature distributionally and then lifts it into a finite element space [GNSW23], [GNSW24], [GN23]. On a piecewise-planar $\Gamma_h$ the smooth (element-interior) curvature vanishes, and the distributional curvature reduces to a sum of facet terms — the conormal jumps $\jump{\boldsymbol{\mu}}$ across the edges of $\Gamma_h$ (in 3D additionally vertex terms):

$\displaystyle \langle \kappa_h, \varphi_h\rangle \;=\; \sum_{E\subset\Gamma_h} \int_E \jump{\boldsymbol{\mu}}\,\varphi_h \,\mathrm{d}\ell .$

Crucially, the construction generalizes directly to higher order: for a curved (isoparametric) $\Gamma_h$ one simply keeps the “inner” terms as well, so the curvature functional becomes

$\displaystyle \langle \kappa_h, \varphi_h\rangle \;=\; \underbrace{\sum_{T\subset\Gamma_h} \int_T \kappa_T\,\varphi_h\,\mathrm{d}s}_{\text{element (smooth) part}} \;+\; \underbrace{\sum_{E\subset\Gamma_h} \int_E \jump{\boldsymbol{\mu}}\,\varphi_h\,\mathrm{d}\ell}_{\text{facet (kink) part}} \;+\; (\text{vertex part in 3D}).$

The discrete curvature $\kappa_h$ is then obtained by a single mass-matrix lifting: find $\kappa_h$ in a finite element space with $\int_{\Gamma_h}\kappa_h\,\varphi_h = \langle\kappa_h,\varphi_h\rangle$ for all $\varphi_h$. This is ordinary finite element assembly (element integrals + facet integrals) — no connectivity graph or per-vertex angle bookkeeping is built by hand — and remarkably it needs no stabilization.

Interactive comparison: naive vs. lifted

Both scenes show the same piecewise-planar (faceted) torus, coloured by mean-curvature magnitude on one shared colour scale — so the comparison isolates the method, not the geometry. On the left, the naive elementwise curvature collapses to (near) zero on every flat facet. On the right, the variational lifting recovers a non-trivial curvature field from the kinks. On this raw faceted surface the lifting is still visibly noisy — exactly the issue that the higher-order geometry and ghost-penalty stabilization discussed above are meant to fix. They are self-contained NGSolve webgui exports and load on demand — click to start the interactive 3D view.

The unfitted / level-set setting

Ghost-penalty stabilization: the open questions

Stabilization in this project plays (at least) two conceptually different roles, and disentangling them is part of the investigation:

• Normal-derivative (volume) stabilization. An anisotropic diffusion in the volume mesh acting orthogonally to the surface. It is standard in trace-FEM [ORG09], [GLR18] and makes sense of the dimensional mismatch between the (volume-defined) trace-FEM functions and the (surface) PDE — it is a mechanism for the unfitted nature of the problem and controls the conditioning.
• Tangential (regularizing) stabilization. The Hansbo–Larson–Zahedi mean-curvature approach [HLZ15] applies an additional ghost penalty that, while also assembled in the surrounding volume, acts tangentially to the surface — genuinely regularizing in the tangential spatial direction to cure the regularity/stability mismatch of the weak formulation.

The Gopalakrishnan–Neunteufel lifting needs no tangential regularization in the body-fitted case. One may therefore hope that in the unfitted setting the volume normal-derivative stabilization alone suffices. Whether that is true is genuinely open: the same regularity/stability mismatch is present, and it is not obvious (and worth understanding) how the distributional-lifting analysis avoids needing a tangential regularization in the first place.

Project roadmap

The project starts in 2D (a 1D curve embedded in the plane), where everything can be visualized and debugged cheaply, and extends to 3D surfaces afterwards.

1. Warm-up: body-fitted distributional FE lifting. Implement the distributional curvature with FE lifting [GNSW23] in NGSolve/Python on a body-fitted mesh, and verify convergence against the exact curvature.
2. Warm-up: level-set geometry in ngsxfem + weak mean curvature. Describe the geometry by a level set and build the cut $\Gamma_h$ in ngsxfem, with the isoparametric mapping [L16] for higher-order geometry; compute the weak mean-curvature vector as in [HLZ15].
3. Variational lifting on the cut curve — stabilization & robustness. On the piecewise-planar surface in 2D, assemble the lifting and compare different stabilization strategies (volume normal-derivative vs. tangential ghost-penalty regularization), with the focus on robustness across favorable and unfavorable cut configurations.
4. Convergence studies. Quantify the accuracy of the (best) variants against exact curvatures and determine the observed convergence orders.
5. Extensions to 3D and higher order. Carry the implementation to 3D surfaces and to the higher-order isoparametric geometry with the “inner” curvature terms.

References