Matrix-Free Ghost Penalty Evaluation via Tensor Product Factorization

TL;DR

The paper tackles ill-conditioning in CutFEM caused by small cut cells by introducing a matrix-free evaluation of the face-based ghost penalty. It exploits a tensor-product factorization of the ghost penalty operator, decomposing it into a sequence of one-dimensional matrix-vector products using precomputed 1D mass and penalty matrices, leading to a robust $O(k^{d+1})$ complexity for degree-$k$ elements in $d$ dimensions on Cartesian meshes. The approach avoids global matrix assembly, reduces memory footprint, and preserves high-order accuracy, with implementation in the deal.II library and validation through numerical experiments demonstrating optimal convergence and practical efficiency. This method enhances the scalability of high-order CutFEM in unfitted geometries by providing an efficient, stable stabilization mechanism compatible with matrix-free solvers.

Abstract

We present a matrix-free approach for implementing ghost penalty stabilization in Cut Finite Element Methods (CutFEM). While matrix-free methods for CutFEM have been developed, the efficient evaluation of high-order, face-based ghost penalties remains a significant challenge, which this work addresses. By exploiting the tensor-product structure of the ghost penalty operator, we reduce its evaluation to a series of one-dimensional matrix-vector products using precomputed 1D matrices, avoiding the need to evaluate high-order derivatives directly. This approach achieves $O(k^{d+1})$ complexity for elements of degree $k$ in $d$ dimensions, significantly reducing implementation effort while maintaining accuracy. The derivation relies on the fact that the cells are aligned with the coordinate axes. The method is implemented within the \texttt{deal.II} library.

Figures (7)

• Figure 1: Illustration of a computational domain discretization for CutFEM: The curved domain boundary $\partial\Omega_h$ (red dashed line) intersects a Cartesian mesh, creating cut cells (blue) and interior cells (green). Ghost penalty stabilization is applied across the faces $\mathcal{F}_h$ (thick blue lines) between cut cells or between cut and interior cells to ensure numerical stability regardless of the boundary position.
• Figure 2: Tensor product numbering of degrees of freedom for two adjacent quadratic elements $K_1$ and $K_2$ in 2D, sharing a face $F$. The numbers indicate the lexicographical ordering of the DoFs within each cell as well as the corresponding multi-indices. The figure shows a 2D case for simplicity of presentation.
• Figure 3: Relative time per application for different evaluation methods on a single cell. The time is normalized by the time of FEEvaluation for $k=1$. In 2D that is 0.077 $\mu$s, while in 3D it is 0.2218 $\mu$s.
• Figure 4: Problem with a single ball: Throughput of matrix-vector multiplication in degrees of freedom per second for different polynomial degrees. All throughput values are reported as DoFs/s, and only DoFs on active cells are included.
• Figure 5: Problem with multiple balls: Throughput of matrix-vector multiplication in DoFs per second for varying numbers of balls.