Unfitted Finite Element Discretization
- Unfitted finite element discretization is a technique where the mesh is independent of the domain geometry, using implicit or explicit representations to define boundaries.
- It employs weak formulations and specialized numerical quadrature over cut cells to achieve high-order accuracy, even for complex interface problems.
- Stabilization methods such as ghost penalties and cell aggregation are integrated to maintain optimal convergence rates and robust performance in simulations.
Unfitted finite element discretization refers to a class of finite element techniques in which the computational mesh does not conform to the geometry of the domain or interfaces present in the problem. Instead, the domain boundary, interfaces, or embedded surfaces are described independently—often implicitly via level set functions or explicitly by surface meshes—while the background finite element mesh (typically structured or quasi-uniform) remains fixed and independent of the geometry. Integration over physical domains or surfaces is handled by numerical quadrature over the intersected ("cut") portions of mesh cells. Unfitted finite element methods (often termed CutFEM, XFEM, embedded, or immersed finite element methods) are designed to decouple geometry handling from mesh generation, enabling efficient simulation on complex or evolving domains and for problems with interfaces or boundaries of arbitrary shape.
1. Geometric Representation and Extension Principles
A standard approach in unfitted finite element methods is to give the domain boundary, interface, or surface Γ implicitly as the zero level set of a function φ, i.e., . Rather than requiring the mesh to align with Γ, one works with a background mesh covering a "narrow band" or encompassing the full domain, using φ to define which regions are active.
For hypersurface PDEs, a common methodology is to "lift" the surface equation to a neighborhood of Γ in the embedding domain, often by extending data and operators along the normal direction. For example, as in (Olshanskii et al., 2014), one extends the PDE to a tubular neighborhood using an extension that ensures solutions remain constant in the normal direction.
Explicit boundary representations are also supported. In the presence of explicit oriented surface meshes (e.g., triangulated STL data from CAD), algorithms robustly intersect the mesh cells with the boundary to define cut cells for integration and enable simulation on non-level-set geometries (Badia et al., 2021).
2. Weak Formulation and Discretization Strategy
Unfitted finite element discretizations proceed by formulating the weak (variational) form of the PDE on the physical domain or surface, then integrating this form over each cut element or intersected cell. For surface or interface problems, the integral is extended to a neighborhood or narrow band, as in the extended Laplace–Beltrami formulation:
which is equivalently extended to the embedding domain as a uniformly elliptic PDE:
where is a determinant of a scaling matrix involving the Hessian of φ (Olshanskii et al., 2014).
The finite element space is typically built from continuous piecewise polynomials defined on the background mesh. The computational mesh is "unfitted," i.e., it does not conform to the boundary or interface. Integrals in the variational form are evaluated only on the intersection of mesh elements with the physical domain or on a narrow band.
High order accuracy can be achieved by approximating the geometry with higher order accuracy. This is realized through parametric mappings that "lift" a low-order (piecewise planar) interface representation to a high-order geometry, using the gradient of the level set and correction functions determined by solving local 1D nonlinear equations (Lehrenfeld et al., 2016, Lehrenfeld et al., 2016, Lederer et al., 2016).
3. Numerical Quadrature over Implicit and Explicit Domains
Numerical integration on implicitly or explicitly defined cut domains or surfaces is a significant technical challenge. Several strategies have been developed (Olshanskii et al., 2016, Badia et al., 2021, Martorell et al., 2023):
- Subdivision (Sub-triangulation): The intersection of the cut cell with the domain is approximated by subdividing into simpler shapes (triangles, tetrahedra), and conventional quadrature is applied.
- Moment-fitting: Quadrature weights are solved for explicitly to satisfy moment conditions for polynomials up to a prescribed degree, enabling optimal-order integration with minimal points.
- Local quasi-parametrization: The cut region's curvilinear boundary is parametrized locally (as a graph), enabling accurate 1D quadrature.
- Monte Carlo integration: Uniform sampling in a small region near the cut interface, though variance and complexity make this less attractive for high accuracy.
- Explicit geometric intersection (B-rep, STL): For domains given by explicit boundary representations, algorithms construct a two-level integration mesh—first computing exact intersections, then further decomposing into convex polytopes or simplex elements, enabling stable explicit quadrature (Badia et al., 2021, Martorell et al., 2023).
The choice of quadrature directly impacts both the computational efficiency and the attainable convergence rate of the unfitted FEM.
4. Error Analysis and Convergence Rates
Error analysis in unfitted finite element methods must account for three primary error sources:
- Finite element approximation (interpolation error): Controlled by the polynomial degree of the finite element space.
- Geometry approximation error: Arising from the error in approximating the boundary/interface, typically determined by the level set representation or the degree of quadrature approximation ().
- Quadrature error: Due to inaccuracies in integrating over cut elements, governed by the chosen integration scheme and its local error order .
Theoretical results show that if the level set function and its Hessian are approximated with sufficient accuracy, and the geometric quadrature is performed with order (with the spatial dimension), then optimal convergence rates in the and norms are achieved (Olshanskii et al., 2014, Olshanskii et al., 2016, Lehrenfeld et al., 2016, Lehrenfeld et al., 2016). The general form of the error estimates is:
Optimal error order can be retained provided the narrow band width (for surface extensions) and the geometry approximation both scale appropriately with the mesh size.
5. Stabilization, Conditioning, and Robustness
Cut cells—elements with small intersection with the physical domain—can cause severe ill-conditioning in the discrete systems. Unfitted finite element methods address this via:
- Aggregation of cut cells: Small cut cells are merged into larger agglomerates, and basis functions are constrained to produce robust FE spaces (AgFEM) (Badia et al., 2018, Badia et al., 2022).
- Ghost penalty stabilization: Additional penalty terms penalize jumps of normal derivatives or function values across facets between cut and interior cells (Lederer et al., 2016).
- Discrete extension operators: Construction of extension operators that propagate stable functions from interior cells to the boundary region, providing intrinsic stability and optimal interpolation properties (Burman et al., 2021).
- Multigrid solvers: Prolongation operators and interface smoothers tailored for non-nested cut spaces yield level-independent convergence, even with large coefficient contrasts in interface problems (Ludescher et al., 2018, Kothari et al., 2021).
Well-posedness, stability, and conditioning are addressed at both formulation and solver levels to guarantee that unfitted methods remain robust for arbitrarily cut configurations and even as the mesh is refined.
6. Applications and Computational Implications
Unfitted finite element discretization is broadly applicable:
- Elliptic PDEs on stationary and evolving surfaces: Via extension formulations and narrow band approaches (Olshanskii et al., 2014).
- Elliptic and parabolic interface problems: Including optimal control and Stokes interface problems, where the method allows for accurate resolution of discontinuities or singularities at interfaces (Lehrenfeld et al., 2016, Lehrenfeld et al., 2016, Lederer et al., 2016, Yang et al., 2018).
- Moving domain and evolving surface PDEs: Space-time and characteristic-Galerkin approaches enable high-order accurate simulation of advection-diffusion and other time-dependent PDEs on domains with moving boundaries (Ma et al., 2021, Lou et al., 2021, Badia et al., 2022, Reusken et al., 2 Jan 2024).
- Explicit CAD geometries: New computational pipelines accommodate B-rep surfaces, providing a direct route to simulation on industrial geometries (Badia et al., 2021, Martorell et al., 2023).
- Contact problems and FSI: Via unfitted coupling strategies with Lagrange multipliers and fictitious domain approaches, supporting independent meshing for fluid and solid regions (Alshehri et al., 5 May 2025).
Unfitted FE methods remove the mesh generation bottleneck for complex or evolving geometries and are especially advantageous in settings where meshing the exact domain is infeasible or prohibitively expensive.
7. Numerical Experiments and Practical Implementation
Numerical results document the empirical convergence rates (H¹ and L²), demonstrate optimal rates for problems on circles, spheres, tori, and complicated domains, and analyze performance for higher order elements (P₂, P₃, etc.) (Olshanskii et al., 2014, Lehrenfeld et al., 2016). Computational complexity studies highlight that quadrature methods based on local quasi-parametrization and moment fitting achieve optimal cost-per-element in 2D, while more expensive subdivision or Monte Carlo approaches are generally not competitive for high-order applications (Olshanskii et al., 2016).
Error constants, volume and surface integration accuracy, and the impact of geometry representation have been quantified for thousands of real-world models in large-scale databases (Badia et al., 2021). High-performance implementations leveraging matrix-free operator evaluation with sum-factorization and SIMD vectorization have demonstrated an order of magnitude speed-up in operator application over matrix-based methods for polynomial degrees ≈3 and higher, while maintaining robustness on large geometries and high cut-cell ratios (Bergbauer et al., 11 Apr 2024).
Future directions involve further development of scalable quadrature schemes in 3D, new mapping approaches for explicit CAD surfaces, and adaptive error estimators for localized mesh refinement near interfaces or singularities. Unfitted finite element discretization remains a primary technical enabler for advancing simulation capabilities in complex, heterogeneous, and evolving domains encountered in multiple areas of computational mathematics and engineering.