Conservation of Forces and Total Work at the Interface Using the Internodes Method

The Internodes method is a general purpose method to deal with non-conforming discretizations of partial differential equations on 2D and 3D regions partitioned into disjoint subdomains. In this paper we are interested in measuring how much the Internodes method is conservative across the interface. If hp-fem discretizations are employed, we prove that both the total force and total work generated by the numerical solution at the interface of the decomposition vanish in an optimal way when the mesh size tends to zero, i.e., like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {O}(h^{p})$\end{document}O(hp), where p is the local polynomial degree and h the mesh-size. This is the same as the error decay in the H1-broken norm. We observe that the conservation properties of a method are intrinsic to the method itself because they depend on the way the interface conditions are enforced rather then on the problem we are called to approximate. For this reason, in this paper, we focus on second-order elliptic PDEs, although we use the terminology (of forces and works) proper of linear elasticity. Two and three dimensional numerical experiments corroborate the theoretical findings, also by comparing Internodes with Mortar and WACA methods.


Introduction
We are interested in the approximation of partial differential equations by non-conforming domain decomposition methods, more specifically in the conservation properties at the interface between non-overlapping subdomains. The non-conformity at the interface can be induced by different independent discretizations inside the subdomains, as well as by non-watertight interfaces.
Is the traction and the work at exact equilibrium across the interface? If not, how well it is, and what is the dependence on the mesh refinement? And maybe even more importantly, which quantities may-be expected to be conserved at the interface?
Among the non-conforming coupling methods, the most widely used and studied is the Mortar method (far from be exhaustive see, e.g., [3-5, 7, 9, 29]), which is optimal in terms of convergence rates and is based on a single L 2 projection operator; as we will see throughout the paper, the Mortar method also conserves traction forces and works along the interface; see also [14] for a discussion in the case of fluid-structure interation problems.
Differently from the Mortar method, the Internodes (INTERpolation for NOnconforming DEcompositionS) method is based on two intergrid interpolation operators, one for transferring the Dirichlet trace across the interface, the others for transferring the Neumann trace, i.e., the fluxes or the stresses. The Internodes method has been proposed in [11] in the context of second order elliptic differential problems, and it has been applied successfully to fluid-structure interaction problems [13], Navier-Stokes equations [11], Stokes-Darcy coupling [18]. Its analysis for 2D and 3D second-order elliptic equations has been carried out in [16] and its generalization to decompositions with more than two subdomains is presented in [16,19]. The Internodes method has also been applied in the context of Isogeometric Analysis [15] to deal with non-conforming multi-patch geometries. For what concerns the theoretical analysis, it has been proved in [16] that, when the mesh sizes h 1 and h 2 of the two subdomains vanish with the same rate and the local polynomial degree p (the same in both subdomain) grows up, the Internodes method features the same convergence order of the Mortar method in the H 1 -broken norm error. More precisely, it holds that the global H 1 -broken norm error behaves like O(h p ), with h = max{h 1 , h 2 }. When this happens, it is said that the convergence order of the multidomain method is optimal. The spectral properties of the algebraic form of the Internodes method have been investigated in [11,Sect. 6] and compared with those of the Mortar method. Numerical results show that the extreme moduli of the eigenvalues of the global Internodes matrix, associated with the discretization of the Laplace operator, behave like the analogous quantities related to the Mortar method.
We point out that the Internodes method works with two independent interpolation matrices, and this choice guarantees to keep the optimal convergence order of the local discretizations inside the two subdomains. On the contrary, the consistent interpolating approach discussed in [14], that coincides with the pointwise interpolation method analysed in [6], features a degradation of the solution of the coupled problem with respect to the mesh-size and to the polynomial degree, although it is conservative across the interface.
Alternatives are Nietsche [1] or discontinuous Galerkin approaches, in which jumps between domains are evaluated, and that needs interpolation or projection operators between subdomains for the computation of integrals of quantities of both domains simultaneously. The Internodes method needs an interpolation operator between the subdomains but there is no need of crossed integrals.
In this paper we show that the Internodes method conserves both traction forces and works across the interface, in the sense that the total force and the total work at the interface vanish at least like the global broken-norm error when the discretizations in both subdomains are refined.
The conservation properties of some coupling methods have been reviewed in [10]. In the present paper, we try to give a clear characterisation of what conservative method means by defining the quantities that a coupling method should conserve (asymptotically or exactly); then, on one hand we provide a mathematical proof of the conservation properties of the Internodes method, on the other hand we yield numerical experiments showing the conservation properties of the Mortar method, the Internodes method, and the weighted average continuity approach (WACA) proposed in [10].
In [22] the authors addressed the problem of conserving the energy across an interface in the framework of non-conforming interfaces between a finite volume and FEM; they consider that the energy is conserved if the work is equal on both sides of the interface.
We observe that the conservation properties of a method are intrinsic to the method itself because they depend on the way the interface conditions are enforced rather than on the problem we are called to approximate. For this reason we focus on second order elliptic partial differential equations (PDEs) and after decompositon of the computational domain into two subregions, we analyse the conservation properties across the common interface; this interface can be either conforming or non-conforming depending on the local space discretisations and on the meshes. In the continuous settings, the total force and the total work across the interface are zero, both when computed as integrals (in strong form) on the interface and as sum of residuals (weak form); indeed they coincide thanks to the Green formula. When using classical conforming methods based on the Galerkin projection [26,28], the sums of residuals are null while the strong forms asymptotically vanish when the mesh-size tends to zero.
In case of non-conforming methods, one has to first establish how to compute the total force and work, since the interface on the two sides may be different (geometric nonconforming, see, e.g., the characterization given in [11]). In this case the integral forms and the residual ones do not coincide anymore. The Mortar method guarantees that the residual forms of total work and total force are zero, however the strong ones are "only" optimally convergent, but not identically null. By optimal convergence we mean here the H 1 -brokennorm convergence order that behaves like the worst best approximation error inside the two subdomains and that is generally used to measure the approximation properties of non-conforming methods.
In Section 4 we prove that Internodes provides optimal convergence for both versions (strong and residual-based) of the total force and work. We have performed two and three dimensional experiments using geometrical conforming and non-conforming decompositions and different polynomial degrees, by using finite element or spectral element methods. We have compared Mortar, Internodes, and WACA methods. Mortar and Internodes provide optimal order convergence (when not exact, as explained above), whilst WACA allows for exact residual versions of the total quantities, but at most for first order convergence for the strong forms, even for high order approximations.
The paper is structured as follows. In Section 2, we briefly recall the mathematical setting of elliptic PDE's in a two subdomains framework in either 2D or 3D domains and, by inheriting the terminology from linear elasticity, we define the quantitites that we want to analyse in order to evaluate the conservation properties of a multidomain method, i.e., the total force and the total work. In Section 3 we recall the Internodes method and its approximation properties. Section 4 is devoted to the analysis of the conservation properties of the Internodes method and in Section 5 we provide numerical evidence to the conservation and approximation properties of the Internodes and Mortar methods by reporting numerical results in both 2D and 3D domains.
The multidomain formulation of (1) reads (see, e.g., [26,28]): for k = 1, 2 look for the functions u k defined on Ω k such that where n k is the outward unit normal vector to ∂Ω k , while ∂Ω k,D = ∂Ω D ∩ ∂Ω k and ∂Ω k,N = ∂Ω N ∩ ∂Ω k are the Dirichlet and Neumann boundaries, respectively, restricted to the domain Ω k . Condition (2) 2 enforces the continuity of the solution across the interface Γ , while (2) 3 enforces the balance of the normal derivatives on Γ . To write the weak form of (2) we define the spaces Λ is the space of traces of the functions of V on the interface.
Provided that f ∈ L 2 (Ω), g N ∈ L 2 (∂Ω N ), and ν, γ ∈ L ∞ (Ω), the weak form of (2) reads: for k = 1, 2 find u k ∈ V k such that where the integrals in the last equation must be interpreted as dualities between the trace space Λ and its dual space Λ . The classical abstract form of (3) reads: where a k (u k , v k ) = Ω k ν∇u k · ∇v k + γ u k v k , F k (v k ) = Ω k f v k + ∂Ω k,N g N v k and R k : Λ → V k denotes any linear and continuous lifting operator from the interface to the domain Ω k . It is well known (see, e.g. [26]) that problem (4) admits a unique solution and that u 1 and u 2 are the restrictions to Ω 1 and Ω 2 , respectively, of the weak solution of (1).
The conditions (3) 3 (that should be interpreted as duality when the normal derivatives are not sufficiently regular) and (4) 3 are the weak counterpart of (2) 3 and, by choosing the test function μ in a suitable way, they express two fundamental balance principles at the interface, that we are going to describe below. We make the following assumptions that we comment at the end of the present section, see Remark 1.

Assumption 1
Let us assume that Neumann boundary conditions are imposed on that part of ∂Ω that matches the interface Γ (i.e., we assume that ∂Γ ∩ ∂Ω D = ∅).
Let Assumption 1 be satisfied. By adopting the terminology of linear elasticity (in which the Laplace operator is replaced by the divergence of the Cauchy stress tensor), if we choose μ ≡ 1 in (3) 3 we obtain a sort of "balance of forces" at the interface, that reads where T F stands for Total Force. This terminology is inherited from linear elasticity, where the normal derivative of the solution on the boundary is the normal component of the stress tensor and expresses a traction or a normal force to the boundary.
If instead we choose μ = u 1|Γ = u 2|Γ (in virtue of (3) 2 ), we obtain a sort of "null total work" on Γ , i.e, here T W stands for Total Work. Also here we refer to linear elasticity, where the product ν ∂u ∂n u is replaced by the product between the normal stress and the displacement, thus giving a work.
Finally, we notice that by setting λ = u 1|Γ = u 2|Γ in (6) and by applying the first Green formula, the identities (5) and (6) can be equivalently written as In fact, they are sum of residuals and are particular instances of the equation (4) 3 with μ = 1 and μ = u |Γ , respectively. Formulas (5) and (6) are what we call strong forms of the total force and total work, respectively, while (7) and (8) are the corresponding residual or weak forms. The balance of forces (5), or (7), and the null total work (6), or (8), express the conservation for the multidomain problem (3).
When the multidomain problem (3) is approximated (by either a conforming or a non-conforming method) it is desiderable that the discrete solution satisfies the discrete counterpart of (5)-(8) up to a term that vanishes with a certain order q with respect to the discretization parameters, e.g. the mesh size of the discretization. If that happens, we say that the multidomain approach is conservative of order q with respect to the mesh size. We refer to Section 4, Definition 1 for a formal definition of this concept.
In this paper, we show that when we combine Internodes with Finite Elements or Spectral Elements, we obtain a conservative multidomain discrete approach and that the order of conservation equals the order of the broken-norm error (with respect to the mesh size).
We notice that conservation of a method is intrinsic to the method itself because conservation depends on the way the interface conditions are enforced rather then on the problem we are called to approximate. Thus, the analysis provided here can be extended to other PDEs for which the conservation at the interface is meaningful.
Remark 1 Assumption 1 is made for two reasons both connected with the theoretical analysis. First of all we notice that, if we set Dirichlet conditions on the part of ∂Ω that matches the interface Γ , then μ should belong to Λ = H 1/2 00 (Λ) and we could not choose μ ≡ 1 in (3) 3 . This would imply that the conservation of forces could not be brought back to the interface condition (3) 3 . The second reason is related to the convergence analysis of the Internodes method that, at the moment, is available provided that Assumption 1 is satisfied.
However, we can infer from the numerical results of Section 5 that the conservation properties of the Internodes method are kept even when Assumption 1 is not satisfied.

Internodes for hp-fem Discretization
We sketch here the idea of Internodes and we refer to [11,16,19] for an exahustive presentation of the method for what concerns theoretical properties, algebraic formulation, and algorithmic aspects.
We consider two a-priori independent families of triangulations T 1,h in Ω 1 and T 2,h in Ω 2 , respectively, characterized by different mesh-sizes h 1 and h 2 . This means that the meshes in Ω 1 and in Ω 2 can be non-conforming on Γ . The elements can be either simplices (triangles if d = 2 or tetrahedra if d = 3) or quads (i.e. quadrilaterals if d = 2 or hexahedra if d = 3). Moreover, different polynomial degrees p 1 and p 2 can be used to define the finite element spaces on T 1,h and T 2,h . Inside each subdomain Ω k we assume that the triangulations T k,h are regular and quasi-uniform [25,Chapter 3]. Moreover, we assume that they are affine when simplices are considered. We denote by Γ 1 and Γ 2 the internal boundaries of Ω 1 and Ω 2 , respectively, induced by the triangulations T 1,h and T 2,h . If Γ is flat, then Γ 1 = Γ 2 = Γ , otherwise Γ 1 and Γ 2 can be different. The finite element approximation spaces (for k = 1, 2) are In order to exchange information between the two independent grids on the interface Γ , we introduce two independent interpolation operators: that are going to interpolate the derivatives and traces from one interface to the other one. These arguments apply both in two and three dimensions. When Γ 1 and Γ 2 coincide, then Π 12 and Π 21 are the classical Lagrange interpolation operators (see [11,16]), otherwise Π 12 and Π 21 can be defined as the Rescaled Localized Radial Basis Function (RL-RBF) interpolation operators introduced in formula (3.1) of [12] (see also [19,Sect. 2
The Internodes method applied to (1) reads: find where is the so-called residual at the interface Γ k , and it is defined starting from the solutions u k,h as follows: -set the real values where i ∈ Y k,h (that is the function in X k,h that is null at all nodes of T k,h not belonging to Γ k and that coincides with μ (k) i on Γ k ), -define the functions where {Φ We remark that the expansion (12) with respect to the dual basis is not suitable to apply the Lagrange interpolation, but in order to interpolate r k,h we have to write it with respect to the primal basis. This is possible because Y k,h and Y k,h are the same algebraic space [8]. The matrix that realizes the change from the primal basis to the dual one is the interface mass matrix M Γ k whose entries are conversely, M −1 Γ k realizes the change from the dual to the primal expansion. Denoting by r (k) Γ the array whose entries are the values r k,i , we define the array whose entries are denoted by z k,j for j = 1, . . . , n k . Then the primal expansion of the residual function r k,h reads and we are going to apply the interpolation to it.
Remark 2 If Assumption 1 are not satisfied and Dirichlet conditions are imposed on the part of ∂Ω matching the boundary of the interface Γ , then formula (11) must be replaced by where Formula (14) is justified as follows. When a Lagrange basis function μ (k) i is associated with a node x i belonging to the boundary of the interface, its extension R k μ (k) i is indeed also not null on a part of the external boundary ∂Ω k \ Γ k . Thus, by integrating by parts a k (u k,h , R k μ (k) i ) and bearing in mind that r k,h must be an approximation of the normal derivative of u k,h on Γ k , we obtain (14).

Algebraic Form of Internodes
Using standard notations for finite elements, we denote by A (k) the stiffness matrix associated with the bilinear form a k and we decompose it into the 2 × 2 block-wise form Γ be the arrays of the d.o.f. in Ω k \ Γ k and on Γ k , respectively, while r (k) Γ is the array of the interface residual introduced before. Finally, let P 12 and P 21 the matrices associated with the interface interpolation operators Π 12 and Π 21 , respectively, introduced in (9).
The algebraic form of the Internodes method reads: find u (1) If the meshes are conforming at the interface, then the interpolation matrices P 12 and P 21 are in fact the identity matrix and M Γ 2 = M Γ 1 . It follows that the Internodes method (15) that is nothing else the algebraic form of a classical substructuring method on two subdomains (see [26,28]).

Accuracy of the Internodes Method
Under the assumptions that problem (1) is well posed (see, e.g., [25]), the convergence analysis of the Internodes method with respect to the mesh sizes h 1 and h 2 is carried out in [16]. More precisely, let denote the Internodes solution and let u ∈ H 1 0,∂Ω D (Ω) be the solution of the weak monodomain formulation We define the broken-norm error Let u λ k be the the weak harmonic extension of λ to Ω k , i.e. the solution of the problem u λ k ∈ V k : a k (u λ k , v) = 0 ∀v ∈ V 0 k , u λ k = λ on Γ, and u k be the solution in Ω k of the elliptic problem with null trace on the interface, i.e.
Thanks to the linearity of the problem, when λ = u |Γ we have u k = u λ k + u k for k = 1, 2.
The following convergence result has been proved in [16,Theorems 8 and 10] in the case that Γ is a flat interface and the interpolation operators Π 12 and Π 21 are based on the Lagrange interpolation.

Theorem 1 Assume that the weak solution u of the monodomain elliptic problem belongs
Moreover, denoting by E the right-hand side of (18), Theorem 8 of [16] ensures that where λ k,h = u k,h | Γ are the traces on Γ k of the discrete (non-conforming) solution.
When the the ratio h 2 /h 1 is uniformly bounded from below and above (i.e., the two mesh-sizes h 2 and h 1 vanish with the same rate), this result guarantees that the Internodes method exhibits optimal accuracy, i.e., the broken-norm error behaves like the maximum between the energy-norm local errors in the two subdomains.
Indeed, by applying standard trace results on polyhedral domains (see, e.g., [21, Theorem 1.4.1]), we have that λ ∈ H s− 1 2 (Γ ) and r 2 ∈ H s− 3 2 (Γ ) and we can take σ = s − 1 2 and τ = s − 3 2 in (18). It follows that ρ k = min{s − 1 2 , p k + 1} and ζ k = min{s − 3 2 , p k + 1} and the minimum exponent among k − 1, ρ k − 1 2 , and ζ k + 1 2 is k − 1. Denoting by C = C(u) a positive constant depending on the exact solution u (and then also on the trace λ and the conormal derivative r 2 ) but independent of h 1 and h 2 , and recalling that we take h 2 /h 1 uniformly bounded from below and above, we have Thus, (18) and (19) can be summarized as Remark 3 Theorem 1 expresses the convergence order of the Internodes method with respect to the mesh sizes h 1 and h 2 , but not with respect to the local polynomial degree p. The convergence with respect to p is showed numerically in Section 5, its analysis is in progress.

Conservation Properties
Let us consider the balance conditions (5)- (8). We are interested in analyzing the discrete counterparts of T F , T W , T F a and T W a in order to measure the conservation properties of a non-conforming domain decomposition method like Internodes. As a matter of fact, it is not guaranteed that all such discrete counterparts are null, even when the meshes are conforming on Γ and the polynomial degrees coincide in the two subdomains.
Let E Γ k be the set of the edges of the elements in T k,h that belong to Γ k (see Fig. 2) and define Although the identities T F = T F a and T W = T W a are guaranteed at the continuous level, analogous identities are no longer valid at the discrete level. Indeed, by counter-integrating by parts both T F a,h and T W a,h , we obtain Fig. 2 The triangulation in Ω k and the sets ω k (light blue region), E k (blue edges), E k,N (green edges), and = w + · n + + w − · n − (following the standard notation of Discontinuous Galerkin methods), E k,N is the set of the edges of the elements in T k,h that belong to ∂Ω k,N , ω k is the set of elements in T k,h having an edge on Γ k , and E k is the set of the edges internal to Ω k (see Fig. 2). Notice that the finite element extension R k μ k,h is null on the blue edges that are not internal to ω k .

Algebraic Form of TF a,h and TW a,h
The terms T F a,h and T W a,h are strictly connected with the residual arrays r (k) Γ introduced in Section 3, and they can be easily computed by algebraic operations as follows [10]. We denote by 1 (k) the array of size n k whose entries are all equal to 1. It holds: where we have exploited definition (11), the fact that the Lagrange basis functions μ (k) i of Y k,h satisfy the unity partition property, i.e., n k i=1 μ (k) i ≡ 1, and that u k = R k (u k,h | Γ k ). It follows that T F a,h = (1 (1) ) T r (1) Γ + (1 (2) Γ .
Thus, it turns out very convenient to measure the conservation properties of a method by evaluating T F a,h and T W a,h by using these algebraic relations.

The Mortar Method
By adopting similar notations used to write the algebraic form (15) of the Internodes method, we write the algebraic counterpart of the Mortar method, which, instead to interpolate the trace and the normal derivative at the interface, is based on a projection process.
To this aim, let us denote by P the matrix implementing the projection of the trace from the interface Γ 1 to Γ 2 . Then we notice that Mortar is a symmetric method, in the sense that the operator used to move from Γ 2 to Γ 1 is the transposed operator of P .
The algebraic form of Mortar reads: Since the projection matrix P satisfies the property P 1 1 = 1 2 (this means that a constant function on Γ 1 is mapped on the same constant function on Γ 2 ), the identities immediately follow from (22) 2,3 .

The Internodes Method
On the contrary the Internodes method is not symmetric, since the two intergrid matrices P 21 and M Γ 1 P 21 M −1 Γ 2 are not one the transpose of the other, thus there is no way that (23) are exactly satisfied by Internodes, but in the next Section we will prove that both T F a,h ans T W a,h provided by Internodes go to zero like the broken-norm error when h 1 , h 2 tend to zero with the same rate. Clearly this result is weaker than (23), nevertheless we know that the broken-norm error of the Internodes method behaves like that of the Mortar method [16], and numerical results (see Section 5) show that T F h and T W h behave in the same manner for Internodes and Mortar methods.
Thus the question is: Are (23) necessary and sufficient conditions to guarantee that a method is conservative and, not less important, accurate?
We state that T F a,h = 0 and T W a,h = 0 alone do not guarantee that the coupling method one is using is convergent. This is the case of the Weighted Average Continuity Approach (WACA) proposed in [10, Sect. 3.5].

WACA
WACA can be formulated like (22), but with the matrix P replaced by M −1 Γ 2 S 2 P 21 S −1 1 M Γ 1 , where M Γ k are the interface mass matrices defined in (13), P 21 is the interpolation matrix associated with the interpolation operator Π 21 (see (9)), while S k is the lumped interface mass matrix 1 on Γ k . In view of the symmetry of WACA (like for Mortar), the identities (23) are satisfied.
However, when the discretization in almost one subdomain is based on Q p Spectral Element Methods with Numerical Integration (SEM-NI), in virtue of the fact the SEM-NI mass matrix is diagonal, we have that S k = M Γ k and P = P 21 , so that WACA method coincides with the so-called pointwise matching method that was presented in the seminal mortar paper [6, eqs. (3.5)-(3.7)] and that is notoriously sub-optimal, as proven in [6,Sect. 3.2] and numerically corroborated in [2] for spectral elements discretizations (see also Fig. 5).

Analysis of Conservation Properties
So far, set h = max{h 1 , h 2 } and let us give the following definition of conservation for a multidomain approach.

Analysis of Conservation Properties of the Internodes Method
The following theorems ensure that the discrete total force T F h and T F a,h and the discrete total work T W h and T W a,h converge to zero as optimally as the broken-norm error (20) and then, that the Internodes method is conservative of the same order of the broken-norm error.
In the whole Section make the following assumption.
Assumption 2 Let us assume that the discretisation of the interface Γ is geometrically conforming, i.e., ∪ e∈E Γ 1 e = ∪ e∈E Γ 2 e = Γ and that the intergrid operators Π 12 and Π 21 are based on the Lagrange interpolation. We also assume that there exist two positive constants c 1 and c 2 such that Remark 4 When the interfaces are not geometrically conforming, we have two difficulties in analysing the conservation properties: first we should be able to quantify the nonconformity of geometric type (and this is not often possible), second we should have convergence estimates of the Internodes method when RBF interpolation instead of Lagrange interpolation is used, and we do not have it. We notice that, for high polynomial degree p (typically when p ≥ 5), the RBF interpolation does not always exploit the same convergence order of the Lagrange interpolation and this can downgrade the accuracy of the Internodes method.

Theorem 2 Let Assumptions 1 and 2 be satisfied. If u ∈ H s (Ω) with s > 3/2, then there exists a positive constant C depending on the data (the domain and the coefficient functions) and on the exact solution u, but independent of h 1 and h 2 such that
where, for k = 1, 2, k = min(s, p k + 1), where p k is the polynomial degree in the domain Ω k .
Proof First we analyse T F h . Since ∪ e∈E Γ k e for k = 1, 2, and n 1 = −n 2 , we have From Cauchy-Schwarz inequality and the trace theorem for polygons and polyhedra [21,27], we have that The proof for T W h needs few more steps and is based on similar arguments: Thanks to the trace theorem for Sobolev spaces and the triangle inequality, we also have that where c is a positive constant depending on Ω k but independent of u. Then, recalling that λ = u| Γ and λ k,h = (u k,h )| Γ , exploiting again the trace inequality and the estimate (18), The proof is completed.
Now we are going to analyse the terms T F a,h and T W a,h . To this aim, we exploit the results proved in [16] after reformulating the problem (4) like a three fields problem as follows.
The following theorem is an immediate consequence of Theorems 2 and 3.
Theorem 4 Set = min{ 1 , 2 } and h = max{h 1 , h 2 }, under the assumptions of Theorem 2, the Internodes method is conservative up to order q = − 1, that is the order of the broken-norm error.
We can also formulate an upper bound for the forms B k defined in (21):

Numerical Results
The numerical results of this section corroborate the theoretical results proved above, that is the Internodes method is conservative of order q = − 1 (see Theorem 4).
We consider here both 2D and 3D test cases with either geometric conformity and nonconformity at the interface. The non-conforming geometric case is not covered by the theory but the numerical results are however consistent with the geometric conforming situation. In all cases we compare the quantities |T F h |, |T F a,h |, |T W h | and |T W a,h | with the broken-norm error (see (17)) and the L 2 error u − u h L 2 (Ω) . Moreover, for the first test case we compare Internodes with the Mortar method in terms of accuracy and conservation properties. We show that T F h and T W h behave like O(h p ) for both Mortar and Internodes (exactly as the broken-norm errors do); that |T F a,h | and |T W a,h | are O(h p ) for Internodes and null up to the machine precision for both Mortar and WACA.
In the whole section we set h = max{h 1 , h 2 }.

Test Case #1: d = 2, Flat Interface (Geometric Conformity)
Let us set the domain Ω = (0, 2)×(0, 1) and the coefficients ν = γ = 1; the right-hand side f and the boundary datum g N are such that the exact solution is u(x, y) = (2x + 1)(2y + 1) sin(xπ/2) sin(πy). Then we decompose Ω into the subdomains Ω 1 = (0, 1) × (0, 1) and Ω 2 = (1, 2) × (0, 1), so that the interface is Γ = {1} × (0, 1). Neumann boundary conditions are imposed on the horizontal sides of Ω and homogeneous Dirichlet boundary conditions are imposed elsewhere. In Fig. 3 we compare Internodes with the Mortar method, more precisely we consider P 1 fem in both subdomains of size h 1 = 1/(2k − 1) and h 2 = 1/k with k = 10, . . . , 20, so that the mesh is finer in the domain Ω 1 . The two approaches provide very similar errors (broken and L 2 -norms), also the quantities |T F h | and |T W h | are very similar and we have As anticipated in Section 4, the quantities |T F a,h | and |T W a,h | behave like O(h) for the Internodes method while they are about the machine precision for the the Mortar method.  In Fig. 4 we compare Internodes with Mortar as before, but now by taking p = p 1 = p 2 = 3, thus Q 3 -SEM-NI are employed in each subdomain. Numerical results show that The quantities |T F a,h | and |T W a,h | behave like O(h 3 ) for the Internodes method while, again, they are about the machine precision for the Mortar method. The different behaviour of T F a,h and T W a,h for the two methods strictly depends on the way the interface conditions are imposed at the interface, as we have explained in Section 4.
In Fig. 5, the broken-norm and the L 2 -norm errors, as well as T F h and T W h , are shown for the WACA method [10]. WACA exhibits the same accuracy of Mortar and Internodes methods when P 1 discretization is adopted in both subdomains, while it is sub-optimal for high-order spectral element discretizations.   from the side k of the interface is approximated accurately, as in the case of the WACA method.
Finally, in Fig. 9 we report the errors and the quantities |T W h |, |T F h | in the case of a decomposition that is conforming at the interface, i.e., with mesh-size h = h 1 = h 2 and polynomial degree p = p 1 = p 2 . These pictures highlight that T W h and T F h are not null also in the case of conforming decomposition, without however detracting from the conservation properties of the multidomain approach; more precisely, T W h and T F h are O(h p ) when h 1 , h 2 → 0 uniformly. In the conforming case, T F a,h = T W a,h = 0 in virtue of the interface conditions (16) 2,3 . The two meshes are regular and uniform with k × k elements in each subdomain, k = 10, 20, . . . , 50 for P 1 and k = 8, . . . , 20 for Q 3 .
For this test case we report numerical results only for the Internodes method (see Fig. 11), since, with the help of RL-RBF interpolation, Internodes does not feature any additional difficulty with respect to the case of geometric matching interfaces. On the contrary, the implementation of Mortar method for curved interfaces is far from trivial: it requires several steps such as projection, intersection, local meshing and numerical quadrature to build up

Test Case #3: d = 2, Discontinuous Coefficients and M > 2 Subdomains
Let us consider now the Kellogg's test case (see, e.g., [17,23]), in which the function ν is piece-wise constant and γ = 0. The exact solution can be written in terms of the polar coordinates r and θ as u(r, θ) = r α μ(θ), where α ∈ (0, 2) is a given parameter, while μ(θ) is a 2π-periodic continuous function (more regular only when α = 1); see Fig. 13. When   Fig. 13 Test case # 3. The Kellogg's test case. On the left, the decomposition of Ω into four subdomains. In the middle, the nonconforming P 1 meshes for k = 10. On the right, the Kellogg's solution with α = 0.4 and γ 1 = 9.472135954999585 computed by INTERNODES and P 1 α = 1, u ∈ H 1+α−ε (Ω), for any ε > 0; the solution features low regularity at the origin and its normal derivatives to the axis are discontinuous. The positive value γ 1 depends on α and on two other real parameters σ and ρ. The set {γ 1 , α, σ, ρ} must satisfy a nonlinear system (see formula (5.1) of [23]). In particular we fixed ρ = π/4. We consider here the domain splitting as well as the discretization described in [17], where the numerical convergence order of the Internodes method has been shown for different values of the parameter α. In [17], the convergence estimate provided by Theorem 1 for two subdomains has been confirmed, although this test case involves four subdomains instead of two.
Here we show that also the conservation properties of the Internodes method are preserved when the problem features discontinuous coefficients. In Fig. 14 we compare the total force and the total work (in both strong and weak form) with the H 1 -broken norm error and the L 2 -norm error, for P 1 (Q 2 , resp) discretization in each subdomain when α = 0.6 (α = 1.8, resp.). Recalling that u ∈ H 1+α−ε (Ω), it is useless to consider higher polynomial degrees.

Test Case #4: d = 2, M > 2 Subdomains and Conservation Properties vs. p
This is another test case with more than two subdomains in which we show both the convergence and the conservation properties of the Internodes method with respect to the local polynomial degree p.
The domain Ω = (0, 3) 2 is split in five subdomains as in a tatami (see Fig. 15). We have considered the same polynomial degree p in each subdomain, while the mesh sizes are all different, as we can evince from the left picture of Fig. 15. In the right picture of Fig. 15 we show the errors between the Internodes solution and the exact one u(x, y) = cos((x + y)π/2)(x − 2y), as well as the weak and strong total forces and total works. The coefficients of the problem are ν = 1 and γ = 0. Dirichlet boundary conditions are imposed on the whole ∂Ω.
The curves plotted in the right picture of Fig. 15 show that the errors (H 1 -broken norm and L 2 -norm), as well as the weak and strong total forces and total works decay more than algebraically with respect to p, as it is typical for spectral elements (or hp-fem) discretizations.   Fig. 15 Test case #4. On the left, the domain splitting (the different colors refer to the different subdomains, the quads are the elements of the mesh, inside each element there are (p +1) 2 nodes) and the non-conforming meshes. On the right, the errors, the total forces and the total works at the interface vs. the polynomial degree p
First, we have set the interface Γ α with α = 3 and we have considered Q 2 discretizations in both the subdomains with variable h 1 = 1/k and h 2 = 1/(k + 1), for k = 2, . . . , 9. RL-RBF interpolation has been used to build the interpolation matrices of Internodes, since the interface cannot be described exactly by Q 2 isoparametric elements, thus we have geometric non-conformity. The errors, as well as the total forces and the total works are shown in the left picture of Fig. 16: they decrease at least like h 2 , i.e., like the H 1 -broken-norm error.   Fig. 16 Test case #5. On the left, Q 2 discretization, cubic interface and RBF interpolation, geometric non-conforming. On the right, Q p discretization, quadratic interface and Lagrange interpolation, geometric conforming Then, we have set the interface Γ α with α = 2 and we have considered Q p Spectral Element discretization in both the subdomains with fixed h 1 = 1/4 and h 2 = 1/3, for p = 3, . . . , 7. Now, Lagrange interpolation has been considered to build the intergrid matrices of Internodes. We can use Lagrange interpolation since the interface can be described exactly by Q p spectral elements 4 with p ≥ 2, thus the interfaces are geometric conforming. The total forces and the total works are shown in the right picture of Fig. 16: they decrease versus p like both the H 1 -broken-norm and the L 2 -norm errors. The numerical results confirm the theoretical estimates of Section 4 also in this 3D test case. This work has been inspired by our long and fruitful friendly collaboration with Alfio Quarteroni, whom we wish to warmly thank. Indeed the contents of this paper have also been subject of common discussions with him.
Funding Open access funding provided by EPFL Lausanne.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.