Ghost stabilisation for cut finite element exterior calculus
October 16, 2025
Abstract
We introduce the cut finite element method in the language of finite element exterior calculus, by formulating a stabilisation – for any form degree – that makes the method robust with respect to the position of the interface relative to the mesh. We
prove that the \(L^2\)-norm on the physical domain augmented with this stabilisation is uniformly equivalent to the \(L^2\)-norm on the “active” mesh that contains all the degrees of freedom
of the finite element space (including those external to the physical domain). We show how this CutFEEC method can be applied to discretize the Hodge Laplace equations on an unfitted mesh, in any dimension and any topology. A numerical illustration is
provided involving a conforming finite element space of \(H^{\text{curl}}\) posed on a filled torus, with convergence and condition number scaling independent of the position of the boundary with respect to the background
mesh.
MSC: 65N30, 14F40
Key words: Cut finite element method, FEEC, Hodge Laplace equation, Hilbert complex
1 Introduction↩︎
A successful strand of the literature on the numerical approximation of partial differential equations (PDEs) has emphasized the relevance of Hilbert complexes [1] in the design of stable numerical schemes. In this context, differential forms provide a natural unifying language to treat mixed problems set in arbitrary space dimension, possibly on manifolds, and where a combination of differential operators appear [2]–[7]. Finite Element Exterior Calculus (FEEC), particularly relevant for the present work, provides a comprehensive view of full and trimmed finite element approximations of the de Rham complex; see the monograph [8] for an introduction. At the same time, the Cut Finite Element Method (CutFEM) has emerged as a versatile technique for handling PDEs on domains with complex geometries, removing the need for the mesh to be compliant with the domain’s boundary; see the recent review [9] and references therein. Through weak imposition of boundary conditions and carefully constructed stabilisation terms over facets close to the boundary, CutFEM leads to stable and well-conditioned discretisations. Moreover, such properties hold independently of the boundary position relative to the unfitted mesh. This work aims at merging these research avenues by proposing a Cut Finite Element Exterior Calculus (CutFEEC) framework that adapts the CutFEM techniques to FEEC. This merger has notably been made possible by the recent extension of the CutFEM technology to mixed problems [10]–[12], which are a natural application of FEEC [8].
The prime example of a Hilbert complex is the de Rham complex, which, for an open connected domain \(\Omega \subset \mathbb{R}^3\), reads \[0\hookrightarrow H^1(\Omega) \xrightarrow{\nabla} \boldsymbol{H}^{\mathop{\mathrm{curl}}}(\Omega) \xrightarrow{\mathop{\mathrm{curl}}} \boldsymbol{H}^{\mathop{\mathrm{div}}}(\Omega) \xrightarrow{\mathop{\mathrm{div}}} L^2(\Omega) \to 0.\] Using the language of differential forms, the de Rham complex can be extended to a domain \(\Omega\) of any dimension \(n\in\mathbb{N}\) as: \[\label{eq:continuous46complex} 0\hookrightarrow H\Lambda^0\Omega \xrightarrow{d^0} H\Lambda^1\Omega \xrightarrow{d^1} \dots \xrightarrow{d^{n-1}} H\Lambda^n\Omega \to 0,\tag{1}\] where \(d^k\) is the exterior derivative and \(H\Lambda^k\Omega\) is the space of \(L^2\)-integrable differential \(k\)-forms on \(\Omega\) with \(L^2\)-integrable exterior derivative. The de Rham complex enters the well-posedness analysis of PDEs through its cohomology spaces \[\mathfrak{H}^k\coloneq \mathop{\mathrm{Ker}}d^k / \mathop{\mathrm{Im}}d^{k-1}.\] These spaces relate to the topology of the domain via their dimension. Preserving these homological structures at the discrete level is crucial for the stability of numerical schemes. A paradigmatic example of PDE problem whose well-posedness hinges on the de Rham complex is the Hodge Laplace equation; see, e.g., [8]. Its mixed formulation naturally leads to a mixed problem where the unknowns are the codifferential, the (scalar- or vector-) potential, and a Lagrange multiplier to enforce \(L^2\)-orthogonality with respect to harmonic forms.
The Finite Element (FE) approximation of mixed problems became a research topic starting from the late 1970s, when the first conforming approximation of vector-valued spaces in the de Rham complex were developed [13]–[15]. Traditional FE methods, such as the ones cited above, require meshes that conform to the geometry of the domain. These can be challenging to generate for intricate geometries or in the presence of evolving interfaces. CutFEM addresses this problem by permitting the interface to cut through the elements of a background mesh, thereby simplifying mesh generation and adaptation [16]. However, ensuring stability and accuracy in such unfitted methods necessitates appropriate stabilisation techniques [17]. Otherwise, when the geometry cuts the mesh in especially nasty ways (\(|T\cap\Omega|\ll |T|\) for an element \(T\) in the mesh), one observes a significant degradation of the condition number of the associated linear system. Recent works have addressed the application of CutFEM techniques to mixed problems. In [10], [12], the authors introduced a stable CutFEM for Darcy flow problems which, thanks to conformity with respect to the tail end of the de Rham complex, preserves the divergence of the velocity field. Using a similar approach, one can also preserve the divergence condition of Stokes flow in an unfitted setting [11].
In this work, we develop a new CutFEM stabilisation to design unfitted robust discrete \(L^2\)-products for any space in the de Rham complex of differential forms. The stabilisation can be seen as a generalisation of the mixed ghost penalty term introduced in [10], [11], [18]. Orthogonalising against the harmonic forms using an \(L^2\)-product including this stabilisation, we construct an arbitrary-order numerical method which is unfitted but still compatible with the continuous de Rham complex, thus ensuring its natural stability and preservation of relevant constraints (e.g., zero divergence of electric field in Maxwell’s equations). We develop facet-based stabilisation terms, but we mention that there exist other equivalent stabilisation terms in the CutFEM literature, though these have not been developed for differential forms, see [9]. We note that the commonly used assumption of a quasi-uniform mesh [9] is not needed in our construction or analysis. It is sufficient for the relevant arguments to use a local quasi-uniformity, which does not entail particular restrictions on the mesh. In [19], the quasi-uniform assumption is used but a similar remark is made [20].
A particularly relevant novelty results from applying this general construction to the special case of \(1\)-forms, which is key to discretize problems such as Maxwell’s equations or the \(H^{\mathop{\mathrm{curl}}}\)-elliptic problem. CutFEM has, indeed, seen little development for problems involving \(H^{\mathop{\mathrm{curl}}}\); see the recent and exhaustive review [9]. There are some previous works which use an unfitted discontinuous Galerkin approach [21], [22]. Otherwise, works which do not use discontinuous/broken finite element spaces are essentially limited to [[23]; [24]; [11]][18]. Of these, [11] and [18] do not contain a stability analysis; [11] considers a vorticity-velocity-pressure formulation for Stokes flow wherein the vorticity is discretized by a \(H^{\mathop{\mathrm{curl}}}\)-conforming variable, while [18] considers numerical experiments of a set of standard formulations of Maxwell’s equations made unfitted. The works [23], [24] focus, respectively, on the time-harmonic Maxwell’s equations and the \(H^{\mathop{\mathrm{curl}}}\)-elliptic interface problem. In [23], the authors mention that their proposed method works numerically when using nodal Lagrange elements, but error analysis is provided only for the case of using piecewise discontinuous polynomials. The numerical scheme in [23] uses the lowest order \(H^{\mathop{\mathrm{curl}}}\)-conforming Nédélec elements of the second kind, but adds diagonal stabilisation terms in the Lagrange multiplier block [24], thus perturbing the saddle point nature of the problem. Similar stabilisation terms also appear in [23]. Moreover, in contrast to our approach, no work mentioned above considers the presence of non-trivial harmonic forms, which appear when considering domains of general topology.
The remainder of this paper is organised as follows. In Section 2 we discuss a motivating example of mixed problem, the Hodge Laplace problem, whose well-posedness relies on the properties of the de Rham complex, and we summarise some notions of FEEC. In Section 3 we recall the notion of unfitted mesh, and we then introduce in Section 4 the concepts of tangential and normal parts of differential forms, relative to a facet of the mesh. In Section 5 we introduce the novel ghost penalty stabilisation operator, and show the uniform equivalence between the \(L^2\)-norm on the physical domain and the \(L^2\)-norm on the extended domain containing all the degrees of freedom. Finally, in Section 6 we validate our results on a numerical example on the filled torus, demonstrating that the use of the stabilisation drastically reduces the condition number of the discrete system.
2 A motivating example↩︎
In this section we briefly discuss the Hodge Laplacian, an archetypal problem whose well-posedness hinges on the properties of the de Rham complex. This example serves as a motivation for the development of a CutFEEC theory, and will be used in Section 6 to numerically demonstrate the efficiency of the ghost penalty stabilisation of Section 5.
2.1 Hodge Laplace equation↩︎
Let \(B\subset\mathbb{R}^n\) be an open set. We denote by \(\Lambda^k B\) the space of smooth differential \(k\)-forms on \(B\). We introduce the \(L^2\)-product \((\omega,\zeta)_B \mathrel{\vcenter{:}}= \int_B \omega\wedge\star \zeta\) on \(\Lambda^k B\), where \(\wedge\) denotes the wedge product and \(\star\) the Hodge star operator [8]. The \(L^2\)-norm over \(B\) is denoted by \(\|{\bullet}\|_B\mathrel{\vcenter{:}}= \sqrt{(\cdot,\cdot)_B}\). We denote by \(L^2\Lambda^k B\) the space of \(k\)-forms over \(B\) that are bounded with respect to this norm. With the exterior derivative \(d^k: L^2\Lambda^kB \to L^2\Lambda^{k+1}B\), viewed as an unbounded operator, we also define \[\begin{align} H\Lambda^k B &\mathrel{\vcenter{:}}= \left\{ \omega\in L^2\Lambda^k B: d^k\omega\in L^2\Lambda^{k+1} B\right\}, \end{align}\] with norm \[\begin{align} \|\omega\|^2_{HB} \mathrel{\vcenter{:}}= \|\omega\|^2_B+\|d^k\omega\|^2_B. \end{align}\]
Let now \(\Omega\subset\mathbb{R}^n\) be a connected and open subset of \(\mathbb{R}^n\) with Lipschitz boundary \(\partial \Omega\) and outward unit normal vector field \({\boldsymbol{n}}\). Given \(f\in L^2\Lambda^k\Omega\), the Hodge Laplace equation consists in seeking \(\eta \in \Lambda^k \Omega\) such that \[\label{eqs:HodgeLaplace} \begin{align} (d^{k-1}\delta^k + \delta^{k+1} d^k)\eta &= f-\pi f &\qquad& \text{in \Omega,} \\ \pi \eta &= 0 &\qquad& \text{in \Omega}, \end{align}\tag{2}\] where the codifferential \(\delta^{k+1}: L^2\Lambda^{k+1}\Omega \to L^2\Lambda^k\Omega\) is the adjoint of \(d\) and \(\pi\mathrel{\vcenter{:}}= \pi_{\mathfrak{H}^k}\) is the orthogonal projection onto the space of harmonic forms \[\begin{align} \mathfrak{H}^k \mathrel{\vcenter{:}}= \left\{ \rho \in L^2 \Lambda^k\Omega: d^k \rho = 0 \text{ and } \delta^k \rho=0 \right\} = \mathop{\mathrm{Ker}}d^k \cap \mathop{\mathrm{Ker}}\delta^k \cong \mathop{\mathrm{Ker}}d^k / \mathop{\mathrm{Im}}d^{k-1}. \end{align}\] Here \(\cong\) means vector space isomorphism. The boundary conditions are implicit in the definition of the Hodge Laplace operator \(d^{k-1}\delta^k + \delta^{k+1} d^k\), and are revealed upon choosing the domain of the exterior derivative. For ease of presentation, let us in this work choose the domain \[\begin{align} D(d^k) &= H\Lambda^k\Omega, \end{align}\] so that the domain of the codifferential becomes (by duality) the space of \(L^2(\Omega)\)-integrable forms \(\omega\) with \(\delta \omega \in L^2(\Omega)\) and \(\gamma(\star \omega) =0.\) After recalling vector proxies in three space dimensions in Table 1, in Table 2, we express problem 2 in terms of vector proxies for different values of \(k\); cf. [5] for further details.
| \(k\) | Form | Proxy | Sobolev space |
|---|---|---|---|
| 0 | \(q\) | \(q\) | \(H^1(\Omega)\) |
| 1 | \(v_1 dx^1 + v_2 dx^2 + v^3 dx^3\) | \(\boldsymbol{v} = (v_1,v_2,v_3)\) | \(\boldsymbol{H}^{\mathop{\mathrm{curl}}}(\Omega)\) |
| 2 | \(w_1 dx^2\wedge dx^3 - w_2 dx^1\wedge dx^3 + w_3 dx^1\wedge dx^2\) | \(\boldsymbol{w} = (w_1,w_2,w_3)\) | \(\boldsymbol{H}^{\mathop{\mathrm{div}}}(\Omega)\) |
| 3 | \(r dx^1\wedge dx^2\wedge dx^3\) | \(r\) | \(L^2(\Omega)\) |
| \(k=0\) | \(k=1\) | \(k=2\) | \(k=3\) |
|---|---|---|---|
| \(\nabla q\cdot {\boldsymbol{n}}= 0,\) | |||
| \(\bar q = 0\) | |||
| \(\boldsymbol{v} \perp \mathfrak{H}^1\) | |||
| \(\boldsymbol{w} \perp \mathfrak{H}^2\) | |||
| \(r = 0\) |
For convenience of notation, we shall for the remainder omit the index and just write \(d\) for the exterior derivative. A variational mixed formulation of 2 is (cf. [8]): Find \((\sigma,\eta,\lambda)\in H\Lambda^{k-1}\Omega \times H\Lambda^k\Omega \times \mathfrak{H}^k\) such that \[\label{eqs:mixed} \begin{align} (\sigma, \tau)_\Omega - (\eta, d\tau)_\Omega &= 0 \qquad &&\forall \tau\in H\Lambda^{k-1}\Omega, \\ (d \sigma, \zeta)_\Omega + (d\eta,d\zeta)_{\Omega} + (\lambda, \zeta)_\Omega &= (f, \zeta)_\Omega \qquad &&\forall \zeta\in H\Lambda^{k}\Omega, \\ (\eta, \rho)_\Omega &= 0 \qquad &&\forall \rho\in \mathfrak{H}^{k}. \end{align}\tag{3}\] The above formulation enforces in a weak manner \(\sigma = \delta \eta\) and \(\lambda = \pi f\). The (natural) boundary conditions are \(\gamma \star \eta = 0\) and \(\gamma \star d\eta = 0\).
Using crucial properties of the de Rham complex such as Hodge decompositions and Poincaré inequalities, it can be proved that this mixed formulation is well-posed [8]. Preserving these properties at the discrete level, as is the case for the compatible discretisations considered here, naturally leads to numerical schemes for which stability can be proved mimicking the arguments for well-posedness of the continuous problem.
2.2 Principles of FEEC↩︎
Assume here that \(\Omega\) is a polytopal domain, and let \(\mathcal{T}_h\) be a conforming simplicial mesh of \(\Omega\). The Finite Element Exterior Calculus (FEEC) method provides the tools for a discrete formulation of 3 which is well-posed [8], but non-conforming in general due to the non-inclusion of the discrete space of harmonic forms into the continuous one (these spaces have the same dimension, so the former being included in the latter would actually mean that one could compute harmonic forms exactly). We make a brief summary of the results which pertain to this work.
Let \(r\ge 1\) be a polynomial degree. For each \(k\in \{0,\dots,n\}\), let \(V_{h}^{k,r}\subset H\Lambda^k\Omega\) be a finite-dimensional space of \(k\)-forms that are piecewise polynomial of degree at most \(r\) on \(\mathcal{T}_h\) and such that \(dV_{h}^{k,r}\subset V_{h}^{k+1,r-1}\). We adopt the convention that \(V_h^{-1,r}\mathrel{\vcenter{:}}= \{0\}\) and \(V_h^{n+1,r}\mathrel{\vcenter{:}}= \{0\}\). These spaces define a sub-complex of the de Rham complex: \[\label{eq:discrete46complex} \begin{tikzcd} 0 \arrow[r,hook] & V_{h}^{0,r} \arrow{r}{d} & V_{h}^{1,r-1} \arrow{r}{d} & \cdots \arrow{r}{d} & V_{h}^{n,r-n} \arrow{r} & 0. \end{tikzcd}\tag{4}\] It can be shown [8], using, e.g., bounded cochain maps between the continuous and discrete de Rham complexes, that the cohomology spaces of 4 are isomorphic to those of 1 . This diagram (in which the polynomial degree decreases along the sequence) corresponds to the case of full finite element spaces, \(V_{h}^{k,r}\mathrel{\vcenter{:}}= P_r\Lambda^k\Omega_{h}\). When using trimmed spaces, \(P^-_r\Lambda^k\Omega_{h}\subset P_r\Lambda^k\Omega_{h}\) in the notation of [8], the degree \(r\) remains the same along the diagram. For application of our results to the trimmed spaces, simply replace any occurrence of \(V_{h}^{k,\bullet}\) with \(P^-_r\Lambda^k\Omega_{h}\).
The discrete exterior derivative acting on a discrete space is just the restriction of the exterior derivative from the corresponding continuous space. When relevant, we shall denote it by \(d_h\mathrel{\vcenter{:}}= d|_{V_{h}^{k,r}}\), otherwise by \(d\) when there is no confusion. Then, one defines the space of discrete harmonic forms as \[\label{eq:Hhk} \mathfrak{H}_h^k \mathrel{\vcenter{:}}= \mathop{\mathrm{Ker}}d_h / dV_{h}^{k-1,r+1} \cong \left\{\rho_h \in V_{h}^{k,r} :\; \text{d\rho_h=0 and (\rho_h,d\tau_h)_\Omega = 0 for all \tau_h\in V_{h}^{k-1,r+1}}\right\}.\tag{5}\] The non-inclusion \(\mathfrak{H}_h^k \not\subset \mathfrak{H}^k\) is due to forms in \(\mathfrak{H}^k\) needing to be orthogonal to all of \(dH\Lambda^{k-1}\Omega\), while being in \(\mathfrak{H}_h^k\) only guarantees orthogonality to the subspace \(dV_{h}^{k-1,r+1}\).
As mentioned in Section 2.1, the Hodge decomposition of the complex spaces is an essential tool to analyse both the continuous Hodge-Laplace equation and its numerical discretisations. For the complex 4 , we have the following discrete Hodge decomposition: \[\begin{align} \label{eq:discrete95hodge95decomp} V_{h}^{k,r} = (\mathop{\mathrm{Ker}}d_h)^\perp \mathbin{\raisebox{-1pt}{\begin{tikzpicture}[scale=0.134] \draw (0,-0.5)--(0,1); \draw (-0.866,-0.5)--(0.866,-0.5); \draw (0,0) circle [radius=1]; \end{tikzpicture}}}dV_{h}^{k-1,r+1} \mathbin{\raisebox{-1pt}{\begin{tikzpicture}[scale=0.134] \draw (0,-0.5)--(0,1); \draw (-0.866,-0.5)--(0.866,-0.5); \draw (0,0) circle [radius=1]; \end{tikzpicture}}}\mathfrak{H}_h^k. \end{align}\tag{6}\]
Given \(f\in L^2\Lambda^k\Omega\) the FEEC scheme for problem 3 reads: Find \((\sigma_h,\eta_h,\lambda_h)\in V_{h}^{k-1,r+1} \times V_{h}^{k,r} \times \mathfrak{H}_h^k\) such that \[\begin{align} {2} (\sigma_h, \tau_h)_\Omega - (\eta_h, d\tau_h)_\Omega &= 0 &&\qquad\forall \tau_h\in V_{h}^{k-1,r+1}, \\ (d \sigma_h, \zeta_h)_\Omega + (d\eta_h,d\zeta_h)_{\Omega} + (\lambda_h, \zeta_h)_\Omega &= (f, \zeta_h)_\Omega &&\qquad\forall \zeta_h\in V_{h}^{k,r}, \\ (\eta_h, \rho_h)_\Omega &= 0 &&\qquad\forall \rho_h\in \mathfrak{H}_h^{k}. \end{align}\] The scheme is well-posed and optimally convergent in the \(H\Lambda^k\Omega\)-norm [8].
3 Unfitted mesh↩︎
Let \(\mathcal{T}_h= \{ T\}\) be an active mesh such that each (open) simplex \(T \in \mathcal{T}_h\) has a nonempty intersection with \(\Omega\), and for which \[\Omega \subset \Omega_{h}\coloneq \left(\bigcup_{T\in\mathcal{T}_h} \overline{T}\right)^\circ.\] The open set \(\Omega_h\) is hereafter referred to as the active domain.
Remark 1 (Construction of active mesh). The term active mesh refers to the following construction. Define a background domain \(\Omega_0\supset \Omega\), which is polytopal of dimension \(n\). Consider a conforming triangular mesh \(\mathcal{T}_{0,h}\) of \(\Omega_0\). Then, the active mesh is defined as \(\mathcal{T}_h \mathrel{\vcenter{:}}= \left\{ T\in\mathcal{T}_{0,h} : T\cap \Omega \neq \emptyset \right\}.\) This construction is illustrated in Figure 1.
We say \(T\in\mathcal{T}_h\) is cut if \(T \not\subset \Omega\) and fully immersed in \(\Omega\) otherwise. We denote the set of cut elements by \[\mathcal{T}_h^{\text{cut}}\mathrel{\vcenter{:}}= \left\{ T\in\mathcal{T}_h : T \not\subset \Omega \right\}.\] We denote by \(\mathcal{F}^\partial_h\) the set of stabilisation facets [11], which are facets internal to \(\Omega_h\) for which at least one neighbouring element is intersected by \(\partial\Omega\): \[\mathcal{F}^\partial_h\mathrel{\vcenter{:}}= \left\{ F \in \bigcup_{T\in\mathcal{T}_h^{\text{cut}}}\mathcal{F}_T \;:\; F \not\subset \partial \Omega_h \right\},\] where \(\mathcal{F}_T\) gathers the (\((n-1)\)-dimensional) facets of \(T \in \mathcal{T}_h\). To each facet we associate a unique normal vector \({\boldsymbol{n}}\). The context will make it clear which normal is meant, that of a facet or that of \(\partial\Omega\).
Assumption 1 (Mesh regularity). For any \(T \in \mathcal{T}_h\), let \(\varrho_T\) and \(h_T\) be, respectively, the diameter of the largest ball contained in \(T\) and that of \(T\). We assume the following:
(Shape regularity) For some constant \(\kappa>0\) independent of the mesh size \(h \mathrel{\vcenter{:}}= \max_{T \in \mathcal{T}_h} h_T\), we have \[\label{eq:shape46regularity} \frac{h_T}{\varrho_T} \leq \kappa \qquad \forall T\in\mathcal{T}_h.\qquad{(1)}\]
(Uniformly bounded cut-to-uncut path) There exists an integer \(N > 0\) independent of \(h\) such that, for all \(T \in \mathcal{T}_h^{\text{cut}}\), there exists a sequence of \(N_T \le N\) elements \(\{T_1 = T, T_2,\ldots,T_{N_T}\}\) such that \(T_i\) and \(T_{i+1}\) share a facet for all \(1\le i\le N_T - 1\), and \(T_{N_T}\in\mathcal{T}_h\setminus \mathcal{T}_h^{\text{cut}}\) is fully immersed in \(\Omega\).
The notation \(a\eqsim b\) means that there exist constants \(C,C'>0\) independent of \(h\), but possibly dependent on \(\kappa\) and \(N\) such that \(a\leq Cb\) and \(b\leq C'a\), and similarly we write \(a\lesssim b\) if and only if \(a\leq Cb\).
Remark 2 (Quasi-uniformity along cut-to-uncut paths). We remark that, by shape regularity and since \(N\) is bounded independently of \(h\), all mesh elements along a given cut-to-uncut path have uniformly comparable diameters: if \(T_1,\ldots,T_{N_T}\) are as in [item:cut-to-uncut], then \(h_{T_i}\eqsim h_{T_j}\) for all \(i,j\in\{1,\ldots,N_T\}\).
4 Tangential and normal traces and parts of differential forms↩︎
Restrictions as well as tangential and normal traces and parts of differential forms are essential concepts to define the stabilisation in 5 below. We consider here an element \(T\in\mathcal{T}_h\) and one of its facets \(F\subset\partial T\).
For \(x\in F\), let \(\mathbb{T}_x C\) be the tangent space of \(C \in \{ F, T\}\) at \(x\). Let \[\iota:F\hookrightarrow T\] be the inclusion map of \(F\) into \(T\). We note that \(\mathbb{T}_xF\) is the vector hyperplane parallel to \(F\) in \(\mathbb{R}^n\), and that \(X\in \mathbb{T}_xF\) can be considered as a vector in \(\mathbb{R}^n\); for this reason, we identify \(\iota_*X\) with \(X\).
For \(\omega\in \Lambda^kT\), by “restriction to the facet \(F\)” we simply mean the restriction of each component of \(\omega\) to \(F\), i.e., writing \(\omega \mathrel{\vcenter{:}}= \sum_{1\leq j_1<j_2< \dots <j_k \leq n} \omega_{j_1 j_2 \dots j_k} \, dx^{j_1}\wedge dx^{j_2}\wedge \dots \wedge dx^{j_k}\), we define \[\label{eq:def46restriction} \omega|_F \mathrel{\vcenter{:}}= \sum_{1\leq j_1<j_2< \dots <j_k \leq n} (\omega_{j_1 j_2 \dots j_k}|_F) \, dx^{j_1}\wedge dx^{j_2}\wedge \dots \wedge dx^{j_k}.\tag{7}\] Notice that the result is a form which formally acts on \(k\)-tuples of vectors in \(T\), although its components will only be evaluated on \(F\). The reason is that we want \(\omega|_F\) as a differential form to also be able to act on vectors normal to \(F\).
The standard trace operator \(\gamma \mathrel{\vcenter{:}}= \iota^* : \Lambda^kT \to \Lambda^kF\) on the facet is defined as the pullback of the inclusion \(\iota\): For all \(\omega \in \Lambda^k T\),\[\label{eqs:trace} \forall x\in F,\qquad (\gamma \omega)_x(X_1,\dots ,X_k) \mathrel{\vcenter{:}}= \omega_x(X_1,\dots,X_k) \qquad \forall X_1,\dots,X_k\in \mathbb{T}_xF.\tag{8}\]
Recall the Hodge star operator \(\star:\Lambda^kT\to \Lambda^{n-k}T\). We can also define a Hodge star \(\star_F:\Lambda^kF\to \Lambda^{n-1-k}F\) on \(F\) satisfying \(\star_F\star_F \omega=(-1)^{k(n-1-k)}\omega\) for all \(\omega\in \Lambda^kF\).
The normal trace operator is defined by “sandwiching” the standard trace between the Hodge star operators on \(T\) and on \(F\), see [25]: \(\gamma_{\boldsymbol{n}}\mathrel{\vcenter{:}}= (-1)^{n(k-1)} \star_F\gamma\star: \Lambda^kT\to \Lambda^{k-1}F\) is therefore such that, for all \(\omega \in \Lambda^k T\), \[\label{eqs:normaltrace} \forall x\in F,\qquad (\gamma_{\boldsymbol{n}}\omega)_x(X_2,\dots ,X_k) \mathrel{\vcenter{:}}= \omega_x({\boldsymbol{n}},X_2,\dots,X_k) \qquad \forall X_2,\dots,X_k\in \mathbb{T}_xF.\tag{9}\] The latter expression is, by definition, the contraction \({\boldsymbol{n}}\neg \omega\) of \(\omega\). The following identities hold, [25]: \[\begin{align} \label{eqs:trace95relations} \star_F \gamma_{\boldsymbol{n}}= \gamma \star, \qquad \gamma_{\boldsymbol{n}}\star = (-1)^k\star_F\gamma. \end{align}\tag{10}\] From now on, we shall denote both \(\star\) and \(\star_F\) by \(\star\), the context making it clear which one is meant. The vector proxy interpretation of the trace and normal trace operators are summarized in Table 3.
| Form degree \(k\) | Proxy | Trace (\(\gamma\)) | Normal trace (\(\gamma_{\boldsymbol{n}}\)) | |
|---|---|---|---|---|
| 0 | \(q\) | \(q\) | 0 | |
| 1 | \(\boldsymbol{v}\) | \({\boldsymbol{n}}\times (\boldsymbol{v} \times {\boldsymbol{n}})\) | \(\boldsymbol{v} \cdot {\boldsymbol{n}}\) | |
| 2 | \(\boldsymbol{w}\) | \(\boldsymbol{w} \cdot {\boldsymbol{n}}\) | \(- {\boldsymbol{n}}\times \boldsymbol{w}\) | |
| 3 | \(r\) | \(0\) | \(r\) |
Sometimes, the standard trace is called the tangential trace, in contrast to the normal trace. This nomenclature can be understood comparing the traces with the restriction for sufficiently smooth differential forms. Consider for the remainder of this section a smooth \(k\)-form \(\omega\in\Lambda^kT\). Using the pullback of the orthogonal projection \(\pi_F: T\to F\), we can define the tangential part of \(\omega\) on \(F\) by setting \[\label{eq:beta46parallel} \omega_\parallel \mathrel{\vcenter{:}}= (\pi_F^*\circ \gamma)\omega.\tag{11}\] At any \(x \in F\), it is also defined as an alternating form in \(\mathop{\mathrm{Alt}}^k \mathbb{R}^n\). Setting then \[\label{eq:beta46perp} \omega_\perp\mathrel{\vcenter{:}}=\omega|_F-\omega_\parallel,\tag{12}\] we have the following decomposition: \[\begin{align} \label{eq:parallel95decomp} \omega|_F = \omega_\parallel + \omega_\perp. \end{align}\tag{13}\]
Lemma 1 (Characterisation of the tangential and normal parts). Let \(\omega\in \Lambda^kT\) be a smooth \(k\)-form. Then the following holds: \[\begin{align} \label{eq:beta46parallel610606162gamma46beta610} \omega_\parallel = 0 &\iff \gamma\omega = 0, \\ \label{eq:beta46perp610606162gamma46n46beta610} \omega_\perp = 0 &\iff \gamma_{\boldsymbol{n}}\omega = 0. \end{align}\] {#eq: sublabel=eq:eq:beta46parallel610606162gamma46beta610,eq:eq:beta46perp610606162gamma46n46beta610} Moreover, \[\label{eq:star46traces} \text{ (\star\omega)_\parallel = \star\omega_\perp \quad and \quad (\star\omega)_\perp = \star\omega_\parallel. }\qquad{(2)}\]
Proof. 1) Proof of ?? . We start by showing that \(\omega_\parallel=0\) if and only if \(\gamma \omega = 0\). It is clear by its definition 11 that \(\omega_\parallel=0\) if \(\gamma\omega = 0\). To prove the converse, we just note that, if \((\omega_\parallel)_x=0\), then, for all \(X_1,\dots,X_k\in \mathbb{T}_xF\subset \mathbb{R}^n\), we have \[(\gamma\omega)_x(X_1,\dots,X_k)=(\gamma\omega)_{\pi_F(x)}(D\pi_F(X_1),\dots,D\pi_F(X_k))= (\omega_\parallel)_x(X_1,\dots,X_k)=0,\] the first equality following from the fact that \(D\pi_F\) is
the orthogonal projection on \(\mathbb{T}_xF\).
2) Proof of ?? . Let \(Y_1,\dots, Y_k\in \mathbb{R}^n\). For an appropriately chosen coordinate system with respect to \(F\), each vector \(Y_j\) can be written as \(Y_j=X_j+Y_j^{\boldsymbol{n}}\), where \(X_j=D\pi_F Y_j\in \mathbb{T}_xF\) and \(Y_j^{\boldsymbol{n}}\)
is the normal component of \(Y_j\). We then have \[(\omega_\perp)_x(Y_1,\dots,Y_k) \begin{align}[t] \overset{\eqref{eq:beta46perp},\,\eqref{eq:beta46parallel}}&= \omega_x(Y_1,\dots,Y_k) -
(\pi_F^*\circ\gamma)\omega_x(Y_1,\dots,Y_k) \\ \overset{\eqref{eqs:trace}}&= \omega_x(Y_1,\dots, Y_k) - \omega_x(X_1,\dots, X_k) \\ &= \left[ \omega_x(Y^{\boldsymbol{n}}_1,X_2,\dots,X_k)+\dots+\omega_x(X_1,\dots, X_{k-1},Y^{\boldsymbol{n}}_k) +
\omega_x(X_1,\dots, X_k) \right] - \omega_x(X_1,\dots,X_k) \\ &= \omega_x(Y^{\boldsymbol{n}}_1,X_2,\dots,X_k)+\dots+\omega_x(X_1,\dots, X_{k-1},Y^{\boldsymbol{n}}_k), \end{align}\] where, in the third step, we have used multiple times the fact
that, for any \(x \in F\), it holds \(\omega_x(X_1,\ldots,Y_j^{\boldsymbol{n}},\ldots,Y_k) = \omega_x(X_1,\ldots,Y_j^{\boldsymbol{n}},\ldots,X_k)\) since \(\omega_x\) is \(k\)-linear and alternating. The above relation shows that there exists permutations \(\sigma_j\), \(1\le j \le
k\), of the indices \(1,\dots,k\) such that, taking \(c_j\) such that \(Y_j^{\boldsymbol{n}}=c_j{\boldsymbol{n}}\), \[(\omega_\perp)_x(Y_1,\dots,Y_k) = \sum_{j=1}^k c_{\sigma_j(k)}(\gamma_{\boldsymbol{n}}\omega)_x(X_{\sigma_j(1)},\dots,X_{\sigma_j(k-1)}).\] This relation shows that \(\gamma_{\boldsymbol{n}}\omega =
0\) implies \(\omega_\perp=0\). To prove the converse statement, we observe that, for any \(x \in F\) and any vectors \(X_2,\ldots,X_k \in \mathbb{T}_x
F\), since \(D\pi_{\mathbb{T}_xF}{\boldsymbol{n}}=0\), we have \((\omega_\parallel)_x({\boldsymbol{n}},X_2,\dots,X_k)=0\) and thus, by 12 , \((\omega_\perp)_x({\boldsymbol{n}},X_2,\ldots,X_k)=\omega_x({\boldsymbol{n}},X_2,\ldots,X_k)=(\gamma_{\boldsymbol{n}}\omega_x)(X_2,\dots,X_{k})\). This shows that \(\omega_\perp = 0\) implies
\(\gamma_{\boldsymbol{n}}\omega = 0\). By 10 , ?? also implies that \(\gamma(\star\omega)=0\) if and only if \(\omega_\perp=0\), as is stated in [4].
3) Proof of ?? . Note that \(\star\omega =\star\omega_\parallel+\star\omega_\perp = (\star\omega)_\parallel + (\star\omega)_\perp\). We have that \[\begin{align}
(\star\omega)_\parallel = \star\omega_\perp, \quad (\star\omega)_\perp = \star\omega_\parallel,
\end{align}\] by applying vectors in \(\mathbb{R}^n\), respectively tangential and normal to \(F\), to the equation \(\star\omega_\parallel+\star\omega_\perp
- (\star\omega)_\parallel - (\star\omega)_\perp=0\). ◻
It can be helpful to look at the decomposition of a \(k\)-form \(\omega\in\Lambda^kT\) in a coordinate system adapted to the facet \(F\). Fix any \(x\in F\). Let \(\{{\boldsymbol{n}},e_1,\dots,e_{n-1}\}\) be a basis of \(\mathbb{R}^n\) such that \(\{e_1,\dots,e_{n-1}\}\) is a basis of \(\mathbb{T}_xF\). Let the dual basis be \(\{ dt,dx^1,\dots, dx^{n-1}\}\), where \(dt\) is the dual to \({\boldsymbol{n}}\). Define the multi-index set \[\mathcal{J}_\ell \mathrel{\vcenter{:}}= \left\{ (i_1,\dots, i_{\ell}) : 1\leq i_1 < i_2 < \dots < i_{\ell} \leq n-1 \right\}.\] Then, for any \(x\in F\), we can write \[\begin{align} \label{eq:normal95tangential95decomp95coord} (\omega)_x = (\omega_\perp)_x + (\omega_\parallel)_x = \sum_{J\in \mathcal{J}_{k-1}} \omega_{t,J}(x) dt\wedge dx^J + \sum_{I\in \mathcal{J}_{k}}\omega_I(x) dx^I, \end{align}\tag{14}\] where \(\omega_{t,J}\) and \(\omega_I\) are coefficients with respect to the chosen basis. One can check (upon identifying \(\gamma dx^{i}=dx^i\)) that \[\label{eq:normal95trace95coord} \text{ (\gamma\omega)_x = \sum_{I\in \mathcal{J}_{k}}\omega_I(x) dx^I\quad and \quad (\gamma_{\boldsymbol{n}}\omega)_x = \sum_{J\in \mathcal{J}_{k-1}} \omega_{t,J}(x) dx^J. }\tag{15}\] These relations shed more light on the relations stated in 1. They are also useful for the next result.
Lemma 2. Let \(\omega\in\Lambda^kT\) and \(\eta\in\Lambda^{n-k}T\). The following relations hold at any \(x\in F\): \[\omega_\parallel \wedge \eta_\parallel = 0,\qquad \omega_\perp \wedge \eta_\perp = 0.\]
Proof. Writing the decomposition 14 for both \(\omega\) and \(\eta\), it is clear that \(\omega_\perp\wedge\eta_\perp\) vanishes due to \(dt\) being present in both expansions. On the other hand, we have \[\begin{align} \omega_\parallel \wedge \eta_\parallel &= \sum_{I\in\mathcal{J}_k,\;L\in\mathcal{J}_{n-k}} \omega_I\eta_L dx^I\wedge dx^L, \end{align}\] which vanishes since one index at least is common to both \(I\) and \(L\), as can be seen by writing (identifying the indices in the families with the sets of these indices): \[|I\cap L| = |I|+|L|-|I\cup L| = k+(n-k)-|I\cup L| \geq n-(n-1)=1.\qedhere\] ◻
5 Ghost penalty↩︎
Using, in FEEC discretisations of weak formulations of PDEs, the inner product \((\bullet,\bullet)_{\Omega_{h}}\) on the active mesh \(\Omega_{h}\) in the discrete formulation is possible, but introduces a geometric error that is at least \(O(h)\) due to the fact that the mesh is not fitted to the boundary \(\partial\Omega\). Avoiding this geometric error requires to only rely on the inner product \((\bullet,\bullet)_\Omega\) in these formulations, which leads to ill-conditioned systems since DOFs attached to elements having a very small intersection with \(\Omega\) have a very small impact on this inner product. The goal of the ghost penalty stabilisation term designed in this section is precisely to restore robustness while preserving optimal convergence rates. The name “ghost” derives from the finite difference terminology of calling “ghost cells” the elements outside of the physical domain. The main result of this section is 1, where we prove that the norm associated to the ghost term is equivalent to the \(L^2(\Omega_{h})\)-norm.
5.1 Jumps of forms↩︎
Let \(F=\overline{T}_1 \cap \overline{T}_2\in \mathcal{F}^\partial_h\) with normal \({\boldsymbol{n}}\) pointing from \(T_2\) into \(T_1\). Set, for \(\omega\in V_{h}^{k,r}\) and \(i=1,2\), \[\omega_i \mathrel{\vcenter{:}}= \omega|_{T_i},\] where \(\omega|_{T_i}\) denotes the standard component-wise restriction of \(\omega\) to \(T_i\) on the canonical basis of \(k\)-alternating multi-linear maps, similar to 7 . We define the jump over \(F\) as \[\begin{align} [\cdot] : \Lambda^k(\overline{T}_1\cup\overline{T}_2) &\to C^\infty(F;\;\mathop{\mathrm{Alt}}^k\mathbb{R}^n),\\ \omega &\mapsto [\omega] \mathrel{\vcenter{:}}= \omega_1|_F- \omega_2|_F. \end{align}\]
Remark 3 (Traces and jumps). Since the normal and tangential traces 8 –9 depend only on the values of the form in \(\mathop{\mathrm{Alt}}^k\mathbb{R}^n\) on the facet \(F\), we can apply these operators to the jump defined above. In other words, \(\gamma[ \omega]\) and \(\gamma_{\boldsymbol{n}}[ \omega]\) make sense.
From the linearity of the jump operator and the relation between the tangential part and the trace operator (cf. 1), it follows that only the normal part of the form (see 4) is involved in the jump.
Proposition 1 (Jump of discrete forms). For all \(\omega\in V_{h}^{k,r}\), it holds \[[\omega] = [\omega_{\perp}],\] where \(\omega_\perp\) is the normal part of the form defined by 12 .
Proof. Let \(F=\overline{T}_1\cap \overline{T}_2\). We have that \(\omega_i \overset{\eqref{eq:parallel95decomp}}=\omega_{i\parallel}+\omega_{i\perp}\) on \(F\) for \(i=1,2\). The discrete \(k\)-forms are exactly defined by their degrees of freedom, so that \(\gamma \omega\) is continuous across inter-element facets, i.e. \(\gamma \omega_1-\gamma \omega_2=0\). Let \(\{Y_1,Y_2,\dots, Y_k\}\) be an arbitrary set of \(k\) vectors in \(\mathbb{R}^n\). By definition 11 of the tangential part of \(\omega\), it holds \[(\omega_{i\parallel})_x(Y_1,\dots,Y_k) = (\gamma \omega_i)_x(D\pi_F Y_1,\dots,D\pi_F Y_k) \qquad \forall x\in F,\] where \(\pi_F\) is the orthogonal projection onto \(F\). By linearity of the jump operator, we have: \[(\omega_{1\parallel}-\omega_{2\parallel})_x(Y_1,\dots,Y_k) = (\gamma \omega_1-\gamma \omega_2)_x(D\pi_F Y_1,\dots,D\pi_F Y_k)=0 \qquad \forall x\in F. \qedhere\] ◻
The result of Proposition 1 is translated in terms of vector proxies in 3.
| \(k\) | Proxy | Jump \([\omega_\perp]\) | |
|---|---|---|---|
| 0 | \(q\) | \(0\) | |
| 1 | \(\boldsymbol{v}\) | \([\boldsymbol{v} \cdot {\boldsymbol{n}}]\) | |
| 2 | \(\boldsymbol{w}\) | \([{\boldsymbol{n}}\times (\boldsymbol{w} \times {\boldsymbol{n}})]\) | |
| 3 | \(r\) | \([r]\) |
5.2 Ghost penalty stabilisation↩︎
For an integer \(\ell\geq 0\), we define the \(\ell\)-th order directional derivative \(\nabla_{{\boldsymbol{n}}}^{(\ell)}: V_{h}^{k,r}\to L^2 \Lambda^k \Omega_h\) along \({\boldsymbol{n}}\) as \[\nabla^{(\ell)}_{{\boldsymbol{n}}} \omega \mathrel{\vcenter{:}}= \sum_{1\leq j_1<j_2< \dots <j_k \leq n} (\partial^{(\ell)}_{{\boldsymbol{n}}} \omega_{j_1 j_2 \dots j_k}) \, dx^{j_1}\wedge dx^{j_2}\wedge \dots \wedge dx^{j_k},\] where \(\partial^{(\ell)}_{{\boldsymbol{n}}} \omega_{j_1 j_2 \dots j_k}\) is the \(\ell\)-th order directional derivatives of its component with respect to a chosen coordinate system. The definition does not depend on the choice of coordinate system, and we also identify \(\nabla_{\boldsymbol{n}}^{0}\equiv \mathop{\mathrm{Id}}.\)
For \(\omega\in V_{h}^{k,r}\), it is clear that \(\nabla^{(\ell)}_{\boldsymbol{n}}\omega\) is a polynomial form of degree \(r-\ell\) and that \(\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega=0\) for \(\ell>r\). Notice that the jump of the normal derivatives of a discrete \(k\)-form \(\omega\) does not reduce, in general, to the jump of the normal part \([(\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega)_\perp]\), in contrast to 1.
Recall the definitions 8 of trace and 9 of normal trace, as well as Remark 3 on their applicability to jumps. Letting \(h_F\) be the diameter of \(F\), we can then define the (ghost penalty) stabilisation term \(s : V_{h}^{k,r} \times V_{h}^{k,r} \to \mathbb{R}\) as \[\label{eq:ghost95penalty} s(\omega,\zeta) \mathrel{\vcenter{:}}= \sum_{F\in\mathcal{F}^\partial_h} \sum_{\ell=0}^r \eta h_F^{2\ell+1} \int_{F} \left( \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\zeta] + \gamma[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star\gamma [\nabla_{{\boldsymbol{n}}}^{(\ell)}\zeta] \right),\tag{16}\] where \(\eta>0\) is a user-defined penalty parameter. The form \(s\) is symmetric and bilinear. Notice that \(s\) can be extended to act on smooth enough \(k\)-forms and it holds, owing to the single-valuedness of traces, \[s(\zeta,\omega) = s(\omega,\zeta) = 0 \qquad \forall (\zeta,\omega) \in H^{r+1}\Lambda^k \Omega_h \times V_{h}^{k,r}.\] We define the ghost inner product and norm respectively as \[\label{eqs:def95stab95inner95product} \begin{alignedat}{2} (\omega,\zeta)_s &\mathrel{\vcenter{:}}= (\omega,\zeta)_{\Omega} + s(\omega,\zeta) &\qquad& \forall (\omega,\zeta) \in V_{h}^{k,r}\times V_{h}^{k,r}, \\ \|\omega\|_s &\mathrel{\vcenter{:}}= \sqrt{(\omega,\omega)_s} &\qquad& \forall \omega\in V_{h}^{k,r}. \end{alignedat}\tag{17}\]
Remark 4 (Formula 16 for vector proxies). Let \(n=3\). The simplifications of the jumps for the case \(\ell = 0\) given in 4 do not extend to the case \(\ell \ge 1\). However, the directional derivative of a \(k\)-form is still a \(k\)-form, so, using the results of 3, we infer that the integrand in 16 is in fact the inner product of the full jumps of the normal derivatives regardless of form order \(k\), see 5.
| Form degree \(k\) | ||
|---|---|---|
| 0 or 3 | \(\displaystyle\sum_{F\in\mathcal{F}^\partial_h}\sum_{\ell=0}^r \eta h_F^{2\ell+1} \int_F [\partial^{\ell}_{\boldsymbol{n}}\bullet][\partial^{\ell}_{\boldsymbol{n}}\bullet]\) | |
| 1 or 2 | \(\displaystyle\sum_{F\in\mathcal{F}^\partial_h}\sum_{\ell=0}^r \eta h_F^{2\ell+1} \int_F [\partial^{\ell}_{\boldsymbol{n}}\bullet ]\cdot[\partial^{\ell}_{\boldsymbol{n}}\bullet ]\) |
The goal of the remainder of this subsection is to show the following equivalence result.
Theorem 1 (Uniform equivalence for \(L^2\)-norms). The norms \(\|\bullet\|_{\Omega_{h}}\) and \(\|\bullet\|_s\) are equivalent on \(V_{h}^{k,r}\), uniformly in \(h\).
The following technical lemma is what informs the definition of the stabilisation term 16 .
Lemma 3 (Local control through ghost penalty). Let \(r \ge 1\) be an integer and \(\omega \in V_{h}^{k,r}\). Fix a boundary facet \(F\in \mathcal{F}^\partial_h\) shared by the elements \(T_1\) and \(T_2\) of \(\mathcal{T}_h\). Then, the following inequality holds: \[\begin{align} \|\omega\|^2_{T_1} &\lesssim \|\omega\|^2_{T_2} + \sum_{\ell=0}^r h_F^{2\ell+1} \int_{F}\left(\gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] + \gamma[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star\gamma [\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \right). \label{eq:fromT1toT2} \end{align}\qquad{(3)}\]
Proof. Consider the orthogonal projection onto the hyperplane containing \(F\), which we also denote by \(\pi_F\). First, we claim that we can assume to start with a simplex \(T_1\) in which the vertex \(\mathrm{s}_F\) opposite to \(F\) is such that \(\pi_F \mathrm{s}_F \in F\). Indeed, if \(T_1\) is such that \(\pi_F \mathrm{s}_F\not\in F\), then consider the maximal simplex \(T_1'\subseteq T_1\) with facet \(F\) and opposite vertex \(\mathrm{s}_F' \in T_1\) such that \(\pi_F \mathrm{s}_F' \in F\). By a standard scaling argument as in [26], using norm equivalence for polynomial forms, we have that \(\|\omega\|_{T_1}\lesssim \|\omega\|_{T_1'}\).
Consider now a point \(x\in T_1\), and its orthogonal projection \(x_F\) on \(F\). Let \(t\in\mathbb{R}\) be the distance from \(x\) to \(x_F\), such that \(x=x_F+t{\boldsymbol{n}}\). For \(i\in\{1,2\},\) since the restriction \(w_i\mathrel{\vcenter{:}}= w|_{T_i}\) is a polynomial form, there is a canonical extension outside of \(T_i\). The Taylor expansion centered at \(x_F\) of \(\omega_1\), or of the extension of \(\omega_2\) to \(T_1\), is given by \[\begin{align} (\omega_i)_x &= \sum_{\ell=0}^r \frac{(\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega_i)_{x_F}}{\ell!} t^\ell. \end{align}\] Since each \(\omega_i\) is a polynomial form of degree \(r\), the Taylor expansion is exact. Taking the difference of the two expansions gives \[(\omega_1)_x = (\omega_2)_x + \sum_{\ell=0}^r \frac{(\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega_1)_{x_F}-(\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega_2)_{x_F}}{\ell!} t^\ell = (\omega_2)_x + \sum_{\ell=0}^r \frac{[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]_{x_F}}{\ell!} t^\ell.\] Taking the squared \(L^2(T_1)\)-norm and using triangle inequalities, we get \[\begin{align} \|\omega_1\|^2_{T_1} &\lesssim \left(\|\omega_2\|_{T_1} + \frac{|t|^{\ell}}{\ell!} \sum_{\ell=0}^r \|[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]_{x_F}\|_{T_1}\right)^2 \\ &\lesssim \|\omega_2\|^2_{T_1} + h_F^{2\ell} \sum_{\ell=0}^r \|[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]_{x_F}\|^2_{T_1} \\ &\lesssim \|\omega_2\|^2_{T_2} + h_F^{2\ell} \sum_{\ell=0}^r \int_{T_1}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]_{x_F}\wedge \star [\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]_{x_F} , \end{align}\] where, in the second inequality, we have used the Cauchy–Schwarz inequality on the norm together with \(|t|\lesssim h_T\eqsim h_F\), and the conclusion follows from a polynomial transport argument (similar to the one used in the proof of [26]) to replace \(\|\omega_2\|^2_{T_1}\) with \(\|\omega_2\|^2_{T_2}\), and the definition of the \(L^2\)-norm to expand the second term.
Let \(\xi\) be the \(k\)-form on \(T_1\) such that, for all \(x \in T_1\), \(\xi_x\coloneq [\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]_{x_F}\); as above, \(x_F\) denotes the projection of \(x\) onto \(F\). From 13 we can decompose \(\xi\) into its normal and tangential parts inside \(T_1\) as \[\xi = \xi_\perp + \xi_\parallel.\] By definition, \(\xi_x=\xi_{x_F}\) for all \(x\in T_1\). Define on \(T_1\) the \((n-k)\)-form \(\eta\coloneq\star\xi\), for which we also have \(\eta_x=\eta_{x_F}\) for all \(x \in T_1\). By Lemma 1, we have \(\star \xi_\parallel =\eta_\perp\) and \(\star\xi_\perp=\eta_\parallel\) which, combined with 2, gives \[\text{ \xi_\perp \wedge \star \xi_\parallel = \xi_\perp \wedge \eta_\perp = 0 and \xi_\parallel \wedge \star \xi_\perp = \xi_\parallel \wedge \eta_\parallel = 0. }\] Hence, at \(x_F\) we can can decompose the wedge product as \[\xi \wedge \eta = \xi_\perp \wedge \eta_\parallel + \xi_\parallel\wedge\eta_\perp.\] Decomposing \(\xi\) using parallel and normal coordinates to \(F\) as in 14 , we write \[\xi_x = (\xi_\perp)_x + (\xi_\parallel)_x = \sum_{J\in\mathcal{J}_{k-1}} \xi_{t,J}(x) dt\wedge dx^J + \sum_{I\in \mathcal{J}_k}\xi_I(x) dx^I,\] and similarly for \(\eta\). By Fubini’s theorem, denoting by \(\chi\) the characteristic function of \(T_1\) and using a coordinate system \((t,y)\) with \(t\) the coordinate along \({\boldsymbol{n}}\) and \(y\) the coordinate along the affine span of \(F\), we get \[\begin{align} \int_{T_1} \xi_\perp \wedge \eta_\parallel &= \int_{\mathbb{R}^{n-1}}\int_{\mathbb{R}} \chi(t,x_F) \sum_{J\in\mathcal{J}_{k-1},I\in\mathcal{J}_k}\xi_{t,J}(x_F)\eta_I(x_F) \;dt\wedge dx^J\wedge dx^I \\ &= \int_{\mathbb{R}^{n-1}} \sum_{J\in\mathcal{J}_{k-1},I\in\mathcal{J}_k} \xi_{t,J}(x_F)\eta_I(x_F) \left( \int_{\mathbb{R}} \chi(t,x_F) dt \right) dx^J\wedge dx^I \\ &= \int_F L_{x_F}\sum_{J\in\mathcal{J}_{k-1},I\in\mathcal{J}_k} \xi_{t,J}(x_F)\eta_I(x_F) dx^J\wedge dx^I \\ &= \int_F L_{x_F} \sum_{J\in\mathcal{J}_{k-1}} \xi_{t,J} dx^J \wedge \sum_{I\in\mathcal{J}_k} \eta_I dx^I \\ \overset{\eqref{eq:normal95trace95coord}}&= \int_F L_{x_F} \gamma_{\boldsymbol{n}}\xi \wedge \gamma\eta, \end{align}\] where, starting from the third equality, \(L_{x_F}\) denotes the length of the segment at the vertical of \(x_F\) in \(T_1\). Performing similar computations for \(\int_{T_1} \xi_\parallel\wedge\eta_\perp\) and accounting for the fact that \(dx^I\wedge dt= (-1)^{k} dt\wedge dx^I\), we arrive at \[\begin{align} \int_{T_1} \xi \wedge \eta &= \int_F L_{x_F}\;\gamma_{\boldsymbol{n}}\xi \wedge \gamma\eta + (-1)^{k}\int_F L_{x_F}\;\gamma\xi \wedge \gamma_{\boldsymbol{n}}\eta \\ &= \int_F L_{x_F} \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \gamma(\star[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]) + (-1)^{k}\int_F L_{x_F} \gamma[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \gamma_{\boldsymbol{n}}(\star[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]) \\ \overset{\eqref{eqs:trace95relations}}&= \int_F L_{x_F} \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] + \int_F L_{x_F} \gamma[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star\gamma [\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \\ &\lesssim h_F \left(\int_F \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star \gamma_{\boldsymbol{n}}[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] + \int_F \gamma[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \wedge \star\gamma [\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] \right), \end{align}\] where we have expanded \(\xi\) and \(\eta\) according to the respective definitions in the second step while, in the conclusion, we have used the fact that \(L_{x_F}\lesssim h_F\) and that the integrands are both positive multiples of the volume form on \(F\). ◻
The next result shows that the ghost penalty term is enough to control the extension of \(\omega\) to the active mesh.
Lemma 4 (Global control through ghost penalty). For all \(\omega\in V_{h}^{k,r}\), it holds \[\|\omega\|^2_{\Omega_{h}} \lesssim \| \omega \|_s^2.\]
Proof. Let \(\mathcal{T}_h^{\text{im}}\mathrel{\vcenter{:}}= \mathcal{T}_h\setminus \mathcal{T}_h^{\text{cut}}\) be the set of immersed elements. We note that the norm on fully immersed elements is already controlled by the norm on \(\Omega\): \[\begin{align} \sum_{T\in\mathcal{T}_h^{\text{im}}}\|\omega\|^2_{T} \leq \|\omega\|^2_\Omega. \end{align}\] It remains to evaluate the norm on cut elements.
By Assumption 1-(ii), starting from a cut element, the number of elements we must pass through to reach a fully immersed element is at most \(N\), see 2 for an illustration in the case \(N=3\). Let \(g : \mathcal{T}_h^{\text{cut}}\to \mathcal{T}_h^{\text{im}}\) be the map which associates to each cut element \(T \in \mathcal{T}_h^{\text{cut}}\) the first fully immersed element \(g(T) \in \mathcal{T}_h^{\text{im}}\) which is reached by passing through at most \(N\) cut elements. For a given \(T'\in\mathcal{T}_h\), any \(T\in\mathcal{T}_h\) such that \(T'=g(T)\) lies in a chain of elements of cardinality \(\le N\); by Remark 2, all elements along this chain have diameters that are comparable to \(h_{T'}\), and therefore lie in a ball of radius \(\lesssim N h_{T'}\lesssim h_{T'}\). As this ball has measure \(\eqsim h_{T'}^n\) and each \(T\) in it has measure \(\eqsim h_T^n\eqsim h_{T'}^n\), we infer that \[\label{eq:estimate46gT} |g^{-1}(T')|\lesssim \frac{h_{T'}^n}{h_{T'}^n}= 1.\tag{18}\]
By applying 3 at most \(N\) times, we then have \[\begin{align} \sum_{T\in\mathcal{T}_h^{\text{cut}}} \|\omega\|^2_{T} &\lesssim \sum_{T\in\mathcal{T}_h^{\text{cut}}} \|\omega\|^2_{g(T)} + \sum_{F\in\mathcal{F}^\partial_h} \sum_{\ell=0}^r h_F^{2\ell+1} \int_F\left(\gamma_{\boldsymbol{n}}[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]\wedge\star\gamma_{\boldsymbol{n}}[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega] + \gamma[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]\wedge\star\gamma[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]\right) \\ \overset{\eqref{eq:estimate46gT}}&\lesssim \sum_{T'\in\mathcal{T}_h^{\text{im}}} \|\omega\|^2_{T'} + \sum_{F\in\mathcal{F}^\partial_h} \sum_{\ell=0}^r h_F^{2\ell+1} \int_F\left(\gamma_{\boldsymbol{n}}[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]\wedge\star\gamma_{\boldsymbol{n}}[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]+ \gamma[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]\wedge\star \gamma[ \nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]\right). \end{align}\] This completes the proof since \(\sum_{T'\in\mathcal{T}_h^{\text{im}}} \|\omega\|^2_{T'} \le \|\omega\|^2_{\Omega}\). ◻
Proof of 1. Choose an arbitrary \(\omega \in V_{h}^{k,r}\). In view of Lemma 4, we only have to show that \(\| \omega \|_s \lesssim \|\omega\|_{\Omega_{h}}\). That \(\|\omega\|_\Omega \lesssim \|\omega\|_{\Omega_{h}}\) is clear since \(\Omega \subset \Omega_h\). Standard inverse and trace inequalities for polynomials grant \[\|\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega\|^2_{F} \lesssim h_F^{-2\ell} \|\omega\|^2_{F} \lesssim h_F^{-(2\ell+1)} \|\omega\|^2_{T},\] where \(T\) is an element having \(F\) as facet (note that \(h_F \eqsim h_T\)). Hence, using a triangle inequality followed by the above relation, the norm of the jumps can be estimated as follows: denoting by \(T_1\) and \(T_2\) the two elements on each side of \(F\), \[\begin{align} \|[\nabla_{{\boldsymbol{n}}}^{(\ell)}\omega]\|_F^2 \lesssim \sum_{T\in \{T_1,T_2\}} h_F^{-(2\ell+1)} \|\omega\|^2_{T}, \end{align}\] from which \(\| \omega \|_s \lesssim \|\omega\|_{\Omega_{h}}\) follows. ◻
Remark 5 (Macro stabilisation). Notice that the proof of 4 does not require the stabilisation term to act on all facets in \(\mathcal{F}^\partial_h\). It is enough to consider the facets that form a path from each cut element to a fully immersed element, as in the definition of the map \(g\) in the proof of 4. Such a stabilisation term contains fewer summands, and leads naturally to the macro stabilisation procedure outlined in for example [11].
5.3 Hodge decomposition↩︎
Recall the definition 5 of the space \(\mathfrak{H}_h^k\) of discrete harmonic forms for FEEC, replacing \(\Omega\) by \(\Omega_{h}\). This space is constructed using orthogonality with respect to the standard inner product of \(L^2 \Lambda^k \Omega_{h}\). Since we will use the ghost product \((\bullet,\bullet)_s\) (cf. 17 ) to define numerical schemes, we need to consistently modify the space of harmonic forms as follows: \[\begin{align} \mathfrak{H}_s^k \mathrel{\vcenter{:}}= \left\{\rho_h \in V_{h}^{k,r} :\;d\rho_h=0, \;(\rho_h,d\tau_h)_s = 0 \qquad \forall \tau_h\in V_{h}^{k-1,r+1}\right\}. \end{align}\] Since both \((\bullet,\bullet)_s\) and \((\bullet,\bullet)_{\Omega_{h}}\) are inner products on \(V_{h}^{k,r}\), and both \(\mathfrak{H}^k_s\) and \(\mathfrak{H}_h^k\) are orthogonal to the same space, these spaces are isomorphic. We also have \[\begin{align} \mathop{\mathrm{Ker}}d_h = dV_{h}^{k-1,r+1} \mathbin{\raisebox{-1pt}{\begin{tikzpicture}[scale=0.134] \draw (0,-0.5)--(0,1); \draw (-0.866,-0.5)--(0.866,-0.5); \draw (0,0) circle [radius=1]; \end{tikzpicture}}}_s \mathfrak{H}_s^k, \end{align}\] where \(\mathbin{\raisebox{-1pt}{\begin{tikzpicture}[scale=0.134] \draw (0,-0.5)--(0,1); \draw (-0.866,-0.5)--(0.866,-0.5); \draw (0,0) circle [radius=1]; \end{tikzpicture}}}_s\) denotes the orthogonal complement with respect to the ghost product.
Correspondingly, we obtain the Hodge decomposition (compare with 6 ): \[V_{h}^{k,r} = (\mathop{\mathrm{Ker}}d_h)^{\perp_s} \mathbin{\raisebox{-1pt}{\begin{tikzpicture}[scale=0.134] \draw (0,-0.5)--(0,1); \draw (-0.866,-0.5)--(0.866,-0.5); \draw (0,0) circle [radius=1]; \end{tikzpicture}}}_s dV_{h}^{k-1,r+1} \mathbin{\raisebox{-1pt}{\begin{tikzpicture}[scale=0.134] \draw (0,-0.5)--(0,1); \draw (-0.866,-0.5)--(0.866,-0.5); \draw (0,0) circle [radius=1]; \end{tikzpicture}}}_s \mathfrak{H}_s^k.\]
6 Numerical validation↩︎
Consider again the Hodge Laplace equation in mixed weak form 3 . The unfitted method reads: Find \((\sigma_h,\eta_h,\lambda_h)\in V_{h}^{k-1,r+1} \times V_{h}^{k,r} \times \mathfrak{H}_s^k\) such that \[\label{eqs:unfitted95discrete95mixed} \begin{align} (\sigma_h, \tau_h)_s - (\eta_h, d\tau_h)_s &= 0 &&\qquad\forall \tau_h\in V_{h}^{k-1,r+1}, \\ (d \sigma_h, \zeta_h)_s + (d\eta_h,d\zeta_h)_s + (\lambda_h, \zeta_h)_s &= (f, \zeta_h)_\Omega &&\qquad\forall \zeta_h\in V_{h}^{k,r}, \\ (\eta_h, \rho_h)_s &= 0 &&\qquad\forall \rho_h\in \mathfrak{H}_s^k. \end{align}\tag{19}\] We deliberately take the inner product involving \(f\) over the physical domain \(\Omega\), rather than the active domain \(\Omega_{h}\), since we do not want to assume an unphysical extension of the data \(f\).
Remark 6 (Divergence constraint for the Poisson equation in mixed form). The proposed method leads to a scheme with particularly interesting properties in the case \(k = n\), which corresponds to the Poisson equation in mixed form with homogeneous Dirichlet conditions: Find \((\sigma_h,\eta_h)\in V_{h}^{n-1,r+1} \times V_{h}^{n,r}\) such that \[\begin{align} {2} (\sigma_h, \tau_h)_s - (\eta_h, \mathop{\mathrm{div}}\tau_h)_s &= 0 &&\qquad\forall \tau_h\in V_{h}^{n-1,r+1}, \\ (\mathop{\mathrm{div}}\sigma_h, \zeta_h)_s &= (f, \zeta_h)_\Omega &&\qquad\forall \zeta_h\in V_{h}^{n,r}. \end{align}\] Up to interface conditions and a source term, this is the formulation given in [10]. If \(f\in L^2\Lambda^n\Omega\) is regular enough (i.e. has \(r+1\) weak normal derivatives in \(L^2\)), then \(\mathop{\mathrm{div}}\sigma_h\) will be equal, to machine precision, to the \(L^2\)-projection of \(f\) with respect to the \((\bullet,\bullet)_s\)-inner product. Applying the same principle for Stokes flow (\(f=0\)), the corresponding scheme satisfies the incompressibility condition exactly since \(\mathop{\mathrm{div}}\sigma_h=0\). This observation, which is the content of [11], can be seen as a consequence of 1.
6.1 Implementation and results↩︎
We validate in this section the scheme 19 for \(k=1\). We respectively choose the trimmed spaces \(P^-_{1}\Lambda^{0}\Omega_{h}\) (Lagrange) and \(P^-_1\Lambda^1\Omega_{h}\) (Nédélec of the first kind) for \(\sigma_h\) and \(\eta_h\). This choice means that the exterior derivative \(d\) does not reduce the indicated polynomial order, and that the error in the natural stability norm can be expected to converge with order \(O(h)\). The scheme reads as follows. Find \((\sigma_h,\eta_h,\lambda_h)\in P^-_{1}\Lambda^{0}\Omega_{h}\times P^-_1\Lambda^1\Omega_{h}\times \mathfrak{H}_s^1\) such that \[\label{eq:discrete:k611} \begin{alignedat}{2} (\sigma_h, \tau_h)_s - (\eta_h, \nabla\tau_h)_s &= 0 &&\qquad\forall \tau_h\in P^-_1\Lambda^0\Omega_{h}, \\ (\nabla \sigma_h, \zeta_h)_s + (\mathop{\mathrm{curl}}\eta_h,\mathop{\mathrm{curl}}\zeta_h)_s + (\lambda_h, \zeta_h)_s &= (f, \zeta_h)_\Omega &&\qquad\forall \zeta_h\in P^-_1\Lambda^1\Omega_{h}, \\ (\eta_h, \rho_h)_s &= 0 &&\qquad\forall \rho_h\in \mathfrak{H}_s^1, \end{alignedat}\tag{20}\] where we use stabilisation parameter \(\eta=1\) in 16 . Since the trial and test forms are piecewise linear, their second derivatives vanish identically, and since \(\sigma_h,\tau_h\) are continuous, we get from Tables 4 and 5 that \[\begin{align} (\sigma_h, \tau_h)_s &= (\sigma_h, \tau_h)_\Omega + \sum_{F\in \mathcal{F}^\partial_h}h_F^3([\nabla\sigma_h\cdot{\boldsymbol{n}}],[\nabla\tau_h\cdot{\boldsymbol{n}}])_{F}, \\ (\eta_h,\nabla\tau_h)_s &= (\eta_h,\nabla\tau_h)_\Omega + \sum_{F\in \mathcal{F}^\partial_h}h_F([\eta_h],[\nabla\tau_h])_{F},\\ (\mathop{\mathrm{curl}}\eta_h,\mathop{\mathrm{curl}}\zeta_h)_s &= (\mathop{\mathrm{curl}}\eta_h,\mathop{\mathrm{curl}}\zeta_h)_\Omega + \sum_{F\in \mathcal{F}^\partial_h}h_F([\mathop{\mathrm{curl}}\eta_h],[\mathop{\mathrm{curl}}\zeta_h])_{F}, \\ (\lambda_h, \zeta_h)_s &= (\lambda_h, \zeta_h)_\Omega + \sum_{F\in \mathcal{F}^\partial_h}\left( h_F([\lambda_h], [\zeta_h])_{F} + h_F^3([\nabla\lambda_h\cdot{\boldsymbol{n}}], [\nabla \zeta_h\cdot{\boldsymbol{n}}])_{F} \right). \end{align}\] We also use the macro parameter \(\delta=0.25\) in the macro stabilisation procedure discussed in 5.
We solve problem 20 on the following filled torus embedded in \(\mathbb{R}^3\): \(\Omega = \{(x,y,z)\in\mathbb{R}^3: [(\sqrt{x^2+y^2}-0.5)^2 + z^2]^{1/2}\leq 0.25\}\), with major radius \(0.5\) and minor radius \(0.25\). The mesh is cut via the level set function \[\begin{align} \phi = \left[ \left(\sqrt{x^2+y^2}-0.5\right)^2 + z^2\right]^{1/2} - 0.25, \end{align}\] from a uniform background mesh with mesh size \(h\). The toroidal harmonic forms are scalar multiples of \[d\theta = \frac{xdy-ydx}{x^2+y^2} \longleftrightarrow \frac{1}{x^2+y^2}(-y,x,0).\] We pick the right-hand side and the exact solution as \[f= \left(-\frac{3xy}{(x^2+y^2)^{\frac{5}{2}}}, \frac{x^2-2y^2}{(x^2+y^2)^{\frac{5}{2}}}, 0\right), \quad\quad \eta = \left(-\frac{-xy}{(x^2+y^2)^{\frac{3}{2}}}, \frac{x^2}{(x^2+y^2)^{\frac{3}{2}}}, 0\right).\] The C++ code used to compute the example can be found in the fork [27] of the CutFEM library [28]. The code was run on a laptop with Intel i7-8565U x 8 CPU and 16GB RAM. To solve the linear system, we use the sparse direct solver MUMPS. Errors and condition numbers for three mesh sizes, \(h=\frac{1}{13}, \frac{1}{26}, \frac{1}{52},\) are shown in 3. A heat map of the \(x\)-component of the solution \(\eta\) is shown in 4. We note that the expected rate of convergence \(\mathcal{O}(h)\) is indeed achieved for the errors on the curl and gradient, with a slightly better rate for the fields themselves (going possibly towards \(\mathcal{O}(h^2)\)). The condition number for the scheme based on the ghost stabilisation remains at a reasonable magnitude, while for the unstabilized scheme (using the inner product \((\bullet,\bullet)_\Omega\) instead of \((\bullet,\bullet)_s\)) it blows up to a level where solving the linear system becomes infeasible, if not impossible.
Acknowledgements↩︎
Funded by the European Union (ERC Synergy, NEMESIS, project number 101115663). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.