A Broken FEEC Framework for Electromagnetic Problems on Mapped Multipatch Domains

We present a framework for the structure-preserving approximation of partial differential equations on mapped multipatch domains, extending the classical theory of finite element exterior calculus (FEEC) to discrete de Rham sequences which are broken, i.e., fully discontinuous across the patch interfaces. Following the Conforming/Nonconforming Galerkin (CONGA) schemes developed in Campos Pinto and Sonnendrücker (Math Comput 85:2651–2685, 2016) and Campos Pinto and Güçlü (Broken-FEEC discretizations and Hodge Laplace problems. arXiv:2109.02553, 2022), our approach is based on: (i) the identification of a conforming discrete de Rham sequence with stable commuting projection operators, (ii) the relaxation of the continuity constraints between patches, and (iii) the construction of conforming projections mapping back to the conforming subspaces, allowing to define discrete differentials on the broken sequence. This framework combines the advantages of conforming FEEC discretizations (e.g. commuting projections, discrete duality and Hodge–Helmholtz decompositions) with the data locality and implementation simplicity of interior penalty methods for discontinuous Galerkin discretizations. We apply it to several initial- and boundary-value problems, as well as eigenvalue problems arising in electromagnetics. In each case our formulations are shown to be well posed thanks to an appropriate stabilization of the jumps across the interfaces, and the solutions are extremely robust with respect to the stabilization parameter. Finally we describe a construction using tensor-product splines on mapped cartesian patches, and we detail the associated matrix operators. Our numerical experiments confirm the accuracy and stability of this discrete framework, and they allow us to verify that expected structure-preserving properties such as divergence or harmonic constraints are respected to floating-point accuracy.

1. Introduction. Thanks to enlightening research conducted over the last few decades [13,36,33,31,1,7,16], 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 [48,11] to high order curved elements in isogeometric analysis [17,4], via edge Nédélec elements [41,10]. 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 [12,8,27,40]. 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,42,23], or the Hamiltonian structure of the MHD and Vlasov-Maxwell equations [38,34,22]. 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 framework of dual complexes this leads to global discrete Hodge operators mapping the dual de Rham sequence to the primal one, as well as to global commuting projection operators on the dual sequence, as the latter also rely on L 2 projection operators on the finite element spaces.
In this article we follow the broken FEEC approach [20] first developed for Conforming/Nonconforming Galerkin (CONGA) schemes in [24,25]. The principle is to consider local de Rham sequences on subdomains and a global finite element space that is broken, i.e. fully discontinuous at the interfaces. The strong differential operators are then applied by ways of conforming projection operators that enforce the proper continuity conditions at the interfaces. This approach yields block-diagonal matrices for L 2 projections (i.e. mass matrices), Hodge operators, and dual commuting projections, while it satisfies the key properties of discrete structure-preservation, such as primal/dual commuting diagrams and discrete Hodge-Helmholtz decompositions [20]. The coupling between subdomains is encoded by the conforming projection operators, which can be highly local: in practice these may only involve the averaging of degrees of freedom across interfaces, yielding low-rank sparse matrices. Ultimately this approach allows us to construct structure-preserving finite element solvers on complex domains, which are easy to implement and efficient.
Specifically, we extend the theory to the discretization of several boundary value problems arising in electromagnetics, and we describe its application to multipatch mapped spline discretizations [16,4] on general non-contractible domains. As our results show, this approach 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 Section 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 Section 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 [24]. For each problem we state a priori results about the solutions, assuming the well-posedness of the corresponding conforming FEEC discretization. In Section 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 Section 5 we conduct extensive numerical experiments for the electromagnetic problems described in the article, which verify the robustness of our approach. Finally, in Section 6 we summarize the main results of our work and provide an outlook on future research.
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 (2.4) could have been chosen for the primal one and the homogeneous (2.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 (2.1)-(2.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 curl grad = 0 and div curl = 0. This allows us to write an orthogonal Hodge-Helmholtz decomposition for L 2 (Ω) 3 , of the form where H 1 := {v ∈ V 1 ∩ V * 1 : curl v = 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 (2.6) L 1 := − grad div + curl curl seen as an operator L 1 : D(L 1 ) → L 2 (Ω) with domain space Another orthogonal decomposition for L 2 (Ω) 3 is where H 2 := {w ∈ V 2 ∩ V * 2 : curl w = div w = 0} is the space of harmonic 2-forms, which coincides with the kernel of the 2-form Hodge-Laplace operator L 2 := − grad div + curl curl with domain space D(L 2 ) := {w ∈ V 2 ∩ V * 2 : curl w ∈ V * 1 and div w ∈ 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 H 1 = H 2 : in the case considered here where the primal sequence has homogeneous boundary conditions, the harmonic 1-forms have vanishing tangential trace on ∂Ω, while the harmonic 2-forms have vanishing normal trace on ∂Ω. In the symmetric case where the non-homogeoneous de Rham sequence is considered as the primal one, one obtains the same decompositions (2.5) and (2.8) but with opposite order: (H ) non-hom. = H 3− . This isomorphism, known as Poincaré duality, is provided by the Hodge star operator; see [2,Sec. 5.6 and 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, H 1 = {0}. However if the domain Ω 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(H 2 ) = dim(H 1 ) non-hom. = b 1 (Ω) is the first Betti number, which counts the number of "tunnels" through the domain, while dim(H 1 ) = dim(H 2 ) non-hom. = b 2 (Ω) is the second Betti number, which counts the number of "voids" enclosed by the domain.

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 ,c h ⊂ V , 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 Π ,c h : V → V ,c h that commute with the differential operators, in the sense that the relation holds for the different operators in the sequence (2.9), In particular, the stability and the accuracy of several discrete problems posed in the sequence (2.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 or L 2 norms, see e.g. [2, Th. 3.9 and 3.19]. Denoting explicitely by the differential operators restricted to the discrete spaces, we define their discrete adjoints (2.14) where ·, · denotes the L 2 (Ω) scalar product. This yields a compatible discretization of both the primal and dual sequences (2.1), (2.3) in strong and weak form, respectively, using the same spaces (2.9). This framework can be summarized in the following diagram projections to the conforming discrete spaces. We observe that these projections commute with the dual differential operators, as a result of (2.14).
Assumption 1 Throughout the article we assume that the primal projection operators Π ,c h are L 2 stable and satisfy the commuting property (2.10). Remark 2.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. [22]. For this purpose V 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 ⊂ V that require more smoothness (or integrability) as the ones in (2.2). We refer to e.g. [43,39,9,22] for some examples, and to Section 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 (2.5) reads The space 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 H 1 , which corresponds to a Betti number of Ω depending on the boundary conditions, see [2, Sec. 5.6 and 6.2].

Broken FEEC discretizations.
In the case where the domain is decomposed in a partition of open subdomains Ω k , k = 1, . . . , 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 (2.20) is that they are naturally equipped with local basis functions Λ i that are supported each on a single patch Ω k with k = k(i). The corresponding mass matrices are then patch-diagonal, i.e., blockdiagonal with blocks corresponding to the different patches, so that their inversionand hence the L 2 projection on V h -can be performed in each patch independently of the others.
In general the spaces (2.20) are 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 [9]: for H 1 smoothness the fields must be continuous on the interfaces, while H(curl) and H(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 (2.20), we have in general The approach developed for the CONGA schemes in [24,20] extends the construction of Section 2.2 to this broken FEEC setting, by associating to each discontinuous space a projection operator on its conforming subspace, Provided that these conforming spaces form a de Rham sequence (2.9), this allows us to define new primal differential operators on the broken spaces We represent this broken-FEEC discretization with the following diagram where Π h andΠ h are projection operators that commute with the new primal and dual sequences. We postpone their description to the next section. We note that the presence of conforming projections P h in the differential operators (2.22) leads to larger kernels. Specifically, we have where the projection kernels correspond intuitively to "jump spaces" associated with the conforming projections. In [19,25] 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 [20] where the jump spaces (2.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 Π ,c h and Π 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 , that are unisolvent in V ,c h and commute with the differential operators in the sense that there exist coefficients D i,j such that where we have denoted again d 0 = grad, d 1 = curl and d 2 = div. Letting Λ ,c i be the basis of V ,c h characterized by the relations σ ,c i (Λ ,c j ) = δ i,j for i, j = 1, . . . N ,c , one verifies indeed that the operators are projections satisfying the commuting property (2.10), see e.g. [22,Lemma 2] and in passing we note that D is the matrix of d in the bases just defined.
Here, we extend this approach as follows: • Broken degrees of freedom. Each local space V h (Ω k ), k = 1, . . . , K, is equipped with unisolvent degrees of freedom On the full domain we simply set • Broken basis functions. To the above degrees of freedom we associate local basis functions Λ k,µ ∈ V h (Ω k ) (extended by zero outside of their patch), characterized by the relations • Local commutation property. A relation similar to (2.26) must hold on each patch Ω k , i.e., there exist coefficients D k,µ,ν such that holds for all µ ∈ M +1 h (Ω k ) and all v ∈ U (Ω k ). • Inter-patch conformity. The global projection on V h , must map smooth functions to conforming finite element fields, namely 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 Section 4 for a detailed construction of broken degrees of freedom that satisfy the above properties.
The dual projection operators can then be defined following the canonical approach of [20], as "filtered" L 2 projections where Q V h is now the L 2 projection on the broken space V h , By definition of the dual differentials (2.23), these operators indeed commute with the dual part of the diagram (2.24), in the sense that hold on V * 1 , V * 2 and V * 3 respectively. We further observe that (2.32)-(2.33) defines projection operators on the subspaces of V h corresponding to the orthogonal complements of the "jump spaces" (2.25) 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 (2.24). To do so we first observe that the non-conforming differential operators d h := d P h can be reformulated as is the patchwise differential operator which coincides with d on all v such that d v ∈ L 2 . Using some implicit flattening (k, µ) → i ∈ {1, . . . , N } for the multi-indices, we represent these operators as matrices of respective sizes N 1 × N 0 , N 2 × N 1 and N 3 × N 2 . These matrices have a "patchdiagonal" structure, in the sense that they are block-diagonal, with blocks corresponding to the differential operators in each independent patch (namely, the matrices D k in (2.29)). The matrices of the conforming projections, seen as endomorphisms in V h , are In general P is not patch-diagonal, as it maps in the coefficient space of the conforming spaces V ,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 nonconforming differential operators d h : they read where we have denoted (D ) = (G, C, D) for brevity. Therefore, A simple representation of the dual discrete operators (2.23) is obtained in the dual bases {Λ i : i = 1, . . . , N } of the broken spaces V h . These are characterized by the relations where the dual degrees of freedom are defined as It follows from (2.38) and (2.39) that the primal and dual bases are in L 2 duality, i.e. Λ i ,Λ j = δ i,j . The matrices of the dual differential operators in the dual bases read then where again (D ) = (G, C, 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 (2.24) as follows: Here C andC are the spaces of scalar coefficients associated with the primal and dual basis functions described in the previous Section. Both are of the form R N with N = dim(V h ), but we use a different notation to emphasize the different roles played by the coefficient vectors. The interpolation operators I : C → V h andĨ :C → V h , which read are the right-inverses of the respective degrees of freedom σ andσ , and also their left-inverses on V h . The matrices H andH are change-of-basis matrices which allow us to go from one sequence to the other, in the sense that I =Ĩ H andĨ = I H : they correspond to discrete Hodge operators [32]. Here it follows from the duality construction (2.38)-(2.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 (2.30) and (2.32)-(2.33), expressed in terms of primal and dual coefficients. Using the interpolation operators (2.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 (2.27)-(2.31), while the commutation of the dual (lower) diagram follows from the weak definition of the dual differential operators.
Remark 2.2 (conforming case). The conforming FEEC diagram (2.15) (corresponding, e.g., to a single patch discretization) may be extended in the same way as the broken-FEEC one (2.24). Apart from the differences between the discrete differential operators and commuting projection operators already visible in (2.15) and (2.24), the resulting 6-rows diagram would be formally the same as (2.41), but no conforming projection matrices would be involved any longer. Remark 2.3 (change of basis). It may happen that the basis functions Λ i associated with the commuting degrees of freedom σ 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 (2.41), such as the ones described in Section 4.8 for local spline spaces.

Discrete
Hodge-Laplace and jump stabilization operators. In addition to the first order differential operators, an interesting feature of the diagram (2.41) is to provide us with natural discretizations for the Hodge-Laplace operators. On the space V 0 this is the standard Laplace operator, L 0 = −∆, which is discretized as with an operator matrix in the primal basis that reads As for the Hodge-Laplace operator for 1-forms, L 1 = curl curl − grad 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 Section 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 L h,α = L h +αS h where is the symmetrized projection operator on the jump space (2.25), and α is a stabilization parameter. Indeed, it was shown in [20] that the kernel of L h,α coincides with that of the conforming discrete Hodge-Laplacian, for any positive value of α > 0.
Finally we note that the corresponding stabilization matrix is readily derived from the operators in the diagram. It reads 3. Application to electromagnetic problems. In this section we propose broken-FEEC approximations for several problems arising in electromagnetics, and we state a priori results for the solutions. For the sake of completeness we 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 Section 2.4. As pointed out in Remark 2.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 Section 4.8 below.
Throughout this section we will mostly work with the homogeneous spaces V corresponding to (2.2). In the few cases where we will need the inhomogeneous spaces, we shall use a specific notationV for the spaces in (2.4), for the sake of clarity.

Poisson's equation.
For the Poisson problem with homogeneous Dirichlet boundary conditions: we combine the stabilized broken-FEEC discretization described in Section 2.7 for the Laplace operator, and a dual commuting projection (2.32)-(2.33) for the source. The resulting problem reads: and it enjoys the following property.
(Ω) and solves the conforming Poisson problem In particular, φ h is independent of both α and the specific projection P 0 h . Proof. By definition of the different operators, (3.2) reads as long as α = 0, and (3.3) follows by considering test functions ϕ in the conforming subspace V 0,c h . The matrix formulation of (3.2) is easily derived by using the diagram (2.41), writing the equation in the dual basis in order to obtain symmetric matrices as is usual with the Poisson problem. Denoting by φ = σ 0 (φ h ) the coefficient vector of φ h in the primal basis, we thus find with a stabilized stiffness matrix

Time-harmonic Maxwell's equation. Maxwell's equation with homoge-
neous boundary conditions reads: given ω ∈ R and J ∈ L 2 (Ω), find u ∈ H 0 (curl; Ω) such that Here the (complex) electric field corresponds to E(t, x) = iω u(x)e −iωt . When ω 2 is not an eigenvalue of the curl curl operator, equation (3.8) is well-posed: see e.g. [3,Th. 8.3.3]. For this problem we propose a stabilized broken-FEEC discretization with a parameter α ∈ R. We remind that S 1 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.
, and solves the conforming FEEC Maxwell problem In particular, u h is independent of both α and the specific projection as long as α = 0, and (3.11) follows by considering test functions v = v c in V 1,c h . Existence and uniqueness follow from the assumption that ω is not an eigenvalue of the conforming curl curl operator.
The matrix form of (3.9) is easily derived from diagram (2.41). It reads where u := σ 1 (u h ) is the (column) vector containing the coefficients of u h in the primal basis of V 1 h , and A 1 is the resulting stabilized stiffness matrix, (3.13)

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. For a Poisson equation of the form this corresponds to introducing φ g ∈ H 1 (Ω) such that φ g = g on ∂Ω, and characterizing the solution as For an inhomogeneous Maxwell equation of the form we introduce u g ∈ H(curl; Ω) such that n × u g = g on ∂Ω, and characterize the solution as u = u g + u 0 , where u 0 ∈ H 0 (curl; Ω) solves (3.8) 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.2) and (2.4), which amounts to combining different conforming projections. For clarity we denote in this section the respective homogeneous and inhomogeneous spaces by V andV . At the discrete level the broken spaces V h have no boundary conditions, so that the distinction only appears in the conforming projections to the spaces , which we naturally denote by P h andP 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 h that approximates g on ∂Ω, and then we compute the homogeneous part of the solution, h , by solving the homogeneous CONGA problem with modified source for all ϕ ∈ V 0 h , which also reads in matrix form where A 0 is the stabilized stiffness matrix (3.7),P 0 is the matrix of the inhomogeneous projectionP 0 h , and the coefficient vectors involve the broken degrees of freedom φ 0 = σ 0 (φ 0,h ), φ g = σ 0 (φ g,h ) in the full broken space V 0 h . For the Maxwell equation where V 1 = H 0 (curl; Ω) andV 1 = H(curl; Ω), the method consists of first computing u g,h ∈ V 1 h such that n×u g,h approximates g on ∂Ω, and then characterizing the homogeneous part of the solution, u 0,h : where A 1 is the stabilized stiffness matrix (3.13),P 1 is the matrix of the inhomogeneous projectionP 1 h , and the coefficient vectors involve the broken degrees of freedom u 0 = σ 1 (u 0,h ), u g = σ 1 (u g,h ) in the full broken space V 1 h .
and it solves the discrete conforming Maxwell equation inside Ω, in the sense that In particular, u c h is independent of both α and P 1 h . Moreover if the lifted boundary condition u g,h is chosen inV 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 Prop. 3.2. This shows that u c h = u 0,h +P 1 h u g,h and that the tangential trace of u c h coincides with that ofP 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 spacē L u = λu are readily derived using the discrete operators presented in Section 2.7. For a general stabilization parameter α ≥ 0, they take the form: Their convergence can be established under the assumption of a strong stabilization regime and of moment-preserving properties for the conforming projections, see [20]. There it has also been shown that for arbitrary positive penalizations α > 0, the zero eigenmodes coincide with their conforming counterparts, namely the conforming harmonic forms, For the first space V 0 = H 1 0 (Ω) the corresponding operator is the standard Laplace operator L 0 = −∆ and its broken-FEEC discretization reads: where L 0 h and S 0 h are given by (2.45) and (2.49), with stabilization parameter α. For V 1 = H 0 (curl; Ω) the Hodge-Laplace operator is L 1 = curl curl − grad div and its discretization reads: (2.47), and S 1 h is given again by (2.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 (3.25) used for the Hodge-Laplace operator, by dropping the term (− grad h div h ). By testing such an equality with v ∈ grad P 0 h V 0 h we see that the non-zero eigenmodes satisfy div h u h = 0, hence they are also eigenmodes of (3.25) and their convergence follows from the results mentionned above in the case of strong penalization regimes. In the unpenalized case (α = 0) their convergence was also established, as shown in [24].
The main drawback of the approach just described is that the computed eigenvalues do not coincide with those of the conforming operator curl c h curl c h , but only approximate them. (We recall that the knowledge of the eigenvalues of curl c h curl c h is needed in order to verify the hypotheses of Proposition 3.2 for the source problem: ω 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 Proposition 3.4. For λ h = 0, the CONGA eigenvalue problem (3.26) has the same solutions as the conforming eigenproblem h . The relationship (3.29) between the null spaces follows from the equalities ker curl h curl h = ker(curl h ), which was observed in, e.g., [24].
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 Section 3.1 and 3.2. For the Poisson eigenvalue problem (3.24) we obtain . Similarly we can write the Hodge-Laplace eigenvalue problem (3.25) as and the curl-curl eigenvalue problem (3.26) as 3.5. Magnetostatic problems with harmonic constraints. We next consider a magnetostatic problem of the form posed for a current J ∈ L 2 (Ω). 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 curl ∩ ker 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 × B = 0, which is sometimes used to model pseudo-vacuum boundaries [37]. A mixed formulation is then: Find B ∈ H 0 (curl; Ω) and p ∈ H 1 0 (Ω) such that see e.g. [5]. Here the auxiliary unknown p is a Lagrange multiplier for the divergence constraint written in weak form.
The well-posedness of (3.34) is related to the space of harmonic 1-forms H 1 , defined as the kernel of the Hodge-Laplace operator (2.6) in the homogeneous sequence, i.e., (3.35) As discussed in [2], H 1 = H 1 (Ω) is isomorphic to a de Rham cohomology group and its dimension corresponds to a Betti number of the domain Ω. With homogeneous, resp. inhomogeneous boundary conditions its dimension is b 2 (Ω) the number of cavities, resp. b 1 (Ω) the number of handles in Ω. Thus in a contractible domain we have H 1 = {0} and the problem is well posed, but in general H 1 may not be trivial. Additional constraints must then be imposed on the solution, such as B ∈ (H 1 ) ⊥ , which can be associated with an additional Lagrange multiplier z ∈ H 1 . The constrained problem then consists of finding B ∈ H 0 (curl; Ω), p ∈ H 1 0 (Ω) and z ∈ H 1 , such that In our broken-FEEC framework this problem is discretized as: is well-posed for arbitrary non-zero stabilization parameters α 0 , α 1 = 0, and its solution satisfies Moreover the discrete magnetic field can be decomposed according to (2.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 (2.17). In par- 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 ∈ V 1,c h , an arbitrary conforming test function. Remark 3.6. Prop. 3.5 states that the auxiliary unknowns p h and z h may be discarded a posteriori from the system. However they are needed in (3.37) to have a square matrix.
In practice the discrete harmonic fields can be represented by a basis of H 1 h , of the form , computed as the zero eigenmodes of a penalized is the Betti number of order 2, i.e. the number of cavities in Ω. Using these harmonic fields, one assembles the rectangular mass matrix which allows us to rewrite (3.37) as where p, B and z contain the coefficients of p h , B h and z h in the (primal) bases of V 0 h , V 1 h and H 1 h respectively. 3.5.2. Magnetostatic problem with metallic boundary conditions. We next consider a boundary condition of the form n · B = 0 which corresponds to a metallic (perfectly conducting) boundary. In our framework it can be approximated by using the inhomogeneous sequence (2.4) which we denote byV as in Section 3.3, and by modelling the boundary condition in a weak form. The mixed formulation takes a form similar to (3.34) with a few changes. Namely, we now look for p ∈V 0 = H 1 (Ω) and B ∈V 1 = H(curl; Ω), solution to (3.44) grad q, B = 0 Compared to (3.34), 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| ∂Ω = 0 which enforces div B = 0, and a second set with q| ∂Ω = 0 which enforces n · 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 (2.4) involvingV 1 = H(curl; Ω) andV * 1 = H 0 (div; Ω). 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 ∈ (ker grad) ⊥ . In practice this constraint may either be enforced in a similar way as the harmonic constraint, or by a socalled regularization technique where a term of the form ε q, p , with ε = 0, is added to the first equation of (3.44). The constrained problem reads then: Find B ∈ H(curl; Ω), p ∈ H 1 (Ω) and z ∈H 1 , such that We note that here the result does not depend on ε (indeed by testing with v = grad p we find that p is a constant, and with q = p we find that p = 0), hence we may set ε = 1.
At the discrete level we write a system similar to (3.37) whereP 0 andP 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 Section 3.5.1. The resulting system reads: h . An operator form similar to (3.39) can also be written.
Finally the matrix form of our broken-FEEC magnetostatic equation with metallic boundary conditions is similar to (3.43), with the following differences: • the conforming projection operators map in the full conforming spacesV (boundary dofs are not set to zero), • the discrete harmonic space is nowH 1 h the kernel of theL 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 M 0 p involving the mass matrix in V 0 h is added to the first equation to fully determine p.
Proposition 3.7. System (3.50) is well-posed for arbitrary stabilization parameters α 0 > 0, α 1 = 0, and its solution satisfies Moreover the discrete magnetic field can be decomposed according to (3.48), as with discrete grad, harmonic and curl components characterized by In particular, the solution is independent of the parameters α 0 > 0, α 1 = 0, and the conforming projection operators.
Proof. Taking v = z h (inH 1 h ⊂V 1,c h ) and using the orthogonality (3.48) shows that z h = 0. Taking next v = gradP 0 h p h and using again (3.48) shows thatP 0 h p h is a constant, and with q = 1 we see that 1, p h = 0. It follows that with q = (I −P 0 h )p h the first equation becomes p h Since it is orthogonal to bothH 1 h and gradV 0,c h , it belongs to curl c h V 2,c h , and it is characterized by curl v, curl B h = curl v, J for all v ∈V 1,c h , which completes the proof. 3.6. Time-dependent Maxwell's equations. We conclude this section by recalling the CONGA discretization of Maxwell's time-dependent equations, with a discrete source approximated with the dual commuting projection This discretization has been proposed and studied in [24,25,26], 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 div hΠ 1 h J =Π 0 h div J satisfied by the dual projection operators, see (2.34). Moreover, the filtering by (P 1 h ) * involved in the source approximation operator (2.32) avoids the accumulation of approximation errors in the non-conforming part of the kernel of the CONGA curl curl operator [24].
Remark 3.8. The commuting diagram property (2.34) is actually not needed for the preservation of the discrete Gauss law. One can indeed verify that (3.56) still holds for a current source defined by a simple L 2 projection, so that (3.56) follows from the compatibility relation div J + ∂ t ρ = 0 satisfied by the exact sources. The additional filtering of the source by (P 1 h ) * involved in the dual commuting projection (3.55) is nevertheless important for long term stability, as analyzed in [24] and verified numerically in Section 5.7 below.
An attractive feature of the CONGA discretization is that for an explicit timestepping method such as the standard leap-frog scheme h each time step is purely local in space. This is easily seen in the matrix form of (3.58) where all the matrices except P 1 are patch-diagonal, and 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 (3.56) 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 valid for all n ≥ 0. In practice we set ∆t = C cfl (2/ | curl h |) with C cfl ≈ 0.8, and we evaluate the operator norm of curl h with an iterative power method, noting that it amounts to computing the spectral radius of a diagonalizable matrix with nonnegative eigenvalues, namely

Geometric broken-FEEC discretizations with mapped patches.
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, For simplicity we consider the case of planar domains Ω k ⊂ R 2 and orientation-preserving mappings, in the sense that the Jacobian matrices DF k (x) have positive Jacobian determinants J F k (x) = det(DF k (x)) > 0 for allx ∈Ω. We also assume that the multipatch decomposition is geometrically conforming, in a sense that is specified in Section 4.6 Examples of such domains are represented on Figure 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. [28], 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.

De Rham sequences in 2D.
In 2D the de Rham sequence (2.1) may be reduced in two different ways, namely where curl and curl denote the vector and scalar-valued curl operators respectively (for clarity we use bold fonts to denote vector-valued fields and operators in the subsequent sections). Here the operators in the second sequence are dual to those in the first one, and in each case we may consider spaces with homogeneous boundary conditions, namely In the numerical examples presented in this article we will consider the first sequence (4.2) as the primal one, both with homogeneous spaces (4.4) or inhomogeneous ones (4.6). As mentionned in Section 3.3 this distinction will only affect the conforming projection operators, indeed the broken spaces consist of local spaces which have no boundary conditions. This point will be further discussed in Remark 4.4 below.

Reference spline complexes.
On each mapped patch the local sequence (2.19) will be built as the push-forward of tensor-product spline spaces defined on the reference patchΩ = ]0, 1[ 2 following [17,16], and we shall equip it with geometric degrees of freedom. In order to simplify the construction of the conforming spaces {V ,c h } over domains with non-trivial inter-patch connectivities, we will define a discretization which has reflectional and rotational symmetries overΩ. This amounts to using the same knot sequence along every dimension ofΩ, and to require such a knot sequence to be left-right symmetric as explained below.
We now recall the construction of the reference tensor-product spline complex with maximal coordinate degree p ≥ 1 on a grid with N cells per dimension. For this we equip the interval [0, 1] with an open knot sequence of the form (4.8) 0 = ξ 0 = · · · = ξ p < ξ p+1 < · · · < ξ n−1 < ξ n = · · · = ξ n+p = 1 where n = N + p, and assumed symmetric for simplicity, in the sense that ξ i = 1 − ξ n+p−i for all i. For i = 0, . . . , n − 1 and q ∈ {p − 1, p} we then let N q i be the normalized B-spline of degree q associated with the knots (ξ i , . . . , ξ i+q+1 ), see The reference spline complex [17] associated with the grad − curl sequence (4.2) reads then

4.3.
Geometric degrees of freedom on the reference patch. Following [6,45,30,22] we next equip the reference complexes (4.10) and (4.11) with geometric degrees of freedom which are known to commute with the differential operators. They are based on an interpolation grid 0 = ζ 0 < · · · < ζ n−1 = 1 for the univariate spline space S p in the sense of [46,Th. 4.61], that is also symmetric, namely ζ i = 1 − ζ n−1−i for all i = 0, . . . , n − 1. Here the usual choice is to consider Greville points associated with the knots (4.8). On the reference domainΩ we then consider the interpolatory nodes, edges and cells Here the square brackets [·] denote convex hulls, e d is the canonical basis vector of R 2 along dimension d and the multi-index sets read Geometric degrees of freedom associated with the grad − curl sequence (4.2) are then defined as the linear forms (4.13) where we observe that the reference edges are oriented according to their natural parametrization, where e ⊥ d := (e 2 , −e 1 ) T corresponds to a π/2 rotation (counterclockwise) of the basis vector e d . As these degrees of freedom are unisolvent we letΛ µ , µ ∈M , be the geometric basis for the spline spaceV , characterized by the duality relations Hence, we need to specify some proper domain spaces.
for the grad − curl sequence (4.2), and for the curl − div sequence (4.3), with anisotropic Sobolev spaces defined as Moreover they commute with the local differential operators, namely we have with graph-incidence Kronecker-product matricesD .
Proof. Nodal degrees of freedomσ 0 µ are well defined onÛ 0 L 1 = W 1 1,2 (Ω) because this space is continuously imbedded in C 0 (Ω): this can be verified using Remarks 9 and 13 from [14, Sec. 9], and a density argument. Next we observe that edge degrees of freedom (4.13) along horizontal edges e 1,i are of the form Similarly we see thatσ 2,i is also well defined onÛ 1 L 1 , and the other degrees of freedom are analyzed with the same reasoning. We also verify easily that in both cases, namely (4.18) and (4.19), the spacesÛ L 1 form a de Rham sequence. Taking the intersection with (4.17) thus yields a subsequenceÛ ⊂V where the degrees of freedom are indeed well-defined. As for the commuting property, it is well known and follows from the Stokes formula, see e.g. [6]. We refer to [22] for a description of the graph-incidence and Kronecker-product structure of the matricesD .

4.4.
Broken FEEC spaces on the mapped patches. Following [33,16,44] we define the local spaces (2.19) on the mapped patches Ω k = F k (Ω) as push-forwards of the reference spline spaces, namely where the push-forward transforms associated with a mapping F k are obtained as the 2D reduction of the usual ones in 3D. For the grad − curl sequence (4.2) they read and for the curl − div sequence (4.3) they read denotes the Jacobian matrix of F k , and J F k its (positive) metric determinant corresponding to the surface measure for twodimensional manifolds in R 3 , see e.g. [33,39,15,35]. The inverse transforms (F k ) −1 are called pull-backs, and a fundamental property of these transforms is that they commute with the differential operators [33,30], namely holds with (d 0 , d 1 ) the differential operators in the sequences (4.2) or (4.3), and (d 0 ,d 1 ) their counterparts on the logical (reference) domain. As these transforms are linear operators, the local spaces are spanned by the push-forwarded basis functions. In particular, setting

Geometric degrees of freedom on the mapped patches.
In connection with the multipatch structure of V h we define broken degrees of freedom using pull-back transforms, which are the inverse of the push-forward operators (4.22) or (4.23), to obtain By construction these degrees of freedom are in duality with the local basis functions (4.25), i.e.
A key feature of the pull-back operators is to carry the geometric nature of the reference degrees of freedom to the mapped elements. For = 0 the degrees of freedom simply consist of pointwise evaluations on the mapped nodes, i.e.
(4.29) σ 0 k,i (v) = v| Ω k (n k,i ) with n k,i := F k (n i ). For = 1 with the grad − curl sequence (4.2), they correspond to line integrals along the mapped edges e k,µ := F k (ê µ ) with µ = (d, i) ∈M 1 . Specifically, using the pullbackv k : (4.22) and the parametrization (4.14), we have where τ is the unit vector tangent to e k,µ with the orientation inherited from that of e µ , as mapped by F k : at x = x e k,µ (s) = F k (x e µ (s)) it reads For = 1 in the curl − div sequence (4.3), the pull-back corresponding to (4.23) readŝ where τ ⊥ is the result of a π/2 rotation (counterclockwise) of the unit tangent vector (4.31), and reads Finally for = 2 the degrees of freedom are the integrals on the mapped cells, v| Ω k with c k,i := F k (ĉ i ).  Proof. The first statement is a direct consequence of the definition (4.27) and Prop. 4.1. In particular the fact that the spaces (4.35) form a sequence follows from the similar property of the reference spaces and the commutation (4.24) of the pullbacks. As for the local commuting property, it also follows from (4.24) and the similar property (4.20) of the reference degrees of freedom, indeed 4.6. Conformity of the multipatch geometry. We next assume that the patches are fully conforming in the sense that any interface Γ k,l = ∂Ω k ∩ ∂Ω l between any two distinct patches (i) is either a vertex, or a full edge of both patches, (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 wherex 0 andŷ 0 are vertices of the reference domainΩ, condition (ii) means that holds with θ k,l an affine bijection on [0, 1], that is θ k,l (s) = s or 1 − s.
Since the reference knot sequence was supposed symmetric in Section 4.2, these conditions imply that the patches are fully matching in the sense of Assumption 3.3 from [18] (applied to the scalar spline spaces in the sequence).

Conforming projection operators.
An important property of the geometric degrees of freedom is that they provide us with a simple characterization of the discrete fields v h in the broken spaces V h which actually belong to the conforming spaces V ,c h := V h ∩ V . In order to belong to the space V appearing in the continuous de Rham sequence (2.1), a field must satisfy regularity conditions which are well known for piecewise-smooth functions. For = 0, the condition for H 1 (Ω) 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 = 1 in the grad − curl sequence (4.2), the interface constraint for H(curl; Ω) regularity consists in the continuity of the tangential traces, which amounts to requiring that edge degrees of freedom (4.30) associated to the same edge coincide, up to the orientation: Here ε e (k, µ; l, ν) = ±1 is the relative orientation of the mapped edges e k,µ and e l,ν , characterized by the relation where x e k,µ (s) = ∂F k (x e µ (s)) and similarly for x e l,ν , see (4.31), and θ k,l is the affine bijection involved in the interpatch conformity assumption (4.36)-(4.37). For = 1 in the curl − div sequence (4.3), the interface constraint for H(div; Ω) 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 (4.39) with the corresponding definition (4.32). Finally as V 2 = L 2 (Ω) there are no constraints for = 2, i.e., the spaces V 2 h and V 2,c h coincide. We gather these constraints in a single formula where g k,µ denotes the geometric element (node, edge or cell) of dimension associated with the multi-index (k, µ) ∈ M h , ε 1 denotes the relative orientation of two edges, and ε 0 = ε 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 ,c h . In particular we may define simple conforming projection operators P by averaging the broken degrees of freedom associated with interface elements. Using the broken basis functions (4.25) this yields an expression similar to the one given in [20] for tensor-product polynomial elements, with additional relative orientation factors due to the general mapping configurations, where M h (g l,ν ) contains the broken indices corresponding to the same geometric element. The entries of the corresponding operator matrix (2.36) read then 0 otherwise for all (k, µ), (l, ν) ∈ M h . We may summarize our construction with the following result.
Proposition 4.3. The broken degrees of freedom defined in Section 4.5 satisfy the properties listed in Section 2.4, with local domain spaces U (Ω k ) defined in (4.35) and local differential matrices independent of the mappings F k . The operators P h : V h → V h defined by (4.41) are projections on the conforming subspaces V ,c h = V h ∩ V . Proof. The boundedness of σ k,µ on the local spaces U (Ω k ), as well as the commuting properties, have been verified in Prop. 4.2. The inter-patch conformity property (2.31) follows from the fact that piecewise smooth functions that belong to H 1 (Ω), resp. H(curl; Ω) and H(div; Ω), admit a unique trace, resp. tangential and normal trace on patch interfaces [9]. The same property, together with the interpolation nature of the geometric basis functions, allows verifying that the operators P h : V h → V h are characterized by the relations for all (k, µ) ∈ M h . Using the geometric condition (4.40), this allows verifying that P h is indeed a projection on the conforming subspace V ,c h . Remark 4.4 (boundary conditions). As stated, the above construction actually corresponds to the inhomogeneous sequences (4.6) and (4.7). If the conforming spaces are homogeneous, as in (4.4) or (4.5), then P h should further set the boundary degree of freedom to 0, which amounts to restricting the non-zeros entries in (4.42) to the geometrical elements g k,µ that are inside Ω. We note that this does not affect the broken spaces V h , as already observed in Section 3.3, since they are not required to have boundary conditions. Remark 4.5 (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 [16,4], and the same geometric degrees of freedom as in [22,Sec. 6.1] for the primal sequence.
4.8. Primal-dual matrix diagram with B-splines. In practice, a natural approach is to work with B-splines. Indeed the geometric (interpolatory) splines Λ k,µ introduced in Sections 4.3 and 4.4 are not known explicitely: they depend on the degrees of freedom, i.e. on the interpolation grids (4.12). They are also less local than the B-splines, as they are in general supported on their full patch Ω 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 [4,22] and used tensor-product splinesB µ on the reference domainΩ, composed of normalized B-splines N p i in the dimensions of degree p and of Curry-Schoenberg B-splines D p−1 (also called M-splines) in the dimensions of degree (p − 1). On the mapped patches Ω k = F k (Ω) the basis functions are defined again as push-forwards B k,µ : The discrete elements in the matrix diagram (2.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 the second equality uses (4.27). From B k,ν = µ K (k,µ),(k,ν) Λ k,µ it follows that the B-spline coefficients of the geometric projection (2.30), namely read (using again some implicit flattening M h (k, µ) → i ∈ {1, . . . , N }) Accordingly, we obtain the matrices P B = (β i (P h B j )) 1≤i,j≤N of the conforming projection operator (4.41) in the B-spline bases by combining the matrices (4.42) with the above change of basis: this gives (4.47) Here we note that the collocation matrices K are Kronecker products of univariate banded matrices, see [22,Sec. 6.2], so that (4.46) and (4.47) 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
, these are patch-wise incidence matrices, just as the ones in the geometric bases [4,22]. For the dual sequence we use bi-orthogonal splines characterized by the relations B k,µ ∈ V h (Ω k ),β k,ν (B k,µ ) = δ µ,ν ∀µ, ν ∈M with dual degrees of freedom This leads to defining again the primal Hodge matrices as patch-diagonal mass matrices, H B = M B and the dual Hodge ones as their patch-diagonal inversesH B = (M B ) −1 . These matrices can be computed on the reference domain according to Specifically, for the 2D grad − curl sequence (4.2), the transforms (4.22) yield and for the 2D curl − div sequence (4.3) we find from (4.23) that while M 0 B and M 2 B take the same values as in (4.50). In 3D the formulas can be derived from the pull-backs given in [4].
Finally the coefficients of the dual commuting projections (2.32)-(2.33) in the dual B-spline bases readβ using the dual degrees of freedom (4.48) and the same implicit flattening as in (4.46), (4.47).
Gathering the above elements we obtain a new version of diagram (2.41), where the coefficient spaces (still defined as R N ) are now denoted C B andC B to indicate that they contain coefficients in the B-spline and dual B-spline basis, respectively. (4.52) Note that here the primal and dual interpolation operators are simply (4.53)

Numerical results.
In this section we illustrate several of the numerical schemes described in Section 3, using the 2D multipatch spline spaces presented above. The results have been obtained with the Psydac library [47].

5.1.
Poisson problem with homogeneous boundary conditions. We first test our CONGA scheme for the homogeneous Poisson problem presented in Section 3.1. For this we consider an analytical solution on the pretzel-shaped domain shown in Figure 1, given by Here, s(x) =x−ỹ, t(x) =x+ỹ withx = x−x 0 ,ỹ = y −y 0 and we take x 0 = y 0 = 1.5, a = (1/1.7) 2 , b = (1/1.1) 2 and σ = 0.11 in order to satisfy the homogeneous boundary condition with accuracy ≈ 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 Figure 2 we plot the numerical solutions corresponding to spline elements of degree 3 × 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.     Table 1.

Poisson problem with inhomogeneous boundary conditions.
Our second test is with an inhomogeneous Poisson-Dirichlet problem (3.14), using the lifting method described in Section 3.3 for the boundary condition and the solver tested just above for the homogeneous part of the solution. Specifically, we define φ g,h ∈ V 0 h by computing its boundary degrees of freedom from the data g on ∂Ω (this is straightforward with our geometric boundary degrees of freedom (4.29)), and we compute φ 0,h = φ h − φ g,h by solving (3.17). 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 φ(x) = sin(πx) cos(πy) is the exact solution.
In Figure 3 we plot the convergence curves corresponding to various N × N grids for the 18 patches of the domain, and various degrees p = 2, . . . 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, which confirms the numerical accuracy of our stabilized CONGA formulation for the Poisson problem with smooth solutions. where φ and τ are as in (5.1). We plot this source in Figure 4. We then consider two values for the time pulsation, namely ω = √ 50 and ω = √ 170, and for each of these values we use as reference solutions the numerical solutions computed using our method on a mesh with 20 × 20 cells per patch (i.e. 7200 cells in total) and elements of degree 6 × 6 in each patch. These reference solutions are shown in Figure 5. Interestingly, we observe that for the higher pulsation ω = √ 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 ω, and we also plot in Figure 6 the different solutions corresponding to spline elements of degree 3 × 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 ω = √ 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            Table 2.
the lifting method described in Section 3.3. Specifically, we define u g,h ∈ V 1 h by computing its boundary degrees of freedom from the data g = n × u on ∂Ω (this is again straightforward with our geometric boundary degrees of freedom (4.30)), and we compute u 0,h = u h − u g,h by solving (3.19). Here we take ω = π and consider (3.15) with the source-solution pair (5.5) J = −π 2 sin(πy) cos(πx) 0 , u = sin(πy) sin(πx) cos(πy) and a boundary condition given by g := n × u on ∂Ω.
In Figure 7 we plot the convergence curves corresponding to spline elements of various degrees p × p and N × 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. 5.5. Eigenvalue problems. We next assess the accuracy of our CONGA approximation (3.26) for the curl-curl eigenvalue problem.
We test our discretization on the two domains shown in Figure 1, and plot in Figure 8 the amplitude of the first five eigenmodes, together with their positive eigenvalues. Here the eigenmodes are computed using spline elements of degree 6 × 6 and 56 × 56, resp. 20 × 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 22692 degrees of freedom for the broken space V 1 h , while on the latter domain it corresponds to 7200 cells and 23400 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 Figure 9 we then plot the relative eigenvalue errors as a function of the eigenvalue index i, for degrees p = 3 and 5, and N × 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 [28,29], 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 × 6 and N × N cells per patch with N ≤ 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 Figure 8. 5.6. Magnetostatic test-cases. We next study the CONGA discretizations of the magnetostatic problems presented in Section 3.5: either the one for the problem with pseudo-vacuum boundary conditions (3.36) or the one with metallic boundary conditions (3.46).
To this end we consider a scalar dipole current source, in the pretzel-shaped domain, with σ = 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 Section 3.5.1 and 3.5.2 for the problems with pseudo-vacuum and metallic boundary conditions, respectively. In Figure 10 we plot the scalar source J z together with the vector-valued curl J z field, and for each of the boundary conditions, we plot in Figure 11 fine solutions computed using the CONGA scheme on a mesh with 20 × 20 cells per patch (i.e. 7200 cells in total) and elements of degree 6 × 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 Figure 11. We also plot in Figure 12 the solutions corresponding to spline elements of degree 3 × 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.

Time-dependent Maxwell equation.
We finally assess the qualitative and quantitative properties of our mapped spline-based CONGA scheme for the timedependent Maxwell system, using a leap-frog time stepping (3.58)-(3.59). 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 σ = 0.02. The initial magnetic field and the source is (5.9) B(t = 0) = 0, and J (t, x) = 0.
In Figure 13 we compare successive snapshots of two solutions corresponding to different numbers N × N of cells per patch and degrees p × 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 Section 3.6. In Figure 14 we plot the time evolution    , 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 (3.56), indeed we have ρ = div E(t = 0) = 0 in this test case. On Figure 14 we verify that it is zero up to machine accuracy and small quadrature errors.
We next study the quality of the source approximation operator J → J h which is involved in the discrete Ampere equation (3.54). For this we use a second timedependent test-case with a zero initial condition         In Figures 15-16 we compare three different approximation operators for the current source, namely: (i) the primal finite element projection:  (ii) the L 2 projection on the broken space: J h = P V 1 h J , (iii) the dual projection: J h =Π 1 h J . We note that each of these projection operators are local, in the sense that none requires solving a global problem on the computational domain Ω. We also remind that the primal projection Π 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 (3.56). In contrast, both the L 2 projection on the broken space V 1 h and the dual projectionΠ 1 h guarantee the preservation of the discrete Gauss laws, however we expect an increased stability for the latter one as discussed in Remark 3.8.
In Figure 15 we first show some snapshots of the three numerical solutions on a time range t ∈ [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 Figure 16 two error indicators for each one of the numerical strategies described above, namely the Gauss law errors associated with the broken solution G n h (E n h ), and that of its conforming projection G n h (P 1 h 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 3.8: 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 Π 1 h , with a linear time growth (which is somehow consistent with the strong growth of the former error indicator G n h (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Π 1 h , and it oscillates but does not seem to grow. These results tend to indicate that the growing field visible in Figure 15 corresponds to a numerical error, and they highlight the improved stability of the CONGA scheme with a proper source approximation.
6. Conclusions. 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 [20], 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 twodimensional 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.
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.  h E n h | are shown at t = 2, 5, 10 and 20 (from top to bottom). In the left panels the source projection J → J h is a primal projection operator Π 1 h which does not commute with the dual differential operators (note the time-varying color scale). The middle panels use an orthogonal projection P V 1 h , and the ones on the right a dual commuting projectionΠ 1 h . See the text for more details.