About
CutFEM (Cut Finite Element Method) is a framework that allows the finite element method to be applied even when boundaries or interfaces do not align with the background mesh. It is valuable in applications such as the level set method for topology optimization, since it enables function evaluation without remeshing. In conventional level set methods, projection functions like the Heaviside function were used to create smooth functions, but more recent approaches appear to make use of CutFEM.
Weak Form
a(u,v) := \int_\Omega \varepsilon(v):\big(\mathbb C:\varepsilon(u)\big)\,d\Omega,
\qquad \\
\ell(v) := \int_\Omega \mathbf f\cdot v\,dx + \int_{\Gamma_N} t\cdot v\,ds.Nitsche method
What is this?
To solve FEM under CutFEM condition, Nitsche method is one of the promising approach. As you can see below article, condensation is generally used for ordinal FEM.
But since this approach eliminate degree of freedom, it Dirichlet BCs are strictly enforced. In contrast, the Nitsche method incorporates boundary conditions directly into the weak formulation by adding extra terms to the variational integral, rather than eliminating degrees of freedom explicitly. Because of that, it is able to solve equation even when it is under a condition where it is difficult to be applied. This approach offers several advantages:
- It does not require large penalty parameters, unlike the penalty method.
- It does not introduce additional unknowns, unlike the Lagrange multiplier method.
- It yields a solution that is both stable and consistent.
The formulation
Add “symmetrized boundary term + stabilization term” to the weak form:
\begin{align}
a(u,v)
- \int_{\Gamma} (\sigma(u)\cdot n)(v)\,ds
- \int_{\Gamma} (\sigma(v)\cdot n)(u-g)\,ds
+ \frac{\gamma}{h}\int_{\Gamma} (u-g)(v)\,ds
= l(v)
\end{align}where
\sigma(u)\cdot n : normal component of the stress vector
\gamma/h : penalty stabilization term (proportional to the mesh width)
l(v) = \int_\Omega f v \, dx + \int_{\Gamma_N} t v \, dsSee each terms in (1) in detail.
a(u,v)
\;-\;\underbrace{\int_{\Gamma_D} (\sigma(u)\cdot n)\,v \,ds}_{\text{(2) boundary flux term}}
\;-\;\underbrace{\int_{\Gamma_D} (\sigma(v)\cdot n)\,(u-g) \,ds}_{\text{(3) symmetry term}}
\;+\;\underbrace{\frac{\gamma}{h}\int_{\Gamma_D}(u-g)\,v \,ds}_{\text{(4) stabilization/penalty term}}
= l(v)(2)
\int_{\Gamma_D} (\sigma(u)\cdot n) v- This comes from integration by parts as the natural boundary flux term.
- For Neumann boundaries, such a term would simply move to the right-hand side as $\int_{\Gamma_N} t v$.
- But on Dirichlet boundaries, instead of prescribing a traction, we want to enforce $u=g$.
- If we only keep this term, the Dirichlet condition is not properly enforced.
(3)
\int_{\Gamma_D} (\sigma(v)\cdot n)(u-g)- This is the key of Nitsche’s method.
- Term (2) alone ensures consistency, but lacks stability.
- By adding this “symmetric counterpart,” where the roles of $u$ and $v$ are exchanged, the formulation becomes symmetric and stable.
- This ensures that the Dirichlet condition is properly reflected in the weak form.
(4)
\tfrac{\gamma}{h}\int_{\Gamma_D}(u-g)v- This is the stabilization (penalty) term.
- Unlike in the pure penalty method, the coefficient does not need to be extremely large.
- With a suitable choice of $\gamma$, the method is both stable and accurate, even for finite element spaces of limited richness.
Ghost Penalty
However, in unfitted FEM there is the problem that cut elements can become extremely small.
To stabilize this, Nitsche’s method is often combined with a ghost penalty.
In particular, in the level set method, boundary conditions are imposed weakly using Nitsche’s method, while the ghost penalty is employed to stabilize the ill-conditioning of the cut elements.
As you can see in [3], Ghost Penalty is introduced in Weak Form,
\begin{align}
A(d,s) = a(d,s) + k_\psi(d,s) + j(d,s) = l(s), \quad \forall s \in T
\end{align}(1)
The standard bilinear form for elasticity. This represents the integral of the stress σ\sigmaσ and the strain ε\varepsilonε. It is the standard stiffness term in the finite element method
a(d,s) = \int_{\Omega(\varphi)} \sigma(d) : \varepsilon(s)\, dx(2)
A penalty term in the level set method, introduced to enforce zero displacement inside isolated material regions.
Here, \psi is a function that marks such regions.
k_\psi(d,s) = \int_{\Omega(\varphi)} \psi d \cdot s \, dx(3)
This is a stabilization term known as ghost penalty stabilization.
- \mathcal{T}_G: ghost skeleton (the collection of cells near the interface)
- F: facets (faces) of cells cut by the interface
- [n_F \cdot \nabla d]: jump of the gradient in the normal direction
- h_F: size of the facet
- \gamma(\lambda+\mu): scaling factor based on Lamé constants
j(d,s) = \gamma(\lambda+\mu) \sum_{F \in \mathcal{T}_G} \int_F h_F^3 \,[n_F \cdot \nabla d]\cdot[n_F \cdot \nabla s]\, ds
- In CutFEM, the mesh is cut by the interface, which may result in very small cut cells that cause numerical instabilities.
- To address this, a penalty is imposed on the jump of the normal gradient across elements cut by the interface, thereby weakly enforcing continuity of the solution.
- This prevents ill-conditioning and allows for stable computations.
Reference
[1] ChatGPT 5
[2] Burman “CutFEM: Discretizing geometry and partial differential equations”
[3] Z.J. Wegert, J. Manyer, C. Mallon, S. Badia, V.J. Challis, “Level-set topology optimisation with unfitted finite elements and automatic shape differentiation”
