A Broken FEEC Framework for Electromagnetic Problems on Mapped Multipatch Domains

Thanks to enlightening research conducted over the last few decades [1, 6, 11, 15, 33, 35, 38], it is now well understood that preserving the geometrical de Rham structure of physical problems is a key tool in the design of good finite element methods. A significant field of success is electromagnetics, where this principle has produced stable and accurate discretization methods, from simplicial Whitney forms [9, 49] to high order curved elements in isogeometric analysis [14, 27], via edge Nédélec elements [8, 42]. The latter, in particular, have been proven to yield Maxwell solvers that are free of spurious eigenmodes in a series of works dedicated to this issue [7, 10, 26, 43]. Here the existence of stable commuting projection operators plays a central role, as highlighted in the unifying analysis of finite element exterior calculus (FEEC) [1, 2].

A central asset of structure-preserving finite elements is their ability to reproduce discrete Hodge–Helmholtz decompositions which, in combination with proper commuting projection operators, allow them to preserve important physical invariants such as the divergence constraints in Maxwell’s equations [21, 22, 44], or the Hamiltonian structure of the MHD and Vlasov–Maxwell equations [24, 36, 40]. In the latter application the aforementioned commuting projection operators couple structure-preserving finite element fields with numerical particles.

As they primarily involve strong differential operators, structure-preserving finite elements have been essentially developed within the scope of conforming methods, where the discrete spaces form a sub-complex of the continuous de Rham sequence. In practice this imposes continuity conditions at the cell interfaces which strongly degrade the locality of key operations such as \(L^2\) projections, as sparse finite element mass matrices have no sparse inverses in general. In the natural framework of dual complexes composed of weak differential operators, this leads to Galerkin Hodge operators which are global, in the sense that on a given cell their values depend a priori on the function values on the whole domain. As the Hodge operators allow to map from the dual spaces to the primal ones, this results in discrete coderivatives operators that share this undesirable globality property. Furthermore the canonical commuting projection operators for the dual sequence are also global, since they coincide with the \(L^2\) projections on the discrete finite element spaces. This feature makes their application potentially expensive on fine meshes, and it represents a serious hurdle in parallel codes where communications between distant cells should be avoided as much as possible.

In this article we follow the broken FEEC approach [23] first developed for Conforming/Nonconforming Galerkin (CONGA) schemes in [18, 19]. The principle is to consider local de Rham sequences on subdomains (which we shall call patches) and discretize the global problems with broken finite element spaces, which are fully discontinuous at the patch interfaces. Strong differential operators are then obtained by composing the local (broken) ones with conforming projection operators that enforce the proper continuity conditions at the interfaces. Being fully discontinuous at the patch interfaces, the broken finite element spaces have mass matrices which are block diagonal (with blocks corresponding to patches) so that the inverse mass matrices are also block diagonal. This readily yields \(L^2\) projections and Galerkin Hodge operators which are patch-local.

Since the coupling between neighboring patches is encoded by discrete conforming projection operators which in practice involve the averaging of degrees of freedom across interfaces, the corresponding derivative and coderivative operators are local in the sense that on a given patch their values only depend on the function values on the contiguous patches. In addition our broken FEEC sequences admit dual commuting projection operators that involve \(L^2\) projections on the broken spaces and transposed conforming projection matrices [23], so these are also local.

The good news is that this approach does not sacrifice structure preservation, indeed our broken FEEC sequences satisfy primal/dual commuting diagrams, uniform stability estimates and discrete Hodge–Helmholtz decompositions [18, 23]. In particular, broken FEEC sequences should allow one to construct structure-preserving finite element solvers on complex domains with enhanced locality properties.

In the present work we follow this principle and propose specific numerical schemes for several boundary value problems arising in electromagnetics. We describe their precise form in the case of multipatch mapped spline spaces [15, 27] on general non-contractible domains, and we perform several numerical experiments to study their practical behaviour. As our results show, this approach indeed allows us to preserve most properties of the conforming FEEC approximations, such as stability and accuracy of the solutions, topological invariants, as well as harmonic and divergence constraints of the discrete fields.

The outline is as follows: We first recall in Sect. 2 the main lines of FEEC discretizations using conforming spaces and we describe their extension to broken spaces, with a detailed description of the fully discrete diagrams involving the primal (strong) and dual (weak) de Rham sequences and their respective broken commuting projection operators. In Sect. 3 we apply this discretization framework to a series of classical electromagnetic problems, namely Poisson’s and harmonic Maxwell’s equations, curl-curl eigenvalue problems, and magnetostatic problems; we also recall the approximation of the time-dependent Maxwell equations from Campos Pinto and Sonnendrücker [18]. For each problem we state a priori results about the solutions, assuming the well-posedness of the corresponding conforming FEEC discretization. In Sect. 4 we detail the construction of a geometric broken FEEC spline discretization on mapped multipatch domains. Here we consider a 2D setting for simplicity, but the same method applies to 3D domains. In Sect. 5 we conduct extensive numerical experiments for the electromagnetic problems described in the article, which verify the robustness of our approach. Finally, in Sect. 6 we summarize the main results of our work and provide an outlook on future research.

Some necessary “standard” definitions are specified in the Appendix, which provides a self-contained notation for the material of Sect. 4: Appendix Sect. A describes the construction of broken multipatch spline complexes, while Appendix Sect. B defines the geometric degrees of freedom and the corresponding commuting projection operators.

2.1 De Rham Sequences and Hodge–Helmholtz Decompositions

In this work we consider discretizations of Hilbert de Rham sequences. At the continuous level these are of the form

with infinite-dimensional spaces such as

Following the well established analysis of Hilbert complexes by Arnold, Falk and Winther [1, 2] we also consider the dual sequence

involving the adjoint differential operators (denoted with their usual name) and their corresponding domains, namely

Here the construction is symmetric, in the sense that the inhomogeneous sequence (4) could have been chosen for the primal one and the homogeneous (2) for the dual one. Since the symmetry is broken in the finite element discretization, to fix the ideas in this article we mostly consider the choice (1)–(4), except for a few places where we adopt a specific notation. A key property of these sequences is that each operator maps into the kernel of the next one, i.e., we always have \({{\,\textrm{curl}\,}}{{\,\textrm{grad}\,}}= 0\) and \({{\,\textrm{div}\,}}{{\,\textrm{curl}\,}}= 0\). This allows us to write an orthogonal Hodge–Helmholtz decomposition for \(L^2(\Omega )^3\), of the form

(5)

where \(\mathcal {H}^1:= \{ v \in V^1 \cap V^*_1: {{\,\textrm{curl}\,}}v = {{\,\textrm{div}\,}}v = 0\}\), see e.g. [2, Eq. (15)]. This space corresponds to harmonic 1-forms, as it coincides with the kernel of the 1-form Hodge–Laplace operator

seen as an operator \(\mathcal {L}^1: D(\mathcal {L}^1) \rightarrow L^2(\Omega )\) with domain space

Another orthogonal decomposition for \(L^2(\Omega )^3\) is

(8)

where \(\mathcal {H}^2:= \{w \in V^2 \cap V^*_2: {{\,\textrm{curl}\,}}w = {{\,\textrm{div}\,}}w = 0\}\) is the space of harmonic 2-forms, which coincides with the kernel of the 2-form Hodge–Laplace operator \(\mathcal {L}^{2}:= -{{\,\textrm{grad}\,}}{{\,\textrm{div}\,}}+ {{\,\textrm{curl}\,}}{{\,\textrm{curl}\,}}\) with domain space \(D(\mathcal {L}^2):= \{w \in V^2 \cap V^*_2: {{\,\textrm{curl}\,}}w \in V^*_1 \text { and } {{\,\textrm{div}\,}}w \in V^3\}\).

We point out that, while the 1-form and 2-form Hodge–Laplace operators are formally identical, their domain spaces differ in the boundary conditions and hence \(\mathcal {H}^1 \ne \mathcal {H}^2\): in the case considered here where the primal sequence has homogeneous boundary conditions, the harmonic 1-forms have vanishing tangential trace on \(\partial \Omega \), while the harmonic 2-forms have vanishing normal trace on \(\partial \Omega \). In the symmetric case where the non-homogeoneous de Rham sequence is considered as the primal one, one obtains the same decompositions (5) and (8) but with opposite order: \((\mathcal {H}^\ell )_\text {non-hom.} = \mathcal {H}^{3-\ell }\). This isomorphism, known as Poincaré duality, is provided by the Hodge star operator; see [2, Sects. 5.6, 6.2].

On contractible domains the above sequences are exact in the sense that the image of each operator coincides exactly with the kernel of the next one, and the harmonic space is trivial, \(\mathcal {H}^1 = \{0\}\). However if the domain \(\Omega \) is non-contractible, this is no longer the case and there exist non trivial harmonic forms. Indeed, the dimensions of these two harmonic spaces depend on the domain topology: \(\dim (\mathcal {H}^2) = \dim (\mathcal {H}^1)_\text {non-hom.} = b_1(\Omega )\) is the first Betti number, which counts the number of “tunnels” through the domain, while \(\dim (\mathcal {H}^1) = \dim (\mathcal {H}^2)_\text {non-hom.} = b_2(\Omega )\) is the second Betti number, which counts the number of “voids” enclosed by the domain.

2.2 Conforming FEEC Discretizations

Finite Element Exterior Calculus (FEEC) discretizations consist of Finite Element spaces that form discrete de Rham sequences,

Here, the superscript c indicates that the spaces are assumed conforming in the sense that \(V^{\ell ,c}_h \subset V^{\ell }\), and the subscript h loosely represents some discretization parameters, such as the resolution of an underlying mesh. A key tool in the analysis of FEEC discretizations is the existence of projection operators \(\Pi ^{\ell ,c}_h: V^\ell \rightarrow V^{\ell ,c}_h\) that commute with the differential operators, in the sense that the relation

holds for the different operators in the sequence (9),

In particular, the stability and the accuracy of several discrete problems posed in the sequence (9), relative to usual discretization parameters such as the mesh resolution h or the order of the finite element spaces, follow from the stability of the commuting projections in \(V^\ell \) or \(L^2\) norms, see e.g. [2, Theorems 3.9, 3.19].

Denoting explicitly by

the differential operators restricted to the discrete spaces, we define their adjoints by discrete \(L^2\) duality, i.e.,

where \(\langle \cdot ,\cdot \rangle \) denotes the \(L^2(\Omega )\) scalar product. Owing to the homogeneous boundary conditions, we denote the resulting coderivative operators by

This yields a compatible discretization of both the primal and dual sequences (1) and (3), in strong and weak form, respectively, using the same spaces (9). Note that in (13) the discrete spaces \(V_h^{\ell ,c}\) are implicitely identified with their dual spaces when equipped with the \(L^2\) norm, which is standard practice, see e.g. [2, Sec 3.3]. This framework can be summarized in the following diagram

(15)

where the operators \(Q_{V^{\ell ,c}_h}: L^2(\Omega ) \rightarrow V^{\ell ,c}_h\) represent \(L^2\) projections to the conforming discrete spaces. We observe that these projections commute with the dual differential operators, as a result of (14).

Assumption 1

Throughout the article we assume that the primal projection operators \(\Pi ^{\ell ,c}_h\) are \(L^2\) stable and satisfy the commuting property (10).

Remark 1

In practice commuting projection operators also play an important role as they permit to approximate coupling or source terms in a structure-preserving way, see e.g. [24]. For this purpose \(V^\ell \) or \(L^2\)-stable projections may not be the best choices as they can be difficult to apply, and simpler commuting projections, defined through proper degrees of freedom, are often preferred. These projections are then usually defined on sequences

involving spaces \(U^\ell \subset V^\ell \) that require more smoothness (or integrability) as the ones in (2). We refer to e.g. [5, 24, 41, 45] for some examples, and to Sect. 2.4 below for more details on such constructions.

Another asset of FEEC discretizations is to provide structure-preserving Hodge–Helmholtz decompositions for the different spaces. For 1-forms, the discrete analog of the continuous decomposition (5) reads

(17)

where \(\mathcal {H}^1_h:= \{v \in V^{1,c}_h: {{\,\textrm{curl}\,}}^c_h v = \widetilde{{\,\textrm{div}\,}}^c_h v = 0\}\) is the kernel of the discrete Hodge–Laplace operator \(\mathcal {L}^{1,c}_h: V^{1,c}_h \rightarrow V^{1,c}_h\) defined as

The space \(\mathcal {H}^1_h\) may thus be seen as discrete harmonic 1-forms and under suitable approximation properties of the discrete spaces, its dimension coincides with that of the continuous harmonic forms \(\mathcal {H}^1\), which corresponds to a Betti number of \(\Omega \) depending on the boundary conditions, see [2, Sects. 5.6, 6.2].

2.3 Broken FEEC Discretizations

In the case where the domain is decomposed in a partition of open subdomains \(\Omega _k\), \(k = 1, \dots , K\), we now consider local FEEC sequences

and global spaces obtained by a simple juxtaposition of the local ones

One attractive feature of the broken spaces (20) is that they are naturally equipped with local basis functions \(\Lambda ^\ell _i\) that are supported each on a single patch \(\Omega _k\) with \(k = k(i)\). The corresponding mass matrices are then patch-diagonal, i.e., block-diagonal with blocks corresponding to the different patches, so that their inversion – and hence the \(L^2\) projection on \(V^\ell _h\) – can be performed in each patch independently of the others.

However, the spaces (20) are in general not subspaces of their infinite-dimensional counterparts, indeed it is well-known that piecewise smooth fields must satisfy some interface constraints in order to be globally smooth [5]: for \(H^1\) smoothness the fields must be continuous on the interfaces, while \(H({{\,\textrm{curl}\,}})\) and \(H({{\,\textrm{div}\,}})\) smoothness of vector-valued fields require the continuity of the tangential and normal components, respectively, on the interfaces. As these constraints are obviously not satisfied by the broken spaces (20), we have in general

The approach developed for the CONGA schemes in [18, 23] extends the construction of Sect. 2.2 to this broken FEEC setting, by associating to each discontinuous space a projection operator on its conforming subspace,

In practice \(P^\ell _h\) can be defined by a local averaging of interface degrees of freedom.

Provided that these conforming spaces form a de Rham sequence (9), this allows us to define new primal differential operators on the broken spaces

and new dual ones \(\widetilde{{\,\textrm{div}\,}}_h:V^{1}_h \rightarrow V^{0}_h\), \(\widetilde{{\,\textrm{curl}\,}}_h: V^{2}_h \rightarrow V^{1}_h\) and \(\widetilde{{\,\textrm{grad}\,}}_h: V^{3}_h \rightarrow V^{2}_h\) as discrete \(L^2\) adjoints, characterized by the relations

We represent this broken FEEC discretization with the following diagram

(24)

where \(\Pi ^\ell _h\) and \({\tilde{\Pi }}^\ell _h\) are projection operators that commute with the new primal and dual sequences. Here \(\Pi ^\ell _h\) are seen as unbounded operators with dense domains \(U^\ell \subset V^\ell \). For the dual sequence the commuting projections \({\tilde{\Pi }}^\ell _h\) will be bounded in \(L^2\): We postpone their description to the next sections.

We note that the presence of conforming projections \(P^\ell _h\) in the differential operators (22) leads to larger kernels. Specifically, we have

where the projection kernels

correspond intuitively to “jump spaces” associated with the conforming projections. In [17, 19] these extended kernels motivated a modification of the discrete differential operators in order to retain exact sequences on contractible domains. Here we follow the approach of [23] where the jump spaces (25) are handled by stabilisation and filtering operators in the equations, through extended Hodge–Helmholtz decompositions.

2.4 Commuting Projections with Broken Degrees of Freedom

A practical approach for designing commuting projection operators \(\Pi ^{\ell ,c}_h\) and \(\Pi ^\ell _h\) for the above diagrams is to define them via commuting degrees of freedom on the discrete Finite Element spaces. In the standard case of conforming spaces these are linear forms

defined on infinite-dimensional spaces \(U^\ell \subset V^\ell \) satisfying \(d^\ell U^\ell \subset U^{\ell +1}\), that are unisolvent in \(V^{\ell ,c}_h\) and commute with the differential operators in the sense that there exist coefficients \(D^\ell _{i,j}\) such that

where we have denoted again \(d^0 = {{\,\textrm{grad}\,}}\), \(d^1 = {{\,\textrm{curl}\,}}\) and \(d^2={{\,\textrm{div}\,}}\). Letting \(\Lambda ^{\ell ,c}_i\) be the basis of \(V^{\ell ,c}_h\) characterized by the relations \(\sigma ^{\ell ,c}_i(\Lambda ^{\ell ,c}_j) = \delta _{i,j}\) for \(i,j = 1, \dots N^{\ell ,c}\), one verifies indeed that the operators

are projections satisfying the commuting property (10), see e.g. [24, Lemma 2] and in passing we note that \(D^\ell \) is the matrix of \(d^\ell \) in the bases just defined.

Here, we extend this approach as follows:

• Broken degrees of freedom. Each local space \(V^{\ell }_h(\Omega _k)\), \(k = 1, \ldots , K\), is equipped with unisolvent degrees of freedom

  $$\begin{aligned} \sigma ^{\ell }_{k,\mu }: U^\ell (\Omega _k) \rightarrow \mathbb {R}, \quad \mu \in \mathcal {M}^\ell _h(\Omega _k) \end{aligned}$$

  (27)

  defined on local spaces \(U^\ell (\Omega _k) \subset V^\ell (\Omega _k)\) satisfying \(d^\ell U^\ell (\Omega _k) \subset U^{\ell +1}(\Omega _k)\), with multi-index sets with cardinality \(\#\mathcal {M}^\ell _h(\Omega _k) = \dim (V^{\ell }_h(\Omega _k))\). On the full domain we simply set

  $$\begin{aligned} \sigma ^{\ell }_{k,\mu }(v):= \sigma ^{\ell }_{k,\mu }(v|_{\Omega _k}). \end{aligned}$$

  (28)

• Broken basis functions. To the above degrees of freedom we associate local basis functions \(\Lambda ^\ell _{k,\mu } \in V^\ell _h(\Omega _k)\) (extended by zero outside of their patch), characterized by the relations

  $$\begin{aligned} \sigma ^{\ell }_{k,\mu }(\Lambda ^{\ell }_{k,\nu }) = \delta _{\mu ,\nu } \quad \text { for } \quad \mu , \nu \in \mathcal {M}^\ell _h(\Omega _k). \end{aligned}$$

• Local commutation property. A relation similar to (26) must hold on each patch \(\Omega _k\), i.e., there exist coefficients \(D^\ell _{k,\mu ,\nu }\) such that

  $$\begin{aligned} \sigma ^{\ell +1}_{k,\mu }(d^\ell v) = \sum _{\nu \in \mathcal {M}^\ell _h(\Omega _k)} D^\ell _{k,\mu ,\nu } \sigma ^{\ell }_{k,\nu }(v) \end{aligned}$$

  (29)

  holds for all \(\mu \in \mathcal {M}^{\ell +1}_h(\Omega _k)\) and all \(v \in U^\ell (\Omega _k)\).

• Inter-patch conformity. The global projection on \(V^{\ell }_h\),

  $$\begin{aligned} \Pi ^{\ell }_h: v \mapsto \sum _{k=1}^K \Pi ^{\ell }_{h,k} v \quad \text { with } \quad \Pi ^{\ell }_{h,k}: v \mapsto \sum _{\mu \in \mathcal {M}^{\ell }_h(\Omega _k)} \sigma ^{\ell }_{k,\mu }(v) \Lambda ^{\ell }_{k,\mu } \end{aligned}$$

  (30)

  must map smooth functions to conforming finite element fields, namely

  $$\begin{aligned} \Pi ^\ell _h U^\ell \subset V^{\ell ,c}_h \quad \text { where } \quad U^\ell := \{v \in V^\ell : v|_{\Omega _k} \in U^\ell (\Omega _k), ~ \forall k = 1, \dots K\}. \end{aligned}$$

  (31)

This setting guarantees that the commutation properties of the local projection operators extend to the global ones, since one has then

and it yields simple expressions for the latter in the broken bases. We refer to Sect. 4 for a detailed construction of broken degrees of freedom that satisfy the above properties.

The dual projection operators (namely, the projection operators on the finite element spaces which commute with the differential operators from the dual sequences) can then be defined following the canonical approach of Campos-Pinto and Güçlü [23], as “filtered” \(L^2\) projections

where \(Q_{V^\ell _h}\) is now the \(L^2\) projection on the broken space \(V^\ell _h\),

By definition of the dual differentials (23), these operators indeed commute with the dual part of the diagram (24), in the sense that

hold on \(V^*_1\), \(V^*_2\) and \(V^*_3\) respectively. We further observe that (32)–(33) defines projection operators on the subspaces of \(V^\ell _h\) corresponding to the orthogonal complements of the “jump spaces” (25). Indeed, \(({\tilde{\Pi }}^\ell _h)^2 = {\tilde{\Pi }}^\ell _h\) holds with

2.5 Differential Operator Matrices in Primal and Dual Bases

Using broken degrees of freedom and basis functions as described in the previous section, we now derive practical representations for the operators involved in the diagram (24). To do so we first observe that the non-conforming differential operators \(d^\ell _h:= d^\ell P^\ell _h\) can be reformulated as \(d^\ell _h = d^\ell _{\textrm{pw}} P^\ell _h\) where \(d^\ell _{\textrm{pw}}: v \rightarrow \sum _{k=1}^K {\mathbb {1}}_{\Omega _k} d^\ell |_{\Omega _k} v\) is the patch-wise differential operator which coincides with \(d^\ell \) on all v such that \(d^\ell v \in L^2\). Using some implicit flattening \((k,\mu ) \mapsto i \in \{1, \dots , N^\ell \}\) for the multi-indices, we represent these operators as matrices

of respective sizes \(N^1 \times N^0\), \(N^2 \times N^1\) and \(N^3 \times N^2\). These matrices have a “patch-diagonal” structure, in the sense that they are block-diagonal, with blocks corresponding to the differential operators in each independent patch [namely, the matrices \(D^\ell _k\) in (29)]. The matrices of the conforming projections, seen as endomorphisms in \(V^\ell _h\), are

In general \(\mathbb {P}^\ell \) is not patch-diagonal, as it maps in the coefficient space of the conforming spaces \(V^{\ell ,c}_h\), where the global smoothness corresponds to some matching of the degrees of freedom across the interfaces.

With these elementary matrices we can build the operator matrices of the non-conforming differential operators \(d^\ell _h\): they read

where we have denoted \(({\mathbb {d}}^\ell )_\ell = (\mathbb {G}, \mathbb {C}, \mathbb {D})\) for brevity. Therefore,

A simple representation of the dual discrete operators (23) is obtained in the dual bases \(\{{\tilde{\Lambda }}^\ell _i: i = 1, \dots , N^\ell \}\) of the broken spaces \(V^\ell _h\). These are characterized by the relations

with dual degrees of freedom defined as

where \(\langle \cdot ,\cdot \rangle \) denotes again the \(L^2\) scalar product. It follows from (38) and (39) that the primal and dual bases are in \(L^2\) duality, i.e. \(\langle \Lambda ^\ell _i,{\tilde{\Lambda }}^\ell _j\rangle = \delta _{i,j}\). The matrices of the dual differential operators in the dual bases read then

where again \((\mathbb {D}^\ell )_\ell = (\mathbb {G}, \mathbb {C}, \mathbb {D})\) for brevity. Therefore,

The change-of-basis matrices are described in the next section.

2.6 Primal-Dual Diagram in Matrix Form

Using the matrix form of the differential operators just described, we extend the broken FEEC diagram (24) as follows:

(41)

Here \(\mathcal {C}^\ell \) and \({\tilde{\mathcal {C}}}^\ell \) are the spaces of scalar coefficients associated with the primal and dual basis functions described in the previous Section. Both are of the form \(\mathbb {R}^{N^\ell }\) with \(N^\ell = \dim (V^\ell _h)\), but we use a different notation to emphasize the different roles played by the coefficient vectors. The interpolation operators \(\mathcal {I}^\ell : \mathcal {C}^\ell \rightarrow V^\ell _h\) and \({\tilde{\mathcal {I}}}^\ell : {\tilde{\mathcal {C}}}^\ell \rightarrow V^\ell _h\), which read

are the right-inverses of the respective degrees of freedom \({\varvec{\mathsf {\sigma }}}^\ell \) and \({\tilde{{\varvec{\mathsf {\sigma }}}}}^\ell \), and also their left-inverses on \(V^\ell _h\).

The matrices \(\mathbb {H}^\ell \) and \({\tilde{\mathbb {H}}}^\ell \) are change-of-basis matrices which allow us to go from one sequence to the other, in the sense that \(\mathcal {I}^\ell ={\tilde{\mathcal {I}}}^\ell \mathbb {H}^\ell \) and \({\tilde{\mathcal {I}}}^\ell = \mathcal {I}^\ell {\tilde{\mathbb {H}}}^\ell \): they correspond to discrete Hodge operators [34]. Here it follows from the duality construction (38)–(39) that they are given by the mass matrices and their inverses,

In our framework, the Hodge matrices have a patch-diagonal structure due to the local support of the broken basis functions. In particular, the dual basis functions are also supported on a single patch.

The remaining operators are the primal and dual commuting projections, respectively defined by (30) and (32)–(33), expressed in terms of primal and dual coefficients. Using the interpolation operators (42) we may write them as

Finally we remind that here the commutation of the primal (upper) diagram follows from the assumed properties of the primal degrees of freedom (27)–(31), while the commutation of the dual (lower) diagram follows from the weak definition of the dual differential operators.

Remark 2

(Conforming case) The conforming FEEC diagram (15) (corresponding, e.g., to a single patch discretization) may be extended in the same way as the broken FEEC one (24). Apart from the differences between the discrete differential operators and commuting projection operators already visible in (15) and (24), the resulting 6-rows diagram would be formally the same as (41), but no conforming projection matrices would be involved any longer.

Remark 3

(Change of basis) It may happen that the basis functions \(\Lambda ^\ell _i\) associated with the commuting degrees of freedom \(\sigma ^\ell _i\) are not the most convenient ones when it comes to actually implement the discrete operators. In such cases a few changes are to be done to the diagram (41), such as the ones described in Sect. 4.2 for local spline spaces.

2.7 Discrete Hodge–Laplace and Jump Stabilization Operators

In addition to the first order differential operators, an interesting feature of the diagram (41) is to provide us with natural discretizations for the Hodge–Laplace operators. On the space \(V^0\) this is the standard Laplace operator, \(\mathcal {L}^0 = -\Delta \), which is discretized as

with an operator matrix in the primal basis that reads

As for the Hodge–Laplace operator for 1-forms, \(\mathcal {L}^1 = {{\,\textrm{curl}\,}}{{\,\textrm{curl}\,}}- {{\,\textrm{grad}\,}}{{\,\textrm{div}\,}}\), its discretization is

with an operator matrix (again in the primal basis)

The Hodge–Laplace operators for 2 and 3-forms are discretized similarly.

As discussed in Sect. 2.2, a key asset of FEEC discretizations is their ability to preserve the exact dimension of the harmonic forms, defined here as the kernel of the Hodge–Laplace operators. In our broken FEEC framework this property is a priori not preserved by the non-conforming operators due to the extended kernels of the conforming projection operators, but it is for the stabilized operators \(\mathcal {L}^\ell _{h,\alpha } = \mathcal {L}^\ell _h + \alpha S^\ell _h\) where

is the symmetrized projection operator on the jump space (25), and \(\alpha \) is a stabilization parameter. Indeed, it was shown in [23] that the kernel of \(\mathcal {L}^\ell _{h,\alpha }\) coincides with that of the conforming discrete Hodge–Laplacian, for any positive value of \(\alpha > 0\). Finally we note that the corresponding stabilization matrix is readily derived from the operators in the diagram. It reads

In this section we propose broken FEEC approximations for several problems arising in electromagnetics, and we state a priori results for the solutions. Let us begin by listing the main differences between broken and conforming FEEC approximations:

• Differential and projection operators from the commuting diagram (41) apply to broken (fully discontinuous) spaces and involve conforming projection operators which in practice couple interface degrees of freedom and can be applied matrix free.

• For problems that involve solving a system (e.g., Poisson, time-harmonic Maxwell and Magnetostatic equations), symmetric stabilization operators must be added to make the system non singular. This is not needed for the time-dependent Maxwell equations.

• For the former problems, our stabilized broken FEEC schemes yield the same solution as the corresponding conforming FEEC method. This is not the case for the time-dependent Maxwell equations.

• For eigenvalue problems several formulations can be used (with or without stabilization terms), depending on whether zero or positive eigenmodes are searched for.

For the sake of completeness we will provide several formulations (all equivalent) for the discrete problems: one in operator form, one in weak form and one in matrix form, using the primal and dual bases introduced in Sect. 2.4. As pointed out in Remark 3, practical implementation may be more efficient with other basis functions. In this case the matrices need to be changed, as will be described in Sect. 4.2 below.

Throughout this section we will mostly work with the homogeneous spaces \(V^\ell \) corresponding to (2). In the few cases where we will need the inhomogeneous spaces, we shall use a specific notation \({\bar{V}}^\ell \) for the spaces in (4), for the sake of clarity.

3.1 Poisson’s Equation

For the Poisson problem with homogeneous Dirichlet boundary conditions: given \(f \in L^2(\Omega )\), find \(\phi \in V^0 = H^1_0(\Omega )\) such that

we combine the stabilized broken FEEC discretization described in Sect. 2.7 for the Laplace operator, and a dual commuting projection (32)–(33) for the source. The resulting problem reads: Find \(\phi _h \in V^0_h\) such that

and it enjoys the following property.

Proposition 1

For all \(\alpha \ne 0\), Eq. (52) admits a unique solution \(\phi _h\) which belongs to the conforming space \(V^{0,c}_h = V^0_h \cap H^1_0(\Omega )\) and solves the conforming Poisson problem

In particular, \(\phi _h\) is independent of both \(\alpha \) and the specific projection \(P^0_h\).

Proof

By definition of the different operators, (52) reads

which may be reformulated using test functions \({\varphi }\in V^{0}_h\), as

Taking \({\varphi }= (I-P^0_h)\phi _h\) yields \(\phi _h = P^0_h \phi _h \in V^{0,c}_h\) as long as \(\alpha \ne 0\), and (53) follows by considering test functions \({\varphi }\) in the conforming subspace \(V^{0,c}_h\). \(\square \)

The matrix formulation of (52) is easily derived by using the diagram (41), writing the equation in the dual basis in order to obtain symmetric matrices as is usual with the Poisson problem. Denoting by \({\varvec{\mathsf {\phi }}} = {\varvec{\mathsf {\sigma }}}^0(\phi _h)\) the coefficient vector of \(\phi _h\) in the primal basis, we thus find

with a stabilized stiffness matrix

3.2 Time-Harmonic Maxwell’s Equation

The time-harmonic Maxwell equation with homogeneous boundary conditions reads: given \(\omega \in \mathbb {R}\) and \(J \in L^2(\Omega )\), find \(u \in H_0({{\,\textrm{curl}\,}};\Omega )\) such that

Here the (complex) electric field corresponds to \(E(t,x) = i\omega \, u(x) e^{-i \omega t}\). When \(\omega ^2\) is not an eigenvalue of the \({{\,\textrm{curl}\,}}{{\,\textrm{curl}\,}}\) operator, Eq. (58) is well-posed: see e.g. [3, Theorem 8.3.3]. For this problem we propose a stabilized broken FEEC discretization where \(u_h \in V^1_h\) solves

with a parameter \(\alpha \in \mathbb {R}\). We remind that \(S^1_h:= (1-P^1_h)^* (1-P^1_h)\) is the jump stabilization operator, according to (49) and (25) for \(\ell =1\). With test functions \(v \in V^{1}_h\), the discrete problem (59) writes

Here, we note that the zeroth-order term is filtered by the symmetric operator \((P^1_h)^* P^1_h\): this allows us to obtain a conforming solution in the broken space.

Proposition 2

Let \(\alpha \ne 0\), and \(\omega \) such that \(\omega ^2\) is not an eigenvalue of the conforming discrete \(\widetilde{{\,\textrm{curl}\,}}^c_h {{\,\textrm{curl}\,}}^c_h\) operator defined in (12)–(14). Then Eq. (59) admits a unique solution \(u_h\) which belongs to the conforming space \(V^{1,c}_h = V^{1}_h \cap H_0({{\,\textrm{curl}\,}};\Omega )\), and solves the conforming FEEC Maxwell problem

In particular, \(u_h\) is independent of both \(\alpha \) and the specific projection \(P^1_h\).

Proof

Taking \(v = (I-P^1_h)u_h\) in (60) gives \(u_h = P^1_h u_h \in V^{1,c}_h\) as long as \(\alpha \ne 0\), and (61) follows by considering test functions \(v = v^c\) in \(V^{1,c}_h\). Existence and uniqueness follow from the assumption that \(\omega \) is not an eigenvalue of the conforming curl curl operator. \(\square \)

The matrix form of (59) is easily derived from diagram (41). It reads

where \({\varvec{\textsf{u}}}:= {\varvec{\mathsf {\sigma }}}^1(u_h)\) is the (column) vector containing the coefficients of \(u_h\) in the primal basis of \(V^1_h\), and \(\mathbb {A}^1\) is the resulting stabilized stiffness matrix,

3.3 Lifting of Boundary Conditions

In the case of inhomogeneous boundary conditions a standard approach is to introduce a lifted solution and solve the modified problem in the homogeneous spaces. This approach applies seamlessly to broken FEEC approximations. For a Poisson equation of the form

this corresponds to introducing \(\phi _g \in H^1(\Omega )\) such that \(\phi _g = g\) on \(\partial \Omega \), and characterizing the solution as \(\phi = \phi _g + \phi _0\), where \(\phi _0 \in V^0 = H^1_0(\Omega )\) solves (51) with a modified source, namely

For an inhomogeneous Maxwell equation of the form

we introduce \(u_g \in H({{\,\textrm{curl}\,}};\Omega )\) such that \(n \times u_g = g\) on \(\partial \Omega \), and characterize the solution as \(u = u_g + u_0\), where \(u_0 \in H_0({{\,\textrm{curl}\,}};\Omega )\) solves (58) with a modified source, namely

The lifting approach can be applied in broken FEEC methods by combining the discretizations of the homogeneous and inhomogeneous sequences (2) and (4), which amounts to combining different conforming projections. For clarity we denote in this section the respective homogeneous and inhomogeneous spaces by \(V^\ell \) and \({\bar{V}}^\ell \). At the discrete level the broken spaces \(V^\ell _h\) have no boundary conditions, so that the distinction only appears in the conforming projections to the spaces \(V^{\ell ,c}_h = V^\ell _h \cap V^\ell \) and \({\bar{V}}^{\ell ,c}_h = V^\ell _h \cap {\bar{V}}^\ell \), which we naturally denote by \(P^\ell _h\) and \({\bar{P}}^\ell _h\) respectively. In practice the latter projects on the inhomogeneous conforming spaces, while the former further sets the boundary degrees of freedom to zero.

For the Poisson equation where \(V^0 = H^1_0(\Omega )\) and \({\bar{V}}^0 = H^1(\Omega )\), we first compute \(\phi _{g,h} \in V^{0}_h\) that approximates g on \(\partial \Omega \), and then we compute the homogeneous part of the solution, \( \phi _{0,h}:= \phi _h - \phi _{g,h} \in V^{0}_h, \) by solving the homogeneous CONGA problem with modified source

for all \({\varphi }\in V^{0}_h\), which also reads in matrix form

where \(\mathbb {A}^0\) is the stabilized stiffness matrix (57), \({\bar{\mathbb {P}}}^0\) is the matrix of the inhomogeneous projection \({\bar{P}}^0_h\), and the coefficient vectors involve the broken degrees of freedom \({\varvec{\mathsf {\phi }}}_0 = {\varvec{\mathsf {\sigma }}}^0(\phi _{0,h})\), \({\varvec{\mathsf {\phi }}}_g = {\varvec{\mathsf {\sigma }}}^0(\phi _{g,h})\) in the full broken space \(V^0_h\).

For the Maxwell equation where \(V^1 = H_0({{\,\textrm{curl}\,}};\Omega )\) and \({\bar{V}}^1 = H({{\,\textrm{curl}\,}};\Omega )\), the method consists of first computing \(u_{g,h} \in V^{1}_h\) such that \(n \times u_{g,h}\) approximates g on \(\partial \Omega \), and then characterizing the homogeneous part of the solution, \( u_{0,h}:= u_h - u_{g,h} \in V^{1}_h, \) by

for all \(v \in V^{1}_h\). In matrix terms, this reads

where \(\mathbb {A}^1\) is the stabilized stiffness matrix (63), \({\bar{\mathbb {P}}}^1\) is the matrix of the inhomogeneous projection \({\bar{P}}^1_h\), and the coefficient vectors involve the broken degrees of freedom \({\varvec{\textsf{u}}}_0 = {\varvec{\mathsf {\sigma }}}^1(u_{0,h})\), \({\varvec{\textsf{u}}}_g = {\varvec{\mathsf {\sigma }}}^1(u_{g,h})\) in the full broken space \(V^1_h\).

Proposition 3

Let \(\alpha \) and \(\omega \) be as in Proposition 2. Then (69) admits a unique solution \(u_{0,h}\) which belongs to the conforming space \(V^{1,c}_h = V^{1}_h \cap H_0({{\,\textrm{curl}\,}};\Omega )\). The projection of the full solution in \({\bar{V}}^{1,c}_h = V^{1}_h \cap H({{\,\textrm{curl}\,}};\Omega )\), \(u^c_h:= {\bar{P}}^1_h u_h\), approximates the boundary condition, \(n \times u^c_h \approx g\) on \(\partial \Omega \), and it solves the discrete conforming Maxwell equation inside \(\Omega \), in the sense that

In particular, \(u^c_h\) is independent of both \(\alpha \) and \(P^1_h\). Moreover if the lifted boundary condition \(u_{g,h}\) is chosen in \({\bar{V}}^{1,c}_h\), then \(u_h = u^c_h\).

Proof

The conformity and unicity of \(u_{0,h}\) can be shown with the same arguments as in Proposition 2. This shows that \(u^c_h = u_{0,h} + {\bar{P}}^1_h u_{g,h}\) and that the tangential trace of \(u^c_h\) coincides with that of \({\bar{P}}^1_h u_{g,h}\) on the boundary, which itself should coincide with that of \(u_{g,h}\) as there are no conformity constraints there for the inhomogeneous space \({\bar{V}}^1 = H({{\,\textrm{curl}\,}};\Omega )\). Equation (71) and the additional observations follow easily. \(\square \)

3.4 Eigenvalue Problems

Broken FEEC approximations to Hodge–Laplace eigenvalue problems of the form

are readily derived using the discrete operators presented in Sect. 2.7. For a general stabilization parameter \(\alpha \ge 0\), they take the form:

with \(u_h \in V^\ell _h \setminus \{0\}\). Their convergence can be established under the assumption of a strong stabilization regime and of moment-preserving properties for the conforming projections, see [23]. There it has also been shown that for arbitrary positive penalizations \(\alpha > 0\), the zero eigenmodes coincide with their conforming counterparts, namely the conforming harmonic forms,

see (17)–(18).

For the first space \(V^0 = H^1_0(\Omega )\) the corresponding operator is the standard Laplace operator \(\mathcal {L}^0 = -\Delta \) and its broken FEEC discretization reads: Find \(\lambda _h \in \mathbb {R}\) and \(\phi _h \in V^0_h {\setminus } \{0\}\) such that

where \(\mathcal {L}^0_h\) and \(S^0_h\) are given by (45) and (49), with stabilization parameter \(\alpha \). For \(V^1 = H_0({{\,\textrm{curl}\,}};\Omega )\) the Hodge–Laplace operator is \(\mathcal {L}^1 = {{\,\textrm{curl}\,}}{{\,\textrm{curl}\,}}-{{\,\textrm{grad}\,}}{{\,\textrm{div}\,}}\) and its discretization reads: Find \(\lambda _h \in \mathbb {R}\) and \(u_h \in V^1_h {\setminus } \{0\}\) such that

where \(\mathcal {L}^1_h = \widetilde{{\,\textrm{curl}\,}}_h {{\,\textrm{curl}\,}}_h - {{\,\textrm{grad}\,}}_h \widetilde{{\,\textrm{div}\,}}_h\), see (47), and \(S^1_h\) is given again by (49). Eigenvalue problems in \(V^2_h\) and \(V^3_h\) are discretized in a similar fashion.

Another important eigenvalue problem is that of the curl-curl operator, for which one could simply adapt the broken FEEC discretization (75) used for the Hodge–Laplace operator, by dropping the term \((-{{\,\textrm{grad}\,}}_h \widetilde{{\,\textrm{div}\,}}_h)\). By testing such an equality with \(v \in {{\,\textrm{grad}\,}}P^0_h V^0_h\) we see that the non-zero eigenmodes satisfy \(\widetilde{{\,\textrm{div}\,}}_h u_h = 0\), hence they are also eigenmodes of (75) and their convergence follows from the results mentionned above in the case of strong penalization regimes. In the unpenalized case (\(\alpha = 0\)) their convergence was also established, as shown in [18].

The main drawback of the approach just described is that the computed eigenvalues do not coincide with those of the conforming operator \(\widetilde{{\,\textrm{curl}\,}}^c_h {{\,\textrm{curl}\,}}^c_h\), but only approximate them. (We recall that the knowledge of the eigenvalues of \(\widetilde{{\,\textrm{curl}\,}}^c_h {{\,\textrm{curl}\,}}^c_h\) is needed in order to verify the hypotheses of Proposition 2 for the source problem: \(\omega ^2\) should not concide with one of them.) To overcome this difficulty we propose a novel broken FEEC discretization, which consists of the generalized eigenvalue problem

With test functions \(v \in V^1_h\), the discrete problem (76) writes

and its matrix formulation is provided below.

Proposition 4

For \(\lambda _h \ne 0\), the (generalized) CONGA eigenvalue problem (76) has the same solutions as the conforming eigenproblem

For \(\lambda _h = 0\), the corresponding eigenspace is

(79)

Proof

For \(\lambda _h \ne 0\) we take \(v = (I - P^1_h) u_h\) in (77), and find that \(u_h = P^1_h u_h \in V^{1,c}_h\). Hence we can restrict the test space to functions \(v = P^1_h v \in V^{1,c}_h\), and obtain the conforming eigenproblem \(\langle {{\,\textrm{curl}\,}}v, {{\,\textrm{curl}\,}}u_h \rangle = \lambda _h \langle v, u_h \rangle \). Conversely, we verify that if \(u_h \in V^{1,c}_h\) solves (78) then it also solves (77) for all \(v \in V^1_h\). The relationship (79) between the null spaces follows from the equalities \(\ker \big (\widetilde{{\,\textrm{curl}\,}}_h {{\,\textrm{curl}\,}}_h \big ) = \ker ({{\,\textrm{curl}\,}}_h)\), and

which was observed in, e.g., [18]. \(\square \)

These eigenvalue problems may be expressed in symmetric matrix form by expressing the eigenmodes in the primal bases and the equations in the dual bases, as we did for the Poisson and Maxwell equations in Sects. 3.1 and 3.2. For the Poisson eigenvalue problem (74) we obtain

with \(\mathbb {A}^0 = \mathbb {H}^0(\mathbb {L}^0 + \alpha \mathbb {S}^0) = (\mathbb {G}\mathbb {P}^0)^T \mathbb {H}^1 \mathbb {G}\mathbb {P}^0 + \alpha (\mathbb {I}-\mathbb {P}^0)^T \mathbb {H}^0 (\mathbb {I}-\mathbb {P}^0)\) as in (57). Similarly we can write the Hodge–Laplace eigenvalue problem (75) as

with matrix

and the curl-curl eigenvalue problem (76) as

3.5 Magnetostatic Problems with Harmonic Constraints

We next consider a magnetostatic problem of the form

posed for a current \(J \in L^2(\Omega )\). Unlike the Poisson and Maxwell source problems above, this problem requires an additional constraint in non-contractible domains. Indeed it can only be well-posed if \(\ker {{\,\textrm{curl}\,}}\cap \ker {{\,\textrm{div}\,}}= \{0\}\), i.e., if the harmonic forms are trivial. To formulate this more precisely we need to take into account particular boundary conditions.

3.5.1 Magnetostatic Problem with Pseudo-Vacuum Boundary Condition

We first consider a boundary condition of the form \(n \times B = 0\), which is sometimes used to model pseudo-vacuum boundaries [39]. A mixed formulation is then: Find \(B \in H_0({{\,\textrm{curl}\,}};\Omega )\) and \(p \in H^1_0(\Omega )\) such that

see e.g. [28]. Here the auxiliary unknown p is a Lagrange multiplier for the divergence constraint written in weak form.

The well-posedness of (84) is related to the space of harmonic 1-forms \(\mathcal {H}^1\), defined as the kernel of the Hodge–Laplace operator (6) in the homogeneous sequence, i.e.,

As discussed in [2], \(\mathcal {H}^1 = \mathcal {H}^1(\Omega )\) is isomorphic to a de Rham cohomology group and its dimension corresponds to a Betti number of the domain \(\Omega \). With homogeneous, resp. inhomogeneous boundary conditions its dimension is \(b_2(\Omega )\) the number of cavities, resp. \(b_1(\Omega )\) the number of handles in \(\Omega \). Thus in a contractible domain we have \(\mathcal {H}^1 = \{0\}\) and the problem is well posed, but in general \(\mathcal {H}^1\) may not be trivial. Additional constraints must then be imposed on the solution, such as \( B \in (\mathcal {H}^1)^\perp , \) which can be associated with an additional Lagrange multiplier \(z \in \mathcal {H}^1\). The constrained problem then consists of finding \(B \in H_0({{\,\textrm{curl}\,}};\Omega ), p \in H^1_0(\Omega )\) and \(z \in \mathcal {H}^1\), such that

In our broken FEEC framework this problem is discretized as: find \(p_h \in V^0_h\), \(B_h \in V^1_h\) and \(z_h \in \mathcal {H}^1_h\), such that

for all \(q \in V^0_h\), \(v \in V^1_h\) and \(w \in \mathcal {H}^1_h\). Here, \(\alpha ^\ell \in \mathbb {R}\) are stabilization parameters for the jump terms in \(V^0_h\) and \(V^1_h\), and \(\mathcal {H}^1_h = \ker (\mathcal {L}^1_h + \alpha ^1S^1_h)\) may be computed as the kernel of the stabilized discrete Hodge–Laplace operator (47). As discussed above it coincides with that of the conforming Hodge–Laplace operator, hence

In operator form, this amounts to finding \(p_h \in V^0_h\), \(B_h \in V^1_h \cap (\mathcal {H}^1_h)^\perp \) and \(z_h \in \mathcal {H}^1_h\) such that

Proposition 5

System (87) [or equivalently (89)] is well-posed for arbitrary non-zero stabilization parameters \(\alpha ^0\), \(\alpha ^1 \ne 0\), and its solution satisfies

Moreover the discrete magnetic field can be decomposed according to (17), as

with discrete grad, harmonic and curl components characterized by

In particular, the solution is independent of the non-zero stabilization parameters and the conforming projection operators.

Proof

The result follows by using appropriate test functions and the orthogonality of the conforming discrete Hodge–Helmholtz decomposition (17). In particular, taking \(v = z_h\) (in \(\mathcal {H}^1_h \subset V^{1,c}_h\)) yields \(z_h=0\), taking \(q = (I-P^0_h) p_h\) and \(v = {{\,\textrm{grad}\,}}{\bar{P}}^0_h p_h\) shows that \(p_h = 0\), and the characterization of \(B_h\) follows by taking first \(v = (I-P^1_h) B_h\) and then \(v \in V^{1,c}_h\), an arbitrary conforming test function. \(\square \)

Remark 4

Proposition 5 states that the auxiliary unknowns \(p_h\) and \(z_h\) may be discarded a posteriori from the system. However they are needed in (87) to have a square matrix.

In practice the discrete harmonic fields can be represented by a basis of \(\mathcal {H}^1_h\), of the form \(\Lambda ^{1,H}_i\), \(i = 1, \ldots , \dim (\mathcal {H}^1_h) = b_2(\Omega )\), computed as the zero eigenmodes of a penalized \(\mathcal {L}^1_h + \alpha S^1_h\), where \(\alpha > 0\) is arbitrary and \(b_2(\Omega )\) is the Betti number of order 2, i.e. the number of cavities in \(\Omega \). Using these harmonic fields, one assembles the rectangular mass matrix

which allows us to rewrite (87) as

where \({\varvec{\textsf{p}}}\), \({\varvec{\textsf{B}}}\) and \({\varvec{\textsf{z}}}\) contain the coefficients of \(p_h\), \(B_h\) and \(z_h\) in the (primal) bases of \(V^0_h\), \(V^1_h\) and \(\mathcal {H}^1_h\) respectively.

3.5.2 Magnetostatic Problem with Metallic Boundary Conditions

We next consider a boundary condition of the form

which corresponds to a metallic (perfectly conducting) boundary. In our framework it can be approximated by using the inhomogeneous sequence (4) which we denote by \({\bar{V}}^\ell \) as in Sect. 3.3, and by modelling the boundary condition in a weak form. The mixed formulation takes a form similar to (84) with a few changes. Namely, we now look for \(p \in {\bar{V}}^0 = H^1(\Omega )\) and \(B \in {\bar{V}}^1 = H({{\,\textrm{curl}\,}};\Omega )\), solution to

Compared to (84), the only change is in the spaces chosen for the test functions. In particular we may formally decompose the second equation into a first set of test functions with \(q|_{\partial \Omega } = 0\) which enforces \({{\,\textrm{div}\,}}B = 0\), and a second set with \(q|_{\partial \Omega } \ne 0\) which enforces \(n \cdot B = 0\) on the boundary.

Again, this problem is not well-posed in general: on domains with holes (i.e., handles) a constraint must be added corresponding to harmonic forms, which leads to an augmented system with one additional Lagrange multiplier \(z_h\) as we did in the previous section. Here the proper space of harmonic forms is the one associated with the inhomogeneous sequence (4) involving \({\bar{V}}^1 = H({{\,\textrm{curl}\,}};\Omega )\) and \({\bar{V}}^*_1 = H_0({{\,\textrm{div}\,}};\Omega )\). Thus, it reads

We also observe that now p is only defined up to constants, so that an additional constraint needs to be added, such as \(p \in (\ker {{\,\textrm{grad}\,}})^\perp \). In practice this constraint may either be enforced in a similar way as the harmonic constraint, or by a so-called regularization technique where a term of the form \({\varepsilon }\langle q,p\rangle \), with \({\varepsilon }\ne 0\), is added to the first equation of (94). The constrained problem reads then: Find \(B \in H({{\,\textrm{curl}\,}};\Omega ), p \in H^1(\Omega )\) and \(z \in {\bar{\mathcal {H}}}^1\), such that

We note that here the result does not depend on \({\varepsilon }\) (indeed by testing with \(v = {{\,\textrm{grad}\,}}p\) we find that p is a constant, and with \(q = p\) we find that \(p = 0\)), hence we may set \({\varepsilon }= 1\).

At the discrete level we write a system similar to (87) or (89), but we must replace the space (88) with

where we have denoted \({\bar{V}}^{\ell ,c}_h = V^\ell _h \cap {\bar{V}}^\ell \). Observe that the for the inhomogeneous sequence, the decomposition (17) becomes

(98)

where \(\widetilde{{\,\textrm{curl}\,}}^c_h: {\bar{V}}^{2,c}_h \rightarrow {\bar{V}}^{1,c}_h\) is now the discrete adjoint of \({{\,\textrm{curl}\,}}^c_h: {\bar{V}}^{1,c}_h \rightarrow {\bar{V}}^{2,c}_h\).

We compute the space (97) as the kernel of the stabilized Hodge–Laplace operator associated with the inhomogeneous sequence, \({\bar{\mathcal {H}}}^1_h = \ker ({\bar{\mathcal {L}}}^1_h + \alpha ^1 {\bar{S}}^1_h)\), with an arbitrary positive \(\alpha ^1 > 0\). Note that the matrix of \({\bar{\mathcal {L}}}^1_h\) takes the same form as in (48),

where \({\bar{\mathbb {P}}}^0\) and \({\bar{\mathbb {P}}}^1\) are the conforming projection matrices on the inhomogeneous spaces, which leave the boundary degrees of freedom untouched instead of setting them to zero. The Hodge (mass) and differential matrices are those of the broken spaces which have no boundary conditions, hence they are the same as in Sect. 3.5.1. The resulting system reads: find \(p_h \in {\bar{V}}^0_h\), \(B_h \in {\bar{V}}^1_h\) and \(z_h \in {\bar{\mathcal {H}}}^1_h\), such that

for all \(q \in {\bar{V}}^0_h\), \(v \in {\bar{V}}^1_h\) and \(w \in {\bar{\mathcal {H}}}^1_h\). An operator form similar to (89) can also be written.

Finally the matrix form of our broken FEEC magnetostatic equation with metallic boundary conditions is similar to (93), with the following differences:

• the conforming projection operators map in the full conforming spaces \({\bar{V}}^\ell \) (boundary dofs are not set to zero),

• the discrete harmonic space is now \({\bar{\mathcal {H}}}^1_h\) the kernel of the \({\bar{\mathcal {L}}}^1_h\) operator associated with the full sequence, which in practice is also implemented by using conforming projection operators on the inhomogeneous spaces,

• a regularization term \(\mathbb {M}^0 {\varvec{\textsf{p}}}\) involving the mass matrix in \(V^0_h\) is added to the first equation to fully determine p.

Proposition 6

System (100) is well-posed for arbitrary stabilization parameters \(\alpha ^0 > 0\), \(\alpha ^1 \ne 0\), and its solution satisfies

Moreover the discrete magnetic field can be decomposed according to (98), as

with discrete grad, harmonic and curl components characterized by

In particular, the solution is independent of the parameters \(\alpha ^0 > 0\), \(\alpha ^1 \ne 0\), and the conforming projection operators.

Proof

Taking \(v = z_h\) (in \({\bar{\mathcal {H}}}^1_h \subset {\bar{V}}^{1,c}_h\)) and using the orthogonality (98) shows that \(z_h = 0\). Taking next \(v = {{\,\textrm{grad}\,}}{\bar{P}}^0_h p_h\) and using again (98) shows that \({\bar{P}}^0_h p_h\) is a constant, and with \(q = 1\) we see that \(\langle 1,p_h\rangle = 0\). It follows that with \(q = (I-{\bar{P}}^0_h) p_h\) the first equation becomes \(\Vert p_h\Vert ^2 + \alpha ^0 \Vert (I-{\bar{P}}^0_h)p_{h}\Vert ^2 = 0\), hence \(p_h = 0\). Taking next \(v = (I-{\bar{P}}^1_h) B_h\) yields \(\Vert (I-{\bar{P}}^1_h) B_h\Vert ^2 = 0\), hence \(B_h \in {\bar{V}}^{1,c}_h\). Since it is orthogonal to both \({\bar{\mathcal {H}}}^1_h\) and \({{\,\textrm{grad}\,}}{\bar{V}}^{0,c}_h\), it belongs to \(\widetilde{{\,\textrm{curl}\,}}^c_h V^{2,c}_h\), and it is characterized by \(\langle {{\,\textrm{curl}\,}}v,{{\,\textrm{curl}\,}}B_h\rangle = \langle {{\,\textrm{curl}\,}}v,J\rangle \) for all \(v \in {\bar{V}}^{1,c}_h\), which completes the proof. \(\square \)

3.6 Time-Dependent Maxwell’s Equations

We conclude this section by recalling the CONGA discretization of Maxwell’s time-dependent equations,

which consists of computing \(E_h \in V^1_h\) and \(B_h \in V^2_h\) solution to

with a discrete source approximated with the dual commuting projection

This discretization has been proposed and studied in [18,19,20], where it has been shown to be structure-preserving and have long-time stability properties: in addition to conserve energy in the absence of sources, it preserves exactly the discrete Gauss laws

thanks to the commuting diagram property \(\widetilde{{\,\textrm{div}\,}}_h {\tilde{\Pi }}^1_h J = {\tilde{\Pi }}^0_h {{\,\textrm{div}\,}}J\) satisfied by the dual projection operators, see (34). Moreover, the filtering by \((P^1_h)^*\) involved in the source approximation operator (32) avoids the accumulation of approximation errors in the non-conforming part of the kernel of the CONGA \({{\,\textrm{curl}\,}}{{\,\textrm{curl}\,}}\) operator [18].

Remark 5

The full commuting diagram property (34) is actually not needed for the preservation of the discrete Gauss law. One can indeed verify that (106) still holds for a current source defined by a simple \(L^2\) projection,

since one has in this case \(\langle \widetilde{{\,\textrm{div}\,}}_h J_h,{\varphi }\rangle = -\langle J_h,{{\,\textrm{grad}\,}}P^0_h {\varphi }\rangle = -\langle J,{{\,\textrm{grad}\,}}P^0_h {\varphi }\rangle = \langle {{\,\textrm{div}\,}}J,P^0_h{\varphi }\rangle = \langle {\tilde{\Pi }}^0_h {{\,\textrm{div}\,}}J, {\varphi }\rangle \) so that (106) follows from the compatibility relation \({{\,\textrm{div}\,}}J + {\partial _{t}}\rho = 0\) satisfied by the exact sources. The additional filtering of the source by \((P^1_h)^*\) involved in the dual commuting projection (105) is nevertheless important for long term stability, as analyzed in [18] and verified numerically in Sect. 5.7 below.

An attractive feature of the CONGA discretization is that for an explicit time-stepping method such as the standard leap-frog scheme

each time step is purely local in space. This is easily seen in the matrix form of (108)

where all the matrices except \(\mathbb {P}^1\) are patch-diagonal, and \(\mathbb {P}^1\) only couples degrees of freedom of neighboring patches.

Another key property is that for a time-averaged current source,

the discrete Gauss law (106) is also preserved by the fully discrete solution, namely

holds for all \(n > 0\), provided it holds for \(n=0\). We also remind that in the absence of source (\(J=0\)), the leap-frog time scheme requires a standard CFL condition for numerical stability, of the form

(112)

where . Indeed the discrete pseudo-energy

is a constant, \(\mathcal {H}^{n,*}_h = \mathcal {H}^{*}_h\), so that (112) yields a long-time bound

valid for all \(n \ge 0\). In practice we set with \(C_{\text {cfl}} \approx 0.8\), and we evaluate the operator norm of \({{\,\textrm{curl}\,}}_h\) with an iterative power method, noting that it amounts to computing the spectral radius of a diagonalizable matrix with non-negative eigenvalues, namely

In this section we detail the construction of a broken FEEC discretization on a 2D multipatch domain where each patch is the image of the reference square,

with smooth diffeomorphisms \(F_k\), \(k = 1, \dots , K\). For simplicity we consider the case of planar domains \(\Omega _k \subset \mathbb {R}^2\) and orientation-preserving mappings, in the sense that the Jacobian matrices \(DF_k({\hat{{\varvec{x}}}})\) have positive Jacobian determinants \(J_{F_k}({\hat{{\varvec{x}}}}) = \det (DF_k({\hat{{\varvec{x}}}})) > 0\) for all \({\hat{{\varvec{x}}}} \in {\hat{\Omega }}\).

We also assume that the multipatch decomposition is conforming in the sense that any interface \(\Gamma _{k,l} = \partial \Omega _k \cap \partial \Omega _l\) between two distinct patches

1. (i)

   is either a vertex, or a full edge of both patches,

2. (ii)

   admits the same parametrization from both patches, up to the orientation.

In the case where the interface is a vertex, condition (ii) is empty. If it is a full patch edge of the form

where \({\hat{{\varvec{x}}}}_0\) and \({\hat{{\varvec{y}}}}_0\) are vertices of the reference domain \({\hat{\Omega }}\), condition (ii) means that

holds with \(\theta _{k,l}\) an affine bijection on [0, 1], that is \(\theta _{k,l}(s) =s\) or \(1-s\).

With spline spaces defined on each patch with symmetric knot sequences, these conditions imply that the patches are fully matching in the sense of Assumption 3.3 from Buffa et al. [13].

Examples of such domains are represented on Fig. 1 where the left one is a standard curved L-shaped domain with circular patch boundaries that is used in some reference academic test-cases, see e.g. [29], and the right one is a non contractible domain that also involves analytical mappings, shaped as a simplified pretzel for cultural reasons. In the latter domain we remind that the presence of holes leads a priori to non exact de Rham sequences, and hence non trivial harmonic forms.

Two multipatch domains: a curved L-shaped (left) and a non-contractible domain (right). Solid and dashed lines are used to represent the domain boundaries and the patch interfaces, respectively

Full size image

On such multipatch domains, the construction of mapped spline complexes is fairly standard

[15, 16]: for each patch a reference spline complex is first defined on the logical domain \({\hat{\Omega }}\) using tensor-product spaces, and the mapped spline spaces are next obtained by transforming the reference ones with push-forward operators. To build our conforming projection operators we shall use geometric degrees of freedom corresponding to pointwise evaluations and integrations over geometric elements, following the approach of Bochev et al. [4], Robidoux [47] and Gerritsma [31]. Specifically, we consider broken geometric degrees of freedom of the form

associated with local grids with nodes \({\texttt {n}}_{k,{\varvec{i}}}\), edges \({\texttt {e}}_{k,\mu }\) and cells \({\texttt {c}}_{k,{\varvec{i}}}\) in each patch \(\Omega _k\). Here \({\varvec{\tau }}^* = {\varvec{\tau }}\) (a unit vector tangent to \({\texttt {e}}_{k,\mu }\)) or \({\varvec{\tau }}^\perp \) (normal to the edge), depending on which 2D de Rham sequence is used. We refer to Section B in the Appendix for more details.

4.1 Conforming Projection Operators

In addition to yield commuting projection operators, geometric degrees of freedom naturally provide a simple characterization of the discrete fields \(v_h\) in the broken spaces \(V^\ell _h\) which are actually conforming. In order to belong to the space \(V^\ell \) appearing in the continuous de Rham sequence (1) (and hence to the conforming spline space \(V^{\ell ,c}_h = V^{\ell }_h \cap V^\ell \)) a field must satisfy regularity conditions which are well known for piecewise-smooth functions. For \(\ell = 0\), the condition for \(H^1(\Omega )\) regularity amounts to continuity across the patch interfaces, which is equivalent to requiring that the broken degrees of freedom associated to the same interpolation node coincide:

For \(\ell = 1\) in the \(\textbf{grad}-{{\,\textrm{curl}\,}}\) sequence, the interface constraint for \({\varvec{H}}({{\,\textrm{curl}\,}};\Omega )\) regularity consists in the continuity of the tangential traces, which amounts to requiring that edge degrees of freedom associated to the same edge coincide, up to the orientation. Namely,

where \({\varepsilon }^{\texttt {e}}(k,\mu ;l,\nu ) = \pm 1\) is the relative orientation of the edges, see (B.15).

For \(\ell = 1\) in the \(\textbf{curl}-{{\,\textrm{div}\,}}\) sequence, the interface constraint for \({\varvec{H}}({{\,\textrm{div}\,}};\Omega )\) regularity consists in the continuity of the normal traces, which again amounts to requiring that edge degrees of freedom associated to the same edge coincide, up to the orientation: this constraint takes the same form as (118). Finally as \(V^2 = L^2(\Omega )\) there are no constraints for \(\ell =2\), i.e., the spaces \(V^{2}_h\) and \(V^{2,c}_h\) coincide. We gather these constraints in a single formula

where \({\texttt {g}}^\ell _{k,\mu }\) denotes the geometric element (node, edge or cell) of dimension \(\ell \) associated with a multi-index \((k,\mu ) \in \mathcal {M}^\ell _h\), and \({\varepsilon }^1\) denotes the relative orientation of two edges, while \({\varepsilon }^0 = {\varepsilon }^2:= 1\).

Thanks to the conformity assumption of the multipatch geometry and to the symmetry of the interpolatory grid, each broken degree of freedom on an interface can be matched to those of the adjacent patches in such a way that the constraints above are satisfied. Accordingly, the resulting broken discrete field will belong to the conforming space \(V^{\ell ,c}_h\). In particular we may define simple conforming projection operators \(P^\ell \) by averaging the broken degrees of freedom associated with interface elements. Using the broken basis functions \(\Lambda ^{\ell }_{k,\mu }\) associated with the above degrees of freedom, this yields an expression similar to the one given in [23] for tensor-product polynomial elements, with additional relative orientation factors due to the general mapping configurations,

where \(\mathcal {M}^\ell _h({\texttt {g}}^\ell _{l,\nu })\) contains the patch-wise indices corresponding to a given geometric element. The entries of the corresponding operator matrix (36) read

for all \((k,\mu ), (l,\nu ) \in \mathcal {M}^\ell _h\). We may summarize our construction with the following result.

Proposition 7

The broken geometric degrees of freedom (116) satisfy the properties listed in Sect. 2.4, with local domain spaces \(U^\ell (\Omega _k)\) defined in (B.19) and local differential matrices independent of the mappings \(F_k\). Moreover, the operators \(P^\ell _h: V^\ell _h \rightarrow V^{\ell }_h\) defined by (120) are projections on the conforming subspaces \(V^{\ell ,c}_h = V^{\ell }_h \cap V^{\ell }\).

Proof

The boundedness of \(\sigma ^\ell _{k,\mu }\) on the local spaces \(U^\ell (\Omega _k)\), as well as the commuting properties, are verified in Proposition B.2 in the Appendix. The inter-patch conformity property (31) follows from the fact that piecewise smooth functions that belong to \(H^1(\Omega )\), resp. \(H({{\,\textrm{curl}\,}};\Omega )\) and \(H({{\,\textrm{div}\,}};\Omega )\), admit a unique trace, resp. tangential and normal trace on patch interfaces [5]. The same property, together with the interpolation nature of the geometric basis functions, allows verifying that the operators \(P^\ell _h: V^\ell _h \rightarrow V^\ell _h\) are characterized by the relations

for all \((k,\mu ) \in \mathcal {M}^\ell _h\). Using the geometric condition (119), this allows verifying that \(P^\ell _h\) is indeed a projection on the conforming subspace \(V^{\ell ,c}_h\). \(\square \)

Remark 6

(Boundary conditions) As stated, the above construction actually corresponds to the inhomogeneous sequences. If one considers the spaces with homogeneous boundary conditions,

then \(P^\ell _h\) should further set the boundary degree of freedom to 0, which amounts to restricting the non-zeros entries in (121) to the geometrical elements \({\texttt {g}}^\ell _{k,\mu }\) that are inside \(\Omega \). As already observed in Sect. 3.3, this does not affect the broken spaces \(V^\ell _h\) since they are not required to have boundary conditions.

Remark 7

(Extension to 3D) The extension to the 3D setting of the above construction poses no particular difficulty, using the same tensor-product and mapped spline spaces as in [15, 27], and the same geometric degrees of freedom as in [24, Sect. 6.1] for the primal sequence.

4.2 Primal-Dual Matrix Diagram with B-Splines

In practice, a natural approach is to work with B-splines. Indeed the geometric (interpolatory) splines \(\Lambda ^\ell _{k,\mu }\) (dual to the geometric degrees of freedom) are not known explicitly and depend on the interpolation grids. They are also less local than the B-splines, as they are in general supported on their full patch \(\Omega _k\): for patches with many cells, this would lead to an expensive computation of the fully populated mass matrices.

In our numerical experiments we have followed da Veiga et al. [27] and Campos Pinto et al. [24] and used tensor-product splines \({\hat{B}}^\ell _{\mu }\) on the reference domain \({\hat{\Omega }}\), composed of normalized B-splines \(\mathcal {N}^p_i\) in the dimensions of degree p and of Curry-Schoenberg B-splines \(\mathcal {D}^{p-1}_{i} = \Big (\frac{p}{\xi _{i+p+1}-\xi _{i+1}}\Big ) \mathcal {N}^{p-1}_{i+1}\) (also called M-splines) in the dimensions of degree \((p-1)\). On the mapped patches \(\Omega _k = F_k({\hat{\Omega }})\) the basis functions are defined again as push-forwards \(B^\ell _{k,\mu }:= \mathcal {F}^\ell _k({\hat{B}}^\ell _\mu )\) for \((k,\mu ) \in \mathcal {M}^\ell _h\).

The discrete elements in the matrix diagram (41) must then be adapted for B-splines: one first observe that the change of basis is provided by the patch-wise collocation matrices whose diagonal blocks read

where we have used (B.10) from the Appendix. As \(B^\ell _{k,\nu } = \sum _\mu \mathbb {K}^\ell _{(k,\mu ), (k,\nu )} \Lambda ^\ell _{k,\mu }\) it follows that the B-spline coefficients of the geometric projection (30), namely

read (using some implicit flattening \(\mathcal {M}^\ell _h \ni (k,\mu ) \mapsto i \in \{1, \dots , N^\ell \}\))

Accordingly, we obtain the matrices \(\mathbb {P}^\ell _B = (\beta ^\ell _i(P^\ell _h B^\ell _j))_{1 \le i,j \le N^\ell }\) of the conforming projection operator (120) in the B-spline bases by combining the matrices (121) with the above change of basis: this gives

Here we note that the collocation matrices \(\mathbb {K}^\ell \) are Kronecker products of univariate banded matrices, see [24, Sect. 6.2], so that (125) and (126) may be implemented in a very efficient way. The primal sequence is completed by computing the patch-wise differential matrices in the B-spline bases,

Thanks to the univariate relation \(\frac{\, \textrm{d}}{\, \textrm{d}x}\mathcal {N}^p_i = \mathcal {D}^{p-1}_{i-1} - \mathcal {D}^{p-1}_{i}\), these are patch-wise incidence matrices, just as the ones in the geometric bases [24, 27]. For the dual sequence we use bi-orthogonal splines characterized by the relations

with dual degrees of freedom

This leads to defining again the primal Hodge matrices as patch-diagonal mass matrices, \(\mathbb {H}^\ell _B = \mathbb {M}^\ell _B\) and the dual Hodge ones as their patch-diagonal inverses \({\tilde{\mathbb {H}}}^\ell _B = (\mathbb {M}^\ell _B)^{-1}\). These matrices can be computed on the reference domain according to

Specifically, for the 2D \(\textbf{grad}-{{\,\textrm{curl}\,}}\) sequence we obtain

and for the 2D \(\textbf{curl}-{{\,\textrm{div}\,}}\) sequence we find

while \(\mathbb {M}^0_B\) and \(\mathbb {M}^2_B\) take the same values as in (129). Here we have used the explicit form of the pull-back operators (A.12) and (A.13) recalled in the Appendix. In 3D the formulas can be derived from the pull-backs given in [27].

Finally the coefficients of the dual commuting projections (32)–(33) in the dual B-spline bases read

using the dual degrees of freedom (127) and the same implicit flattening as in (125), (126).

Gathering the above elements in 3D we obtain a new version of diagram (41), where the coefficient spaces (still defined as \(\mathbb {R}^{N^\ell }\)) are now denoted \(\mathcal {C}^\ell _B\) and \({\tilde{\mathcal {C}}}^\ell _B\) to indicate that they contain coefficients in the B-spline and dual B-spline basis, respectively.

(131)

Note that here the primal and dual interpolation operators are simply

In this section we conduct numerical experiments to test our novel broken FEEC schemes from Sect. 3 using the multipatch spline spaces described in Sect. 4. These experiments have been performed with the Psydac library [32], where the different operators from the primal/dual commuting diagram (131) have been implemented. In contrast to conforming FEEC appproximations, we remind that these operators are local in the sense explained in the introduction, thanks ot the use of the fully discontinuous spaces \(V^\ell _h\) and the local conforming projections represented by the matrices \(\mathbb {P}^\ell _B\).

5.1 Poisson Problem with Homogeneous Boundary Conditions

We first test our CONGA scheme for the homogeneous Poisson problem presented in Sect. 3.1. For this we consider an analytical solution on the pretzel-shaped domain shown in Fig. 1, given by

Here, \(s({\varvec{x}}) = {\tilde{x}} -{\tilde{y}}\), \(t({\varvec{x}}) = {\tilde{x}} +{\tilde{y}}\) with \({\tilde{x}} = x- x_0\), \({\tilde{y}}= y-y_0\) and we take \(x_0 = y_0 = 1.5\), \(a = (1/1.7)^2\), \(b = (1/1.1)^2\) and \(\sigma = 0.11\) in order to satisfy the homogeneous boundary condition with accuracy \(\approx 1e-10\). The associated manufactured source is

In Table 1 we show the relative \(L^2\) errors corresponding to different grids and spline degrees, and in Fig. 2 we plot the numerical solutions corresponding to spline elements of degree \(3 \times 3\) on each patch. These results show that the numerical solutions converge towards the exact one as the grids are refined. We do not observe significant improvements however when higher order polynomials are used. As the next results will show, this is most likely due to the steep nature of the solution.

Solutions \(\phi _h \in V^0_h\) computed by the scheme (52), with spline elements of degree \(3 \times 3\) and \(N \times N\) of cells per patches as indicated, corresponding to some of the errors shown in Table 1

Full size image

5.2 Poisson Problem with Inhomogeneous Boundary Conditions

Our second test is with an inhomogeneous Poisson–Dirichlet problem (64), using the lifting method described in Sect. 3.3 for the boundary condition and the solver tested just above for the homogeneous part of the solution. Specifically, we define \(\phi _{g,h} \in V^0_h\) by computing its boundary degrees of freedom from the data g on \(\partial \Omega \) (this is straightforward with our geometric boundary degrees of freedom (B.12)), and we compute \(\phi _{0,h} = \phi _{h}-\phi _{g,h}\) by solving (67). As we are not constrained by a specific condition on the domain boundary we consider a smooth solution to assess whether high order convergence rates can be observed despite the singularities in the domain boundaries. Specifically, we use again the pretzel-shaped domain and set the source and boundary condition as

where \(\phi ({\varvec{x}}) = \sin (\pi x)\cos (\pi y)\) is the exact solution.

In Fig. 3 we plot the convergence curves corresponding to various \(N \times N\) grids for the 18 patches of the domain, and various degrees \(p = 2, \dots 5\). They show that the solutions converge with optimal rate \(p+1\) (and even \(p+2\) for \(p=2\)) as the patch grids are refined. Since the stabilized broken FEEC solution coincides with the conforming one according to our analysis in Sects. 3.1 and 3.3, these rates follow from the optimality of conforming FEM approximations to elliptic problems and the approximation power of multipatch spline spaces established in [16].

Convergence study for the inhomogeneous Poisson solver with source and boundary conditions given by (135). The left plot shows the relative \(L^2\) errors corresponding to various grids of \(N \times N\) cells per patch (i.e., \(18 N^2\) cells in total) and spline degrees \(p \times p\) as indicated. The right plot shows one numerical solution \(\phi _h\) of good accuracy

Full size image

5.3 Time-Harmonic Maxwell Problem with Homogeneous Boundary Conditions

We next turn to a numerical assessment of our CONGA solvers for the Maxwell equation, and as we did for the Poisson equation we begin with the homogeneous case presented in Sect. 3.2. Since now the solution depends a priori on the time pulsation \(\omega \) we opt for a physically relevant current source localized around the upper right hole of the pretzel-shaped domain,

where \(\phi \) and \(\tau \) are as in (133). We plot this source in Fig. 4.

Source (136) (vector field and amplitude) for the homogeneous Maxwell solutions plotted in Fig. 5 below

Full size image

Reference numerical solutions for the homogeneous Maxwell problem with source (136) and time pulsation \(\omega = \sqrt{50}\) (top) and \(\omega = \sqrt{170}\) (bottom). These solutions have been obtained using \(20 \times 20\) cells per patch and spline elements of degree \(6 \times 6\). The vector fields are shown on the left while the amplitudes are shown on the right

Full size image

Numerical solutions obtained with spline elements of degree \(3 \times 3\) and different grids as indicated, corresponding to the errors shown in Table 2

Full size image

We then consider two values for the time pulsation, namely \(\omega =\sqrt{50}\) and \(\omega =\sqrt{170}\), and for each of these values we use as reference solutions the numerical solutions computed using our method on a mesh with \(20 \times 20\) cells per patch (i.e. 7200 cells in total) and elements of degree \(6 \times 6\) in each patch. These reference solutions are shown in Fig. 5. Interestingly, we observe that for the higher pulsation \(\omega =\sqrt{170}\) the source triggers a time-harmonic field localized around the upper left hole, opposite to where the source is.

In Table 2 we show the \(L^2\) errors corresponding to different grids and spline degrees for the two values of the pulsation \(\omega \), and we also plot in Fig. 6 the different solutions corresponding to spline elements of degree \(3 \times 3\). These results show that the numerical solutions converge towards the reference ones as the grids are refined, and with a faster convergence in the case of the lower pulsation \(\omega = \sqrt{50}\), due to the higher smoothness of the corresponding solution.

5.4 Time-Harmonic Maxwell Problem with Inhomogeneous Boundary Conditions

As we did for the Poisson problem, we next test our Maxwell solver with a smooth solution and handle the inhomogeneous boundary conditions with the lifting method described in Sect. 3.3. Specifically, we define \({\varvec{u}}_{g,h} \in V^1_h\) by computing its boundary degrees of freedom from the data \(g = {\varvec{n}}\times {\varvec{u}}\) on \(\partial \Omega \) [this is again straightforward with our geometric boundary degrees of freedom (B.13)], and we compute \({\varvec{u}}_{0,h} = {\varvec{u}}_{h}-{\varvec{u}}_{g,h}\) by solving (69). Here we take \(\omega = \pi \) and consider (65) with the source-solution pair

and a boundary condition given by \(g:= {\varvec{n}}\times {\varvec{u}}\) on \(\partial \Omega \).

In Fig. 7 we plot the convergence curves corresponding to spline elements of various degrees \(p\times p\) and \(N \times N\) cells per patch. The results are similar to what was observed for the Poisson problem: the solutions converge with optimal rate \(p+1\) as the grids are refined (again a rate of \(p+2\) is observed for \(=2\)) which confirms the numerical accuracy of our approach for the Maxwell problem combined with a geometric lifting technique for the Dirichlet boundary condition. These convergence rates can be justified by the approximation power of multipatch spline spaces established in [16] and the error analysis from Monk and Demkowicz [43] which establishes optimal convergence under the condition that the finite element spaces satisfy the discrete compactness property: for FEEC spaces the discrete compactness property follows from the existence of \(L^2\) stable commuting projection operators (see Sect. 3.6 [2]), and for multipatch spaces this has been recently proven in [25]. Note that the convergence observed here is actually better than expected, since \(H({{\,\textrm{curl}\,}})\) multipatch spline spaces of maximal degree p should only have p order accuracy [16].

Convergence study for the time-harmonic Maxwell solver with source given by (137) and inhomogeneous boundary conditions \({\varvec{n}}\times {\varvec{u}}= g\). Relative \(L^2\) errors are plotted on the left for various grids of \(N \times N\) cells per patch (corresponding to a total of \(18 N^2\) cells) and spline degrees \(p \times p\) as indicated. The right plot shows the amplitude of a numerical solution \({\varvec{u}}_h\) of good accuracy

Full size image

5.5 Eigenvalue Problems

We next assess the accuracy of our CONGA approximation (76) for the curl-curl eigenvalue problem.

We test our discretization on the two domains shown in Fig. 1, and plot in Fig. 8 the amplitude of the first five eigenmodes, together with their positive eigenvalues. Here the eigenmodes are computed using spline elements of degree \(6 \times 6\) and \(56 \times 56\), resp. \(20 \times 20\) cells per patch in the case of the curved L-shaped, resp. pretzel-shaped domain composed of 3, resp. 18 patches. On the former domain this corresponds to 9408 cells in total and 22,692 degrees of freedom for the broken space \(V^1_h\), while on the latter domain it corresponds to 7200 cells and 23,400 degrees of freedom for \(V^1_h\) (a higher value than for the L-shaped domain despite less cells, because of the duplication of boundary dofs at the patch interfaces).

In Fig. 9 we then plot the relative eigenvalue errors

as a function of the eigenvalue index i, for degrees \(p = 3\) and 5, and \(N \times N\) cells per patch with \(N = 2, 4, 8\) and 16. For the curved L-shaped domain we use as reference the eigenvalues provided as benchmark in [29, 30], and for the pretzel-shaped domain we use the eigenvalues computed using our CONGA scheme, with as many reliable digits as we could find using uniform patches with degree \(6 \times 6\) and \(N \times N\) cells per patch with \(N \le 20\) (this limit being imposed by the fact that we compute the matrix eigenmodes with Scipy’s eigsh solver with a sparse LU decomposition).

This allows us to verify numerically that the discrete eigenvalues converge towards the exact ones, with smaller errors corresponding to the smoother eigenmodes visible in Fig. 8.

Eigenmodes of the curl-curl problem on the curved-L-shaped domain (left) and on the pretzel domain (right), obtained on fine grids with \(20 \times 20\) cells per patch and degree \(6 \times 6\)

Full size image

Relative eigenvalue errors (138) for the curl-curl problem on the curved-L-shaped domain (left) and on the pretzel domain (right), using elements of degree \(p \times p\) as indicated, and \(N \times N\) cells per patch with \(N = 2, 4, 8\) and 16. The resulting total numbers of cells are \(\texttt {nc} = 3N^2\) and \(18N^2\) for the respective multipatch domains

Full size image

5.6 Magnetostatic Test-Cases

We next study the CONGA discretizations of the magnetostatic problems presented in Sect. 3.5: either the one for the problem with pseudo-vacuum boundary conditions (86) or the one with metallic boundary conditions (96).

To this end we consider a scalar dipole current source,

in the pretzel-shaped domain, with \(\sigma = 0.02\). We set the positive current pole at \(x_0 = y_0 = 1\) and the negative one at \(x_1 = y_1 = 2\). We then consider the discrete solvers described in Sects. 3.5.1 and 3.5.2 for the problems with pseudo-vacuum and metallic boundary conditions, respectively.

In Fig. 10 we plot the scalar source \(J_z\) together with the vector-valued \(\textbf{curl}J_z\) field, and for each of the boundary conditions, we plot in Fig. 11 fine solutions computed using the CONGA scheme on a mesh with \(20 \times 20\) cells per patch (i.e. 7200 cells in total) and elements of degree \(6 \times 6\) in each patch.

In Table 3 we then show the relative \(L^2\) errors corresponding to coarser grids and lower spline degrees for both boundary conditions, using as reference the fine solutions shown in Fig. 11. We also plot in Fig. 12 the solutions corresponding to spline elements of degree \(3 \times 3\) on each patch. Again these results indicate that our CONGA solutions converge nicely as the grids are refined, with smaller errors associated with higher polynomial degrees.

Source for the magnetostatic test-cases: the scalar current density (139) is plotted on the left panel and its vector-valued curl is shown on the right

Full size image

Reference numerical solutions obtained with \(20 \times 20\) cells per patch and elements of degree \(6 \times 6\), for the magnetostatic test-cases with pseudo-vacuum (top) and metallic boundary conditions (bottom). The vector fields are shown on the left while the amplitudes are shown on the right

Full size image

Numerical solutions for the magnetostatic test-cases with pseudo-vacuum (top) and metallic boundary conditions (bottom), obtained with elements of degree \(3 \times 3\) and different patch-wise grids, as indicated

Full size image

5.7 Time-Dependent Maxwell Equation

We finally assess the qualitative and quantitative properties of our mapped spline-based CONGA scheme for the time-dependent Maxwell system, using a leap-frog time stepping (108)–(109). We will consider two test-cases.

Our first test-case consists of an initial electric pulse

in the pretzel-shaped domain, with \(x_0 = y_0 = 1\) and \(\sigma = 0.02\). The initial magnetic field and the source is

In Fig. 13 we compare successive snapshots of two solutions corresponding to different numbers \(N \times N\) of cells per patch and degrees \(p \times p\): on the left plots the solution is computed with \(N=8\) and \(p=3\), while that on the right plots uses \(N=20\) and \(p=6\). We observe that these results display a qualitatively correct behaviour for propagating waves with reflecting boundary conditions, and that the profile for both resolutions are similar up to small scale features, which can be seen as a practical indicator of convergence.

Since this test-case is without source the discrete energy should be preserved up to bounded oscillations, as recalled in Sect. 3.6. In Fig. 14 we plot the time evolution of the discrete energy \(\mathcal {H}_h(t^n) = \frac{1}{2} \big (\Vert {\varvec{E}}^n_h\Vert ^2 + \Vert B^{n}_h\Vert ^2\big )\), together with that of the electric and magnetic fields, and we observe that the total energy is very stable. We also plot the amplitude of the discrete divergence as time evolves: assuming all computations exact it should be zero according to (106), indeed we have \(\rho = {{\,\textrm{div}\,}}{\varvec{E}}(t=0) = 0\) in this test case. On Fig. 14 we verify that it is zero up to machine accuracy and small quadrature errors.

Snapshots of numerical solutions to the time-dependent Maxwell test-case (140)–(141) at \(t = 0\), 0.4, 0.8, 1.6 and 3.2 (from top to bottom) discretized with the CONGA method (104) using spline elements of degree \(p \times p\) and \(N \times N\) cells per patch, and a leap-frog time stepping. The solution on the left corresponds to \(N = 8\) and \(p=3\), while the one on the right has been obtained \(N = 20\) and \(p=6\)

Full size image

Energy (top) and discrete divergence amplitude (bottom) of the coarse (left) and fine (right) numerical solutions of the time-dependent Maxwell test-case (140)–(141) shown in Fig. 13

Full size image

We next study the quality of the source approximation operator \({\varvec{J}}\rightarrow {\varvec{J}}_h\) which is involved in the discrete Ampere equation (104). For this we use a second time-dependent test-case with a zero initial condition

and a source of the form

with \(\psi \) as in (140). The associated charge density is then

In Figs. 15 and 16 we compare three different approximation operators for the current source, namely:

1. (i)

   the primal finite element projection: \({\varvec{J}}_h = \Pi ^1_h {\varvec{J}}\),

2. (ii)

   the \(L^2\) projection on the broken space: \({\varvec{J}}_h = {Q_{V^1_h}} {\varvec{J}}\),

3. (iii)

   the dual projection: \({\varvec{J}}_h = {\tilde{\Pi }}^1_h {\varvec{J}}\)

where we remind that \({\tilde{\Pi }}^1_h = (P^1_h)^* Q_{V^1_h}\), see (32). We note that each of these projection operators are local, in the sense that none requires solving a global problem on the computational domain \(\Omega \). We also remind that the primal projection \(\Pi ^1_h\) interpolates the geometric (edge) degrees of freedom and satisfies a commuting diagram property with the primal (strong) differential operators, but not with the dual ones. Hence it does not allow to preserve the discrete electric Gauss law in (106). In contrast, both the \(L^2\) projection on the broken space \(V^1_h\) and the dual projection \({\tilde{\Pi }}^1_h\) guarantee the preservation of the discrete Gauss laws, however we expect an increased stability for the latter one as discussed in Remark 5.

In Fig. 15 we first show some snapshots of the three numerical solutions on a time range \(t \in [0,T]\) with \(T=20\): There we see that the primal projection yields a very strong and steadily growing field in the region of the source (visible from the changing color scale) which points towards a large error. In contrast, the \(L^2\) and the dual projections produce solutions with moderate amplitude with some visible differences, namely an electric field that also builds up in the region where the source is located.

To better analyse the quality of these simulations we next show in Fig. 16 two error indicators for each one of the numerical strategies described above, namely the Gauss law errors associated with the broken solution \(\mathcal {G}^n_h({\varvec{E}}^n_h)\), and that of its conforming projection \(\mathcal {G}^n_h(P^1_h {\varvec{E}}^n_h)\), where we have denoted

Here the former errors are expected to be zero for both the \(L^2\) and dual projections, see Remark 5: The numerical results confirm this value, up to machine accuracy, whereas they show significant errors for the primal projection shown in the left plots. As for the second error, it has no reason to be strictly zero but should remain small for accurate and stable solutions, hence it is also an interesting error indicator. Here the curves show very large values (around 80 at \(t=T\)) for the primal projection operator \(\Pi ^1_h\), with a linear time growth (which is somehow consistent with the strong growth of the former error indicator \(\mathcal {G}^n_h({\varvec{E}}^n_h)\)). For the \(L^2\) projection the error is smaller but it is far from being negligible (close to 3), and also grows linearly in time. In contrast, the error is much smaller (on the order of 0.01) for the dual projection \({\tilde{\Pi }}^1_h\), and it oscillates but does not seem to grow. These results tend to indicate that the growing field visible in Fig. 15 corresponds to a numerical error, and they highlight the improved stability of the CONGA scheme with a proper source approximation.

Snapshots of numerical solutions to the time-dependent Maxwell problem with a time-varying source (143) discretized with the CONGA scheme (104) using spline elements of degree \(3 \times 3\) and \(8 \times 8\) cells per patch. The amplitudes \(|P^1_h{\varvec{E}}^n_h|\) are shown at \(t = 2\), 5, 10 and 20 (from top to bottom). In the left panels the source projection \({\varvec{J}}\rightarrow {\varvec{J}}_h\) is a primal projection operator \(\Pi ^1_h\) which does not commute with the dual differential operators (note the time-varying color scale). The middle panels use an orthogonal projection \({Q_{V^1_h}}\), and the ones on the right a dual commuting projection \({\tilde{\Pi }}^1_h\). See the text for more details

Full size image

Discrete Gauss law errors (144) as a function of time for the test-case (142)–(143), using different approximation operators for the discrete source \({\varvec{J}}_h = P_J {\varvec{J}}\). The top panels show the errors relative to the broken field \({\varvec{E}}_h\), whereas those on the bottom panels show the error relative to its conforming projection \(P^1_h {\varvec{E}}_h\). For each indicator, the left, middle and right plots correspond to the different source approximation operators shown in Fig. 15

Full size image

In this work we have extended the classical theory of finite element exterior calculus (FEEC) to mapped multipatch domains, using finite element spaces that are fully discontinuous across the patch interfaces. We refer to this approach as the “broken FEEC” or “CONGA” (COnforming/Non conforming GAlerkin) method. While the foundational theory relative to the solution of the Hodge–Laplace equation was presented in recent work [23], here we have presented stable broken FEEC formulations for many problems arising in electromagnetic applications, including Poisson’s equation, time-dependent and time-harmonic (source and eigenvalue) Maxwell’s problems, and magnetostatic problems with pseudo-vacuum and metallic boundary conditions. Further, we have detailed a numerical framework based on tensor-product splines on each patch, under the assumption of geometric conformity across the patch interfaces.

For all the electromagnetic problems presented, we have verified our broken FEEC framework through extensive numerical testing in L-shaped and pretzel-like two-dimensional domains. The latter geometry is particularly challenging because of its three holes and sharp reentrant corner. The nominal order of accuracy was achieved in all cases, and the structure-preserving properties (such as divergence of harmonic constraints) were respected to floating-point accuracy. For the time-dependent Maxwell problem with a current source, we could also observe long-term stability of the method, and presented alternative formulations which lack this property.

As explained in the introduction, the main benefit of our approach compared to a standard conforming FEEC discretization is to provide local coderivative operators, as well as local Hodge and commuting projection operators for both the primal and dual complexes. As such, broken FEEC methods allow for parallel algorithms where only communications between contiguous patches are needed in the direct application of every operator appearing in the primal/dual commuting diagram.

Another appealing feature of this approach lies in the modularity of the implementation: as every operator in the commuting diagram involves either a single patch operator or a conforming projection (which performs averages between contiguous patches and can be applied matrix-free), extending a single-patch code to multipatch domains essentially amounts to writing these conforming projections operators. In particular, solving problems in the broken multipatch spaces can be done with a global data structure that just consists of a list of the patch-wise data structures.

Given its solid theoretical bases and convincing numerical results, we are confident that the broken FEEC framework will find practical use in the computational physics community. To this end we plan to relax the grid conformity constraints at the interfaces, allowing for independent refinement of the patches, and to investigate the efficiency of parallel implementations for high-performance computing applications. This will allow us to tackle large problems in three dimensions, including MHD and kinetic models for plasma physics.