Ghost Penalty Bilinear Forms

• Ghost penalty bilinear forms are stabilization techniques in CutFEM that penalize jumps in functions or derivatives to address issues on non-conforming meshes.
• They enable locking-free stabilization by carefully tuning penalty parameters to control ill-conditioning and discretization errors in small or cut cells.
• Recent advances include matrix-free implementations and tailored mass/stiffness penalties that improve computational efficiency in large-scale, high-order simulations.

Ghost penalty bilinear forms are stabilization mechanisms employed in cut finite element methods (CutFEM) and unfitted discretizations to ensure numerical robustness on meshes with elements that are partially inside the domain or poorly cut by the boundary. These forms act by penalizing specific quantities (typically jumps or differences of finite element functions or their derivatives) across mesh entities, such as faces or element pairs, outside or near the domain boundary. Their design addresses issues related to ill-conditioning, locking, and discretization error in situations where the mesh does not conform to the computational domain, particularly in the presence of small cut elements or isolated degrees of freedom. Recent developments have introduced generalizations enabling locking-free stabilization and efficient matrix-free evaluation suited for high-order and large-scale settings.

1. Abstract Framework for Ghost Penalty Bilinear Forms

The general theory of ghost penalty operators, as in "On the Design of Locking Free Ghost Penalty Stabilization and the Relation to CutFEM with Discrete Extension" (Burman et al., 2022), systematizes their construction on active background meshes $\mathcal{T}_h$ partitioned into “large” ($\mathcal{T}_h^L$) and “small” ($\mathcal{T}_h^S$) elements based on the volume fraction relative to the physical domain. Each small element $T \in \mathcal{T}_h^S$ is associated with a nearby large element $S_h(T) \in \mathcal{T}_h^L$. The core structure consists of connecting pairs of elements—either through faces or extension mappings—and penalizing a chosen functional of the jump in the discrete function between those elements.

Let $$b : \mathbb{P}_k(\mathbb{R}^d) \times \mathbb{P}_k(\mathbb{R}^d) \to \mathbb{R}$$ be a symmetric, positive semi-definite bilinear form and $\alpha_m \in \mathbb{R}$ chosen to mirror an $H^m$-seminorm scaling (Assumption A1). For $v, w \in V_h$, a general ghost penalty contribution reads

$$s_{m, b, T, S_h(T)}(v, w) = \tau h^{\alpha_m} b\left([v]_{T, S_h(T)}, [w]_{T, S_h(T)}\right)$$

with the full stabilization aggregated as

$$s_{h, m}(v, w) = \sum_{T \in \mathcal{T}_h^S} s_{m, b_T, T, S_h(T)}(v, w)$$

where $[v]_{T, S_h(T)}$ denotes the canonical polynomial extension jump between $T$ and $S_h(T)$. Variants using functionals, such as degrees of freedom, as the penalized quantity are encompassed within this abstract setting.

2. Principal Variants and Examples

Key practical ghost penalty forms in CutFEM and related methods include:

• Face-based Penalty: For internal faces $F$ between elements $T_1$ and $T_2$, penalizes jumps of normal derivatives,

$$s_{h, m}(v, w) = \sum_{F \in \mathcal{F}_h} \tau h^{3 - 2m} [\nabla_n v]_F [\nabla_n w]_F$$

(standard $m = 1$ variant in CutFEM).

• Element-pair Gradient and $L^2$ Penalties:
  • Gradient penalty:

  $$s_{h, 1}(v, w) = \sum_{T \in \mathcal{T}_h^S} \tau \left( \nabla [v]_{T, S_h(T)}, \nabla [w]_{T, S_h(T)} \right)_T$$ - $L^2$-jump penalty:

  $$s_{h, 0}(v, w) = \sum_{T \in \mathcal{T}_h^S} \tau h^{-2} (v_1 - v_2, w_1 - w_2)_{T \cup S_h(T)}$$

• Nodal (Degree-of-Freedom) Stabilization: Controls individual cut degrees of freedom (DOFs) via rank-one updates,

$$s_{h, 1}(v, w) = \sum_{i \in I^S} \tau h^{d-2} (v(x_i) - v^e_{S_h(T_i)}(x_i))(w(x_i) - w^e_{S_h(T_i)}(x_i))$$

Corresponding penalty parameters and exponents are calibrated to the regularity and norm of interest, e.g., $\alpha_0 = d$ for $L^2$-control, $\alpha_1 = d-2$ for $H^1$-control.

3. Analytical Properties: Stability, Consistency, and Locking

The ghost penalty stabilization achieves uniform coercivity and bounds for the discrete solution, independently of cut configuration, when Assumptions A1 (scaling) and A2 (face-neighbor connectivity) are satisfied (Burman et al., 2022). Stability is demonstrated using discrete Poincaré and chain inequalities; for example,

$$\sum h^{\alpha_m} \|v\|_b^2 \lesssim \|\nabla^m v\|^2_{\Omega} + \|v\|^2_{s_{h,m}}$$

ensures norm equivalence and well-posedness (see (Burman et al., 2022), Lemma 2.1).

A crucial aspect is the behavior of $\ker(s_{h,m})$: for classical face/bulk ghost penalties, the kernel is too restrictive and does not contain adequate approximation subspaces. As a result, these methods exhibit locking—the limiting solution as $\gamma \to \infty$ does not converge to the aggregated FE subspace, and error plateaus or worsens with excessive penalty (Badia et al., 2021). In contrast, locking-free variants penalize the distance to an extension or aggregation operator; their kernel is precisely the robust aggregated FE space ($\mathcal{V}_h^{\rm ag}$), ensuring optimal convergence and absence of locking for large penalties.

4. Ghost Penalty in Mass and Stiffness Bilinear Forms

Ghost penalty bilinear forms are employed not only in stiffness stabilization but also in mass scaling, most notably in explicit dynamics and time-stepping for immersogeometric problems. For wave and shell-type equations, the total mass and stiffness operators receive additional contributions:

• Ghost-mass penalty:

$$M_\gamma(\partial_{tt}u, v) = \sum_{F \in \mathscr{F}_{\rm ghost}} \int_F \rho \gamma_M h_K^{2p+1} [\partial_n^p \partial_{tt} u] [\partial_n^p v]\, dS$$

• Ghost-stiffness penalty:

$$K_\gamma(u, v) = \sum_{F \in \mathscr{F}_{\rm ghost}} \int_F \kappa \gamma_S h_K^{2p-1} [\partial_n^p u] [\partial_n^p v]\, dS$$

The inclusion of the ghost-mass term suppresses high-frequency localized modes on arbitrarily small-cut elements, thereby rendering the maximum eigenvalue of the semi-discrete system and the CFL condition independent of cut-cell geometries (Stoter et al., 2023). This mechanism has been shown to enable time-step sizes that are orders of magnitude larger without loss of accuracy for both membrane and shell cases.

5. Relations to Discrete Extension and Aggregated FE Methods

The kernel of suitable ghost penalty operators, such as those based on discrete extensions, coincides with the so-called aggregated FE space, wherein ill-posed DOFs are linearly dependent on well-posed ones via local extension operators (Badia et al., 2021, Burman et al., 2022). This connection underpins the construction of locking-free ghost penalty stabilizations and explains their robust performance.

Specifically, defining a projector $\mathcal{P}_h^{\rm ag}$ to $\mathcal{V}_h^{\rm ag}$ and penalizing the $L^2$ or $H^1$ distance between a function and its extension,

$$s_h(u, v) = \sum_{U \in \mathcal{T}_h^{\rm ag}} \gamma h_U^{-2} (u - \mathcal{P}_h^{\rm ag} u, v - \mathcal{P}_h^{\rm ag} v)_U,$$

yields stabilization with the aggregated FE space as its nullspace; taking $\gamma \to \infty$ recovers the strong aggregated method, ensuring locking-free behavior and optimal convergence (Badia et al., 2021).

6. Matrix-Free Implementation and Computational Efficiency

Efficient evaluation of high-order ghost penalty terms is essential for large-scale simulations. Recent advances in matrix-free approaches leverage the tensor-product structure of FE shape functions and the face-centric definition of penalties. In this paradigm, the ghost penalty operator is reduced to a sequence of low-dimensional matrix-vector multiplications confined to ghost facets (Wichrowski, 28 Feb 2025).

For polynomial degree $p$ in $d$ dimensions, the ghost penalty takes the form

$$G_h(u, v) = \gamma_A \sum_{F \in \mathcal{F}_h} \sum_{k=0}^p \frac{h_F^{2k - 1}}{(k!)^2} \int_F [\partial_n^k u] [\partial_n^k v]\, ds,$$

with per-face evaluation via Kronecker factorization into 1D penalty and mass matrices. This reduces computational complexity to $O(p^{d+1})$ per face, obviates the explicit computation of high-order derivatives, and enables robust scaling in matrix-free frameworks such as deal.II (Wichrowski, 28 Feb 2025). Benchmarks confirm that these operations contribute a moderate fraction of total computation, even at high polynomial degree and for substantial cut-cell fractions.

7. Comparative Performance, Parameter Selection, and Limitations

Summary of performance, parameter choices, and distinguishing features of principal ghost penalty bilinear forms is provided below:

Variant	Kernel	Locking as $\gamma\to\infty$	Condition Number
Classical face/bulk-ghost	Global polynomials	Locking	$\lesssim h^{-2}$
Aggregated FE (strong)	Aggregated space	Locking-free	$\lesssim h^{-2}$
Locking-free (distance-to-ext)	Aggregated space	Locking-free	$\lesssim h^{-2}$

Classical ghost penalties require careful tuning of the penalty parameter: small values yield inadequate stabilization, while large values cause locking. Locking-free forms—especially those penalizing extension distances or constructed via degree-of-freedom stabilization—exhibit robust error and conditioning characteristics over a wide range of penalty parameters.

In dynamic problems, the addition of ghost-mass terms is critical for ensuring uniform spectral bounds and stable, cut-insensitive time stepping. All stabilization approaches rely on mesh connectivity assumptions (bounded neighbor chains) and scaling conditions to guarantee their analytical and computational properties.

• "On the Design of Locking Free Ghost Penalty Stabilization and the Relation to CutFEM with Discrete Extension" (Burman et al., 2022)
• "Linking ghost penalty and aggregated unfitted methods" (Badia et al., 2021)
• "Critical time-step size analysis and mass scaling by ghost-penalty for immersogeometric explicit dynamics" (Stoter et al., 2023)
• "Matrix-Free Ghost Penalty Evaluation via Tensor Product Factorization" (Wichrowski, 28 Feb 2025)