Trefftz co-chain calculus

We are concerned with a special class of discretizations of general linear transmission problems stated in the calculus of differential forms and posed on Rn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^n$$\end{document}. In the spirit of domain decomposition, we partition Rn=Ω∪Γ∪Ω+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^n=\Omega \cup \Gamma \cup \Omega _{+}$$\end{document}, Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} a bounded Lipschitz polyhedron, Γ:=∂Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma :=\partial \Omega $$\end{document}, and Ω+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{+}$$\end{document} unbounded. In Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document}, we employ a mesh-based discrete co-chain model for differential forms, which includes schemes like finite element exterior calculus and discrete exterior calculus. In Ω+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{+}$$\end{document}, we rely on a meshless Trefftz–Galerkin approach, i.e., we use special solutions of the homogeneous PDE as trial and test functions. Our key contribution is a unified way to couple the different discretizations across Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}. Based on the theory of discrete Hodge operators, we derive the resulting linear system of equations. As a concrete application, we discuss an eddy-current problem in frequency domain, for which we also give numerical results.


Linear partial differential equations in exterior calculus
The calculus of differential forms, also known as exterior calculus, is a powerful tool for expressing a host of apparently different linear boundary value problems in a unifying way through differential forms; see [23, p. 265, Sections 1 and 2]. Concretely, we write Λ k (R n ) for the space of differential forms of order k, 0 ≤ k ≤ n, n ∈ N * , in R n [8,p. 13,Section 2.2].
Each member of the family of partial differential equations we are concerned with in this work can naturally be split into two sets of equations. The first involves the exterior derivative operator d and comprises the equilibrium equations du = (−1) l σ, connecting the differential forms u ∈ Λ l−1 (R n ), σ ∈ Λ l (R n ), j ∈ Λ m (R n ), ψ ∈ Λ m+1 (R n ), and the known source ψ 0 ∈ Λ m+1 (R n ) [25, p. 243, Table 2.1], given l ∈ {1, . . . , n} and m := n − l.
The other set is formed by the constitutive equations The symbols α and γ stand for Hodge operators, which supply linear mappings of l-forms into m-forms [8, p. 12]. These are induced by the material metrics α and γ, which are Riemannian metrics in the simplest case. However, α and γ can be more general (pseudo-Riemannian) metrics. In particular, γ can also be • zero, which is the static case, • negative, which is the case of a Helmholtz-type equation describing wave propagation, or • imaginary, which leads to a dissipation problem, describing, e.g., eddy currents: this is the case of (1.2b).
If R n is equipped with Cartesian coordinates, these metrics can be represented by invertible matrix fields.

Remark 1.
Note that the equations are posed on the unbounded domain R n and, thus, we have to impose suitable decay or radiation conditions for x → ∞ on the forms, or, equivalently, select the forms from weighted Sobolev spaces; see [39,p. 42,Section 1.4], [43, p. 1484, Section 4], [37, p. 602, Section 3]. The concrete conditions will depend on the signs of α and γ.

Concrete case
Let us connect the general problem of Sect. 1.1 to a concrete boundary value problem relevant in computational electromagnetics, namely the eddy-current problem in frequency domain in three dimensions [1]. For (1.1), this amounts to the case n = 3, l = 2, and m = 1. The electromagnetic fields involveddesignated with the customary symbols-are the magnetic vector potential A : R 3 → C 3 , the magnetic flux density B : R 3 → C 3 , the magnetic field H : R 3 → C 3 , the current density J : R 3 → C 3 , and the known source current density J 0 : R 3 → C 3 . These fields are the so-called Euclidean vector proxies [25, p. 243, Table 2.1] of the differential forms u ∈ Λ 1 (R 3 ), σ ∈ Λ 2 (R 3 ), j ∈ Λ 1 (R 3 ), and ψ, ψ 0 ∈ Λ 2 (R 3 ) of Sect. 1

Overview and outline
In this article, we present a domain-decomposition approach for the discretization of (1.1). We introduce a bounded domain Ω ⊂ R n , Γ := ∂Ω, in whose complement Ω + := R n \ Ω we assume constant α, γ ∈ C, that is α and γ are induced by constant multiples of the Euclidean metric. Thus, we can view (1.1) as a transmission problem coupled across the interface Γ. Our main idea is to rely on different discretizations in Ω and Ω + : • In Ω, we want to employ a generic mesh-based discretization, which falls in the class of co-chain-based discretizations. This framework is the generalization of Finite Element Exterior Calculus (FEEC) [8] and Discrete Exterior Calculus (DEC) [22,27,38] and we elaborate it in Sect. 2. • In the complement Ω + of Ω, filled with a homogeneous medium, we use a Trefftz-Galerkin approximation. It relies on basis functions that solve the homogeneous PDE, satisfy the right decay conditions, and span trial and test spaces to be plugged into variational formulations. Details will be given in Sect. 3.
The key issue is how to connect the different discretizations. Our proposals for this coupling represent the main novel contributions of this work. To develop our ideas, in Sect. 1.4 we examine different (variational) formulations of the transmission problem arising from (1.1). Then, Sect. 4 discusses how to couple the discretizations across the interface Γ separating Ω and Ω + . Starting from the natural transmission conditions for differential forms, we derive the fully-discrete coupled system of equations. This generalizes our earlier approaches to couple the finite element method with Trefftz methods: see [15] for Poisson's equation, [16] for 2D Helmholtz equation, [17] for magnetostatic Maxwell's equations, [19] for an eddy-current model as in Sect. 1.2, and [18] for electromagnetic wave propagation.
Finally, we return to the concrete eddy-current problem of Sect. 1.2 in Sects. 5 and 6. In Sect. 5, we show some examples of Trefftz basis functions and corresponding pure Trefftz blocks of the coupled systems. After that (Sect. 6), we also report simulation results for the concrete problem of Sect. 1.2 using the cell method (Sect. 6.1) for the mesh-based discretization.
The focus of this article is on the derivation and presentation of coupled discrete models for (1.1). It does not include a numerical analysis of convergence because rigorous theoretical results are even missing for some variants of the standalone cell method. Theoretical results are available in special cases of the finite element method and we refer the reader to the publications listed above.

Variational formulations
We adopt the domain-decomposition setting introduced in Sect. 1.3: R n = Ω ∪ Γ ∪ Ω + . For the sake of simplicity, we assume that the known source term ψ 0 ∈ Λ m+1 (R n ) is supported inside Ω. We elaborate two ways to arrive at a general transmission problem in R n and establish related variational formulations.
(1.7) forms the basis of the coupling approach for this problem (Sect. 4.1): 1. on the bounded domain Ω, (1.7) is posed on a suitable Sobolev space of (l − 1)-forms with squareintegrable exterior derivative [8,25], while 2. the influence of Ω + is taken into account through the interface term of (1.7).

Option II: hybrid formulation.
With the term "hybrid," we refer to a mixed formulation involving a further potential for the exterior problem. We assume • γ = 0 in Ω + , which implies ψ| Ω+ = 0 and dj| Ω+ = 0 from ( We can then introduce a potential π ∈ Λ m−1 (Ω + ). The equations (1.1b) for the exterior problem become j = dπ, To complete the problem, a suitable decay condition has to be imposed.

Remark 3.
Throughout we rely on ungauged formulations even if γ = 0 locally. The uniqueness of solution components σ and j will not be affected.
The equation (1.1) for γ = 0 allows to write α du = (−1) l j when restricted to Γ. Thus, we can recast the primary variational problem for u ∈ Λ l−1 (Ω) and j ∈ Λ m (Ω) as with j from the solution of (1.8). Concrete Case. For the concrete case of the eddy-current problem presented in Sect. 1.2, the hybrid approach involves a magnetic scalar potential in Ω + . We assume that σ = 0 in Ω + , which implies J| Ω+ = 0 and ∇ × H| Ω+ = 0 from (1.2a) and (1.2b). We then exploit that ∇ × H = 0 in Ω + to switch to a potential representation there; Ω + is supposed to have trivial topology, i.e., first Betti number 2 β 1 (Ω + ) = 0 [25, p. 252, Lemma 2.2]. This can be achieved by choosing Ω + simply connected, if Ω has no handles. We can then introduce a magnetic scalar potential φ : R 3 → C in Ω + . Hence, the eddy-current model in Ω + is converted into with ν + = ν in Ω + and the decay condition (1.11)

Co-chain calculus
The framework of co-chain calculus, explained in [10,23,24,42] and [12,p. 155,Section 16], allows for a unified treatment of a wide class of finite element and finite volume schemes, building on the foundation established by other works like [11]. This framework is the generalization of Finite Element Exterior Calculus (FEEC) [7,8], Discrete Exterior Calculus (DEC) [22,27,38], and the cell method [20,40]. The starting point of co-chain calculus is the linear stationary or time-harmonic elliptic boundary value problem in differential forms (1.1) on Ω. We aim at discretizing those forms 3 by means of a cochain representation. To that end, we introduce 1. a primary mesh M, with F k set of k-dimensional facets of such mesh and N k := #F k , and 2. a secondary mesh M, with N k := # F k (quantities linked to the secondary mesh are tagged by a tilde), on the domain Ω. Note that, despite their names, these meshes can be unrelated; 4 we refer to [23,p. 271,Definition 2] for the definition of a mesh.
We then associate to every k-form occurring in (1.1) a k-co-chain in either M or M, i.e., a mapping F k → R (or F k → R). Thus, every k-co-chain can be identified with its coefficient vector ∈ C N k (or ∈ C N k ), tagged with · (or · ) in the following. In other words, we adopt the following discrete models 5 for the unknown differential forms of (1.1): The entries of each coefficient vector, for example u, can be regarded as an approximation of integrals of the corresponding form, u ∈ Λ l−1 (Ω), over the respective ((l − 1)-dimensional) facets of the (primary) mesh. These integrals serve as degrees of freedom. The representation (2.1) has to be supplemented by introducing matrix representations for all linear operators of (1.1). For instance, the (primary/secondary) exterior derivative matrices D k−1 ∈ {−1, 0, 1} N k ,N k−1 and D k−1 ∈ {−1, 0, 1} N k , N k−1 are the incidence matrices between oriented kand (k − 1)-dimensional facets of the primary/secondary mesh, which allow to write algebraic equations linking primary and secondary co-chain-based vectors of degrees of freedom (2.1) through the discrete equilibrium equations 3 We refer to [23, p. 271, Definition 1] for the conditions that a sequence of spaces of differential forms has to satisfy to be discretizable. 4 Given a numerical method expressed by co-chain calculus, some of its degrees of freedom may be represented on the primary mesh and others on the secondary mesh. In this case, a bijective relationship between the two types of unknowns is needed, which can be achieved by using a secondary mesh dual to the primary mesh. This leads to numerical schemes fitting the framework of DEC. 5 Note that the identification of discrete spaces with mesh facets is restricted to lowest-order methods (for example, lowest-order Whitney forms). In fact, co-chain calculus does not generalize easily to higher-order spaces. The discrete constitutive equations can be written as have to fulfill the following three requirements [24, p. 254]: 6 They are the essential building blocks of discrete Hodge operators. 7 (Ω), respectively. Pairing matrices need to fulfill the algebraic relationship [24, p. 254, (4)] for any l ∈ {1, . . . , n}, m := n − l in R n , n ∈ N * . 3. The third requirement is a discrete counterpart of the integration by parts formula [24, p. 254, (6)] the numbers of (l − 1)-and m-dimensional primary/secondary mesh facets ⊂ Γ, select the degrees of freedom on Γ, while the boundary pairing matrix is only defined on Γ. Next we elaborate the discrete formulations on Ω for both ways of Sect. 1.4 to state the transmission problem of Sect. 1.1.

Option I: discrete primary formulation
We can now deduce an expression for the discrete form of (1.7), also given in [23, p. 276, Primary elimination, (16)]: ψ 0 ∈ C Nm+1 is a known vector determined by integrals of the source (m + 1)-form ψ 0 over (m + 1)dimensional facets of the secondary mesh. Thus, the secondary mesh only comes into play through the right-hand side of (2.5) and does not matter for the ultimate system matrix. Conversely, l−1 is an abstract boundary-energy term related to the Dirichlet-to-Neumann (DtN) operator. We clarify its meaning in Sect. 4.1.
can also be invertible matrices scaled by a negative or imaginary scalar. Imaginary γ is indeed the case of (1.2b). 7 While the computation of the entries of matrices M l α , M l−1 γ is done differently in FEEC and DEC, as discussed in [23, p. 277, Section 5], the discrete matrix forms of co-chain calculus only need to respect the few algebraic requirements listed here, independent of the details of the approximation. Hence, the actual construction of such matrices is only addressed when the framework of co-chain calculus is specialized into a numerical method (e.g., Sect. 6.1). 8 Note that pairing matrices, like mass matrices, depend on the adopted numerical method (i.e., on the choice of the discrete basis functions that approximate the differential forms).

Option II: discrete hybrid formulation
We can also discretize the variational form (1.9) in the general framework of co-chain calculus: where j ∈ C Nm is the coefficient vector whose entries are related to integrals of j ∈ Λ m (Ω) over mdimensional facets of the secondary mesh. Indeed, similarly to the derivation of (1.9) from (1.7), we

Trefftz methods
Trefftz methods seek to approximate the solution of constant-coefficient boundary value problems in (unbounded) domains by means of finitely-many global basis functions that solve the equations of the problem exactly and satisfy suitable conditions at infinity [26]. Hence, the main feature that characterizes a Trefftz method is its own discrete Trefftz space, and the functional form of the corresponding discrete basis functions leads to different types of Trefftz methods. These functions can be: • Plane waves or (generalized) harmonic polynomials, which constitute the most common choice [26].
• Fundamental solutions of the PDE-operators, which yield Trefftz functions with a central singularity placed outside the domain of approximation [30]. • Related to the above, fields generated by point sources solving homogeneous equations; 9 we then get the method of auxiliary sources 10 [44].
In spite of this diversity, all Trefftz methods share a desirable feature and a drawback. The former is the exponential convergence of their approximation error when the unknown possesses an analytic extension beyond the Trefftz approximation domain. 11 This is proven formally for 2D Poisson's equation in [15, p. 3, Section 2.2] and with numerical evidence for 2D Helmholtz equation and Maxwell's equations in [16] and [18], respectively.
The drawback is that, as exact solutions of a PDE are global functions, simple choices for a basis of a Trefftz space may be almost linearly dependent. Hence, Trefftz basis functions typically lead to ill-conditioned dense matrices. Stability may then be an issue.
To overcome the near-linear interdependence, Trefftz basis functions can be made orthogonal by a change of basis [6] or by choosing them orthogonal in the first place. However, the impact of illconditioning is limited due to the low number of degrees of freedom required for Trefftz methods, given their exponential convergence.
Related to the above is the need of heuristic rules to build the discrete Trefftz spaces when the unknown field is difficult to approximate: exponential convergence of Trefftz methods fails when, e.g., the field features singularities. Then, large Trefftz spaces may be needed and instability can have such a large impact that the numerical solution becomes useless.
Coupling a Trefftz method with a volume-mesh-based method can be a way to overcome this issue. As a matter of fact, by truncating the mesh at an artificial boundary that does not coincide with any physical discontinuity, Trefftz methods can be applied to a region where the unknown is sufficiently easy to approximate that the details of choosing the Trefftz space do not have much impact and exponential convergence is preserved. 12

Trefftz spaces for transmission problems
To capture the exterior problem (1.1) subject to a suitable condition at infinity (Option I, Sect. 1.4.1), we want to use functions that belong to the Trefftz space 13 v smooth and satisfying a suitable condition at infinity .
v smooth and satisfying a suitable condition at infinity .
Note that the requirement d ( v) = 0 is just a gauge condition 14 to ensure that we recover a unique potential in Ω + .

Coupling strategy
Given co-chain calculus inside Ω (Sect. 2) and a Trefftz method in the complement Ω + (Sect. 3), we can now devise a general discretization of coupled problems in differential forms (1.1) for both transmission problems of Sect. 1.4. We first derive the corresponding coupled systems in variational form in the continuous case. Afterwards, we rely on the co-chain calculus of Sect. 2 and the Trefftz spaces of Sect. 3 to accomplish discretization.
Our approach weakly imposes transmission conditions across Γ. In a sense, it is inspired by the socalled three-field variational formulation proposed for the sake of domain decomposition in [13]. Owing to the Trefftz approximation, our coupling approach can also be regarded as a generalization of the DtNbased coupling from [15, p.  On one hand, we plug (4.1a) into (1.7) directly. On the other hand, we impose (4.1b) weakly using as test functions the so-called co-normal traces t (−1) l α dv of Trefftz functions v ∈ T l−1 (Ω + ) of (3.1a). This yields the coupled variational problem Another expression for (4.2) can be obtained by noticing that which holds by integration by parts because of v ∈ T l−1 (Ω + ). We then obtain the coupled problem in variational form Seek The bottom-right integral in (4.4) expresses the "energy" associated with (1.5) in Ω + . Discretization Similarly to the discrete counterpart (2.5) of the weak formulation (1.7), we can write the discrete version of (4.2) and (4.4) in abstract algebraic form using, on top of the terms of (2.5), the other following terms: • We call P Γ ∈ C N bnd m ,N+ co-normal trace projection matrix, with N + dimension of the discrete Trefftz space T l−1 h (Ω + ) ⊂ T l−1 (Ω + ). Comparing (4.2) and (4.5), it is clear that P Γ is the matrix representation of a mapping T l−1 (Ω + ) → m-co-chain on M Γ .
• v ∈ C N+ is the vector of coefficients with respect to a basis of the discrete Trefftz space T l−1 h (Ω + ).
• M + ∈ C N+,N+ is the energy matrix in Ω + : To arrive at an alternative expression of (4.5) involving M + , we note that the total number of degrees of freedom of the Trefftz discretization, N + , is generally low because, under certain conditions, Trefftz methods enjoy exponential convergence (see Sect. 3). Thus, M + can easily be inverted by Gaussian elimination, and we can write the Schur complement of (4.5): We can now compare the left-hand side of (4.7) with the original discretization (2.5) and write In this way, M + provides a concrete expression for M l−1 β,Γ . Note that (4.8) is also the discrete version of the integration by parts formula (4.3).

Option II: hybrid formulation
For the hybrid problem (1.1) and (1.8), transmission conditions between Ω and Ω + are To impose these conditions, similarly to what was done for (1.9), we can write that α dv = (−1) l j in Ω + . (4.4) (with γ = 0 in Ω + ) can therefore be rewritten as (4.10) Note that j, ι belong to the Trefftz space T m (Ω + ), which correspond to T l−1 (Ω + ) of (3.1a) after taking co-normal traces t (−1) l α dv of its functions v. Let us now take π, τ ∈ T m−1 (Ω + ) of (3.1b) such that, in Ω + , j = dπ and ι = dτ . System (4.10) then becomes 15 Seek u ∈ Λ l−1 (Ω), π ∈ T m−1 (Ω + ): for the Trefftz space T m−1 (Ω + ) defined in (3.1b). In (4.11), we can finally replace the integral in Ω + with an integral on Γ, similarly to (4.3). Discretization Similarly to (4.4) and (4.5), the variational formulation (4.11) can be linked to a linear system of equations in a discrete setting: Note the occurrence of an additional boundary incidence matrix, D m−1 Also note that, abusing the notation of (4.5), we kept the symbols P Γ ∈ C N bnd m−1 ,N+ and M + ∈ C N+,N+ for the co-normal trace projection matrix and the energy matrix, even though here they have a different construction (and P Γ ∈ C N bnd m ,N+ in (4.5)). The next section will provide some details on the linear systems (4.5) and (4.12) for the specific case of the eddy-current equations from Sect. 1.2.

Specialization
Here, for the eddy-current problem introduced in Sect. 1.2, we present • transmission conditions across Γ, which separates the co-chain-based and Trefftz discretizations, • a basis for the Trefftz spaces of Sect. 3, and • the construction of the matrix of Sect. 4 involving only 16 Trefftz degrees of freedom, i.e., M + .

Option I: discrete primary formulation
Analogously to (4.1), the appropriate transmission conditions across Γ are the tangential continuities of A and ν ∇ × A. For the sake of simplicity, we also assume constant ν = ν + > 0 and σ = σ + > 0 in Ω + . Hence, given a spherical coordinate system in  • ∇ sph denotes the gradient in spherical coordinates and Y lm (θ, ϕ) the spherical harmonics [28, p. 108, (3.53)]. It can be shown that Φ lm , Ψ lm do not depend on r despite its presence in their definitions (5.2b) and (5.2c). Finally, M + of (4.5) has the following entries: where v j , j = 1, . . . , N + , are the Trefftz basis functions (5.1) of T 1 h (Ω + ), n the normal vector on Γ, and the second equality holds because of Trefftz functions in T 1 (Ω + ) being exact solutions of the homogeneous problem and integration by parts.

Option II: discrete hybrid formulation
Now as transmission conditions across Γ we have to impose matching tangential components of ν ∇ × A and ∇φ, and matching normal components of ∇ × A and ν −1 + ∇φ (see (4.9)). For the sake of simplicity, we also assume constant ν = ν + > 0 in Ω + .
To discretize the Trefftz space T 0 (Ω + ) of (3.1b), we can then choose the Trefftz basis functions Hence, for the special transmission problem (1.10) and (1.11) in Ω + , the entries of M + in (5.3) are

Numerical simulation
We again consider the eddy-current problem (1.2) in Ω and (1.10) in Ω + with the radiation condition (1.11). This problem is solved for the magnetic vector potential A in Ω using the cell method (Sect. 6.1) as co-chain-based discretization and the magnetic scalar potential φ in Ω + using Trefftz functions (5.4), hence by the coupling based on the hybrid formulation as discussed in Sect. 4.2 (Option II), which is specialized for the case of the cell method in Sect. 6.2.

Cell method
The Cell Method (from now on, CM), established in [40], relies on a pair of meshes for the spatial discretization of boundary value problems: one mesh being the dual of the other. Specimens of such dual mesh can be built by a Delaunay-Voronoi subdivision, as proposed in [41], whereas barycentric dual meshes are used in [4,20,33,34] for computational advantages [40, p. 13].
CM is an example of a method that can be fitted into the co-chain-based framework presented in Sect. 2. In particular, CM degrees of freedom are integrals of fields on entities of the primary and dual meshes. 18 In the context of electromagnetics, where CM has long been used, examples are fluxes of the magnetic flux density on primary faces or line integrals of the magnetic field on dual edges [41]. This implies 19 that the CM discretization, as a particular case of co-chain calculus, replaces differential operators with incidence matrices. In this way, equilibrium equations involving such operators (e.g., (1.2a)) can be enforced exactly. Basis functions are then required to interpolate fields locally in terms of these integral degrees of freedom and approximate the material laws (constitutive equations, like (1.2b)). With the choice of Whitney basis functions, the matrices of FEEC can be recovered [40,41].
Vector fields can also be interpolated by the piecewise-constant vector basis functions defined in [21] for tetrahedral cells. Specifically, positive-definite mass matrices in the CM discretization scheme can be obtained, for instance, by the so-called energy approach described in [21,36], which makes it possible to exactly reconstruct the energy in a domain Ω with a piecewise-constant field.
Here we give examples of such mass matrices for the eddy-current model of Sect. 1.2 (which corresponds to the case n = 3, l = 2, m = 1). We recognize the general parameters α and γ as ν (reluctivity) and −ıωσ (σ is conductivity), respectively, so that the general mass matrices of Sect. 2, M 2 α and M 1 γ become M ν and −ıωM σ , where M ν ∈ R N faces ,N faces and M σ ∈ R N edges ,N edges (given N faces number of faces and N edges number of edges of the primary mesh M) are defined as: Here, b f i ∈ L 2 (Ω) is the vector face basis function related to the i-th face of the tetrahedral mesh M and b e j ∈ L 2 (Ω) is the vector edge basis function related to the j-th edge of M.

CM-Trefftz coupling
For the eddy-current model introduced in Sect. 1.2, assuming constant ν = ν + > 0 in Ω + , the projection matrix P Γ ∈ C N bnd edges ,N+ (linking the CM and Trefftz discretizations) can be constructed as for v j , j = 1, . . . , N + , chosen Trefftz basis functions of T 1 h (Ω + ), like (5.1) (Option I), where i , i = 1, . . . , N bnd edges , is an edge of M on Γ and n its normal vector. Conversely, for v j chosen Trefftz basis functions of T 0 h (Ω + ), like (5.4) (Option II), the entries of P Γ are with τ the unit tangential vector of i , x i,1 , x i,2 the endpoints of i (which correspond to centroids of faces of M on Γ, using a barycentric dual mesh as M), and i,1 , i,2 being = 1 or −1 depending on the ordering of x i,1 , x i,2 as endpoints of i . Note that, because of the last equality in (6.2b), P Γ for Option II can be reorganized as P Γ ∈ C N bnd nodes ,N+ (degrees of freedom on the boundary nodes of M). Option II, i.e., the discrete system (4.12) with matrices P Γ of (6.2b) and M + of (5.5), is used for the numerical results of Sect. 6.3 for the eddy-current model of Sect. 1.2. Here we specialize the discrete system (4.12) for that eddy-current model, using CM and a Trefftz method for discretization: which has a similar structure to the system described in [36] between CM and the boundary element method, being both complex Hermitian. (6.3) contains the following terms: • D 1 ∈ {−1, 0, 1} N faces ,N edges is the incidence matrix from edges to faces of M in Ω (discrete curl ). • Mass matrices M ν ∈ R N faces ,N faces and M σ ∈ R N edges ,N edges for CM were introduced in (6.1).  [35], which makes it possible to construct a divergence-free source current vector within Ω 0 . • The energy matrix M + ∈ C N+,N+ of (5.5) and the vector of Trefftz coefficients v + ∈ C N+ need no further explanations.

Results
The CM-Trefftz coupling is implemented in a vectorized MATLAB code and tested with a realistic 3D eddy-current problem, which consists of a cylindrical inductor (made of a conductive material) excited by an alternating current-driven coil. To assess the accuracy of the 3D CM-Trefftz method, we also solved this axisymmetric model (on a radial cross-section) with a third-order 2D FEM approximation, computed with COMSOL Multiphysics 5.5, and used that solution as reference. All our numerical tests were run on a laptop with an Intel Core i7-6920HQ processor (2.90 GHz clock frequency) and 16 GB of RAM. Figure 1 shows a sketch of the axisymmetric inductor (radial cross-section). The cylindrical inductor core Ω (10 mm radius, 40 mm high) is made of a conductive material with 25 MS m −1 conductivity and 2 relative permeability. The inductor coil Ω 0 (square cross-section, 4 mm side, centered at r = 15 mm, z = 0) carries a 16 A current at 200 Hz frequency, which excites eddy currents in Ω . Finally, the surrounding domain Ω is a sphere with 40 mm radius, which includes an air domain (null conductivity, 1 relative permeability). Ω is discretized by CM (with piecewise-constant edge and face elements) using a tetrahedral mesh.
With a Trefftz method, meshing the exterior domain Ω + is not required, since the magnetic field there is approximated by the gradients of (5.4). This Trefftz discretization is coupled with CM by suitable transmission conditions at the interface Γ.
The eddy-current distribution in Ω , the magnetic flux density distribution in Ω, i.e., the real and imaginary parts of the rand z-field components, and the aforementioned accurate numerical solution are computed with third-order 2D FEM as reference. To account for the truncation error, the FEM  table 1) are employed to check the convergence of the CM-Trefftz method. We define mesh size as the maximum diameter of spheres circumscribing the tetrahedrons of a mesh. The same fixed Trefftz space is used with different mesh refinements: one expansion (5.4) of fifth order and centered at the origin is used. In fact, it is observed that increasing the order does not provide an improvement in the convergence of the error (i.e., the CM error dominates) and the local magnetic field reconstruction is very good even with a low order.
For the sake of comparison, the following relative errors (in L 2 -norm) for domains Ω and Ω are defined: where J h , B h are the approximated field distributions computed by the 3D CM-Trefftz coupling and B, J the reference field distributions computed by third-order 2D FEM. Table 1 shows the assembling and solving times needed for each mesh. The assembling time encompasses the generation of both incidence and mass matrices, the generation of the right-hand side vector (from source currents), and the assembly of the final matrix system (6.3). The solving time corresponds to the solution of the final matrix system, which is carried out by the BiCGSTAB iterative solver [9, p. 24, Section 2.3.8] with a SSOR preconditioner [9, p. 37, Section 3.3] (given a prescribed tolerance of 10 −10 ). Figures 2 and 3 show the relative errors (6.4), computed for the same mesh refinements as above. It can be observed that the 3D CM-Trefftz coupling attains a first-order convergence rate in the L 2 -norm. In particular, Fig. 2 proves that a very good accuracy in the eddy-current density reconstruction within the core is obtained by using the proposed hybrid approach.   In the following, field lines with 3D CM-Trefftz are considered on the symmetry plane y = 0 (but other radial cross sections can be considered). Moreover, comparisons between field profiles are all taken with the intermediate mesh refinement of h = 2.76 mm. Figure 4 shows the real and imaginary parts of the azimuthal eddy-current density J θ along line A-B. The CM-Trefftz plot is obtained by locally interpolating the electric field (proportional to the eddy-current density) with piecewise-constant edge basis functions. Maximum discrepancies from the third-order 2D FEM solution for the real and imaginary parts of J θ are 18.36% and 13.14%, respectively.
The same comparison is then extended to the magnetic flux density in the domain with CM discretization. In this case, CM-Trefftz plots are obtained by locally interpolating the magnetic flux density with piecewise-constant face basis functions. Figures 5 and 6 show the real and imaginary parts of the xand z-components along line C-D. Maximum discrepancies are 26.06%, 20% for the real and imaginary part of B x , and 8.48%, 11.77% for the real and imaginary part of B z .
In the exterior domain, with Trefftz discretization, the hybrid approach proves to be highly accurate, even though a single expansion to the fifth order is used. For plots shown in Figs. 7 and 8, with field profiles computed along line E-F, all discrepancies are under 1%: 0.81%, 0.86% for the real and imaginary part of B x , and 0.73%, 0.95% for the real and imaginary part of B z .

Conclusions
The authors are not aware of any prior work addressing the coupling of Trefftz methods with a general framework for numerical schemes based on volume meshes, like co-chain calculus. We explain for the first time how to formulate an interface operator for this framework that takes into account the exterior problem (see (4.8)). The particular case of coupling the cell method with Trefftz basis functions is also presented here for the first time.  The other particular case of Trefftz methods coupled with the finite element method as mesh-based discretization can be seen in [18], which also includes numerical results. There one can see that, compared to other hybrid approaches that rely on the boundary element method, using a Trefftz method enjoys advantages [29, p. 51], namely • a simpler assembly process, as there are no singular integrals, and • an exponentially-convergent approximation error when the unknown possesses an analytic extension beyond the Trefftz approximation domain (see Sect. 3). Admittedly, Trefftz methods can suffer from ill-conditioning. However, the impact of ill-conditioning is mitigated due to the low number of degrees of freedom required for Trefftz methods, given their