An a priori error estimate for interior penalty discretizations of the Curl-Curl operator on non-conforming meshes

We prove an a-priori error estimate for regularized Curl-Curl Problems which are discretized by the Interior Penalty/Nitsche’s Method on meshes non-conforming across interfaces. It is shown that the total error can be bounded by the best approximation error which in turn depends on the concrete choice of the approximation space Vh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$V_{h}$\end{document}. In this work we show that if Vh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$V_{h}$\end{document} is the space of edge functions of the first kind of order k we can expect (suboptimal) convergence O(hk−1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$O(h^{k-1})$\end{document} as the mesh is refined. The numerical experiments in (Casagrande et al., SAM Report 2014-32, ETH Zürich, 2014) indicate that this bound is sharp for k=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$k=1$\end{document}. Moreover it is shown that the regularization term can be made arbitrarily small without affecting the error in the |⋅|curl\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\lvert \cdot \rvert_{\mathbf {curl}}$\end{document} semi-norm. A numerical experiment shows that the regularization parameter can be chosen in a wide range of values such that, at the same time, the discrete problem remains solvable and the error due to regularization is negligible compared to the discretization error.


Introduction
In this work we study the D, magnetostatic boundary value problem, which can be used to calculate the magnetic field that originates from a divergence free, stationary current j i .Herein μ denotes the magnetic permeability and g D prescribes Dirichlet boundary data.We seek the magnetic vector potential A that fulfills ()-().The magnetic field is then B = ∇ × A. Note that if g D ≡  on ∂ then () implies (∇ × A) • n = B • n =  on ∂ which reflects the decay of the fields away from the source.Note that the solution of the boundary value problem ()-() is only unique up to a gradient field (if is simply connected), which is not of importance if one is only interested in the magnetic field B. Thus it is possible to solve the ungauged problem ()-() if the current j i lies in the range of the system matrix [].The latter is hard to enforce on nonconforming meshes (cf.Section ) and it is simpler to gauge the formulation ()-() or to add a regularization term to () so that the system matrix has full rank.In case of adding a regularization term to (), one introduces a modeling error which must not dominate the approximation error of the numerical scheme.
In some applications like the simulation of electric machines or magnetic actuators, magnetic fields have to be computed in the presence of moving, rigid parts.Then one may use separate, moving sub-meshes for them in order to avoid remeshing.However, this leads to so-called 'sliding interfaces' , i.e. meshes with hanging nodes (cf. Figure ).
Our goal is to construct a method that approximates the solution of ()-() in a way that is independent of the 'non-conformity' of the sub-meshes at the common interface.This problem has been tackled successfully in the framework of Mortar methods where the continuity constraints are incorporated directly into the trial-space [] or they are enforced by additional Lagrange Multipliers [].However they come at the price of introducing either non-local shape functions or additional unknowns.
Another approach uses the Discontinuous Galerkin (DG) framework to solve problem ()-() in the presence of hanging nodes.In [] problem ()-() is regularized by adding a ∇(∇ • A) term to () and is then solved by the locally discontinuous Galerkin method.However because of the additional regularization term, additional assumptions on the smoothness of the solution have to be made to prove convergence.Alternatively one can use a mixed DG formulation and enforce the gauge condition ∇ • (μ - A) =  explicitly to avoid the introduction of a regularization term [].The stability of this method for arbitrary, sliding meshes remains unclear: In [] it is proven that the mixed method yields the expected rates of convergence on conforming meshes and the experimental results in [] show that it also works on adaptively refined meshes with hanging nodes.However, in light of the results in Section . and in [] it is not clear that the constant in the inf-sup condition of [] is independent of the 'non-conformity' of the sub-meshes at the common interface.
We study a different approach: We apply the Interior Penalty/Nitsche's Method [] to the regularized magnetostatic problem, Here  < ε  is the regularization parameter.In an earlier investigation [] it was shown experimentally that the Interior Penalty Method solves problem ()-() successfully if second order edge functions of the first kind are used.Moreover it was shown that first order edge functions fail to converge to the exact solution as the mesh is refined.In this work we intend to give theoretical explanations of these observations and investigate the error that is introduced by the regularization term in ().
We start our discussion in Section  by introducing Discontinuous Galerkin (DG) notations that were already introduced in [] and which are needed to state the interior penalty formulation of ()-() in Section .Section  also proves an a-priori bound of the total error in terms of the best approximation error for piecewise-polynomial test-and trial spaces V h .In Section  we analyze the particular case where V h is the space of k-th order edge functions, R k .Combining the results of Sections  and  we get rates of convergence for the regularized problem ()-().Section  is devoted to the choice of the local length scale appearing in the Interior Penalty formulation and presents numerical experiments underlining the results of Sections , .Section  discusses the role of the regularization parameter ε and how to choose it.We end our presentation with a short conclusion and outlook in Section .

Preliminaries
Before we can introduce the Symmetric Weighted Interior penalty (SWIP) formulation of ()-() we give some definitions and notations (cf.[]): Subdomains and submeshes: Let us assume that the domain , on which ()-() is posed, is a simply connected polyhedron with Lipschitz boundary.Furthermore we assume to be split into two non-overlapping subdomains,  ∪  = .
We introduce a sequence of simplical meshes T H = (T h ) h∈H on .Here H denotes a countable subset of R + having  as the only accumulation point.For each h ∈ H we let T h ∈ T H denote a particular mesh in the sequence T H and we let T ∈ T h be a mesh element (tetrahedron).The meshwidth is defined as h = max T∈T h h T , where h T is the diameter of element T.
We assume that each mesh T h , which covers , can be split into two conforming, nonoverlapping submeshes, T h = T h, ∪ T h, , that cover  and  , respectively.As before we define T H, = (T h, ) h∈H and T H, = (T h, ) h∈H .
Furthermore we define F T to be the set of the four facets of a tetrahedron T. The intersection of two facets, is called an inner face while the intersection of a facet with the boundary ∂ is called a boundary face.Note that facets are always triangular while inner faces are convex polygons with up to six nodes and boundary faces can have virtually any polygonal shape (cf. Figure ).We denote by F b h the set of all boundary faces, by F i h the set of all inner faces and define F h = F b h ∪ F i h to be the set of all faces.Furthermore, F T stands for the set of all faces which lie on the boundary of element T.

Mesh assumptions:
We assume that the elements are shape regular in the sense of Ciarlet: There is a constant σ max , independent of h, such that for all h ∈ H and for all T ∈ T h we have where ρ T is the radius of the largest ball inscribed in T. It is easy to check that this condition is satisfied if two sequences of static sub-meshes are moved against each other.We will make additional assumptions about the mesh when we discuss choices for the local length scale in Section .

Magnetic permeability:
We assume there exists a partition P = { i,μ } such that each i,μ is a polyhedron and such that the permeability μ >  is constant on each i,μ .Furthermore the mesh sequence T H is compatible with the partition P : For each T h ∈ T H , each element T ∈ T h belongs to exactly one i,μ ∈ P .I.e. the magnetic permeability is constant on each element but it is allowed to jump element boundaries, and in particular over the non-conforming interface :=  ∩  .
Polynomial approximation: Later on we will seek the discrete solution in the piecewise polynomial space (cf.[]) where T h ∈ T H and P k (T) is the usual space of polynomials up to total degree k on mesh element T. L  ( ) is the usual space of square integrable functions on .Note that functions of P k (T h ) are discontinuous across element boundaries.

Mesh faces, jump and average operators:
For each mesh face F and vector valued function A h ∈ P k (T h )  , we define the tangential jump as follows: for The weighted average is defined similarly: Here the normal n F always points from T  to T  and ω  , ω  ∈ [, ] such that ω  + ω  = .Note that the jump and average operators are well defined for all p ∈ P k (T h )  .
The following lemma relates the trace of a polynomial function to its L  norm on the element (cf.[], Lemma .): Lemma  (Discrete trace inequality) Let T H be a sequence of shape regular, possibly nonconforming, simplical meshes.Then for all h ∈ H, all v h ∈ P k (T h ), and all T ∈ T h we have Herein C tr is independent of T and h but depends on σ max , k.
Function spaces: We will use the spaces Herein H s ( ) = W s, ( ) is the Sobolev space of order s with Hölder coefficient p = .The associated norms and semi-norms are

Symmetric Weighted Interior penalty (SWIP) formulation
We chose an arbitrary subspace V h ⊆ P k (T h )  as discrete test and trial space, and use integration by parts (cf.[, ] for details) to arrive at the SWIP formulation of (): Find Here, where η is the penalty parameter.The last four terms of a SWIP h are called consistency, symmetry, penalty, regularization term, respectively.For an inner face F ∈ F i h , F = ∂T  ∩ ∂T  , we chose the weights If F is a boundary face, F ∈ F b h , we choose γ μ,F := μ - .The term a F is the local length scale of face F and can be chosen in different ways (e.g. a F =   (h T  + h T  ) where h T  , h T  are the diameters of the neighboring elements).For now we assume that there exists a constant ς  >  such that for all h ∈ H, all T ∈ T h , and all F ∈ F T : In Section  we will look at concrete choices of a F and discuss the circumstances under which () is fulfilled.It will turn out that depending on the choice of a F we have to make additional assumptions about the mesh regularity to guarantee ().

Remark  If V h ⊆ H(curl;
), then all inner tangential jumps in () will drop out [], Lemma ., and only jumps at the boundary remain, i.e. we are left with a standard FEM formulation where the inhomogeneous boundary conditions () are enforced in a weak sense.

A priori error estimate
In the following we derive an error estimate in the 'energy-norm' for the variational problem ().

Regularity of the exact solution:
We assume that the exact solution A of ()-() (in the sense of distributions) is such that Furthermore, we set V * h := V * + V h .Note that, because A and ∇ × A are in H  (P )  the traces of A and ∇ × A are well defined on the faces of the mesh elements (cf.[], Remark .).Indeed, let T ∈ T h be a mesh element, then by the multiplicative trace inequality H  (T)  , and the same estimate holds for ∇ × A. Therefore we see that a SWIP h : V * h × V h → R is well defined.In order for the right-hand side to be well-posed we assume We begin the proof of the a priori error estimate by showing that the exact solution A fulfills equation (): ).Thus all inner jump terms drop out, Note that the last two sums include only boundary faces.Next we make use of the following identity (which holds for any interior face Here Let us apply the identity to the second term of (): - The second term on the right-hand side vanishes because A is a solution of the strong formulation ().Thus • n F so we can rearrange the face contributions to the element boundaries, - Now substitute () into () and use integration by parts on each mesh element [], Theorem .: Let us introduce the following (semi-)norms on the space V * h : Lemma  (Bound on consistency term) For all A, A ∈ V * h there holds Here h T := max{ x -y |x, y ∈ T} is the diameter of mesh element T ∈ T h .
Proof For an arbitrary inner face F = ∂T  ∩ ∂T  we have by the Cauchy-Schwarz (CS) inequality By using Cauchy-Schwarz again we see that Substitute this back into () to get Similarly, for a boundary face Now use ()-() to bound the sum over all faces, where we have regrouped the face contributions and used that a F ≤ ς  h T in the last step, cf.().
Using Lemma  we can finally prove discrete coercivity.

Lemma  (Discrete coercivity)
The bilinear form a SWIP h is coercive: For all η > C  tr ς  and all h ∈ H there holds The constant C tr stems from the discrete trace inequality () and is independent of h, μ, ε, ς  .
Proof By definition of a SWIP h we have Now let us give a bound on the second term on the right-hand side using Lemma , where we have used the discrete trace inequality () componentwise in the last step.Hence, Now use the inequality x  -βxy + ηy  ≥ η-β  +η (x  + y  ) which holds for arbitrary β, η, x, y as outlined above (it follows from (βxηy)  + (xβy)  ≥ ): Finally, we note that C stab >  if η > C  tr ς  which completes the proof.

Lemma  (Boundedness)
There exists a constant C bnd >  independent of h, μ, and ε such that for all Proof We start by splitting the bilinear form a SWIP h into five terms, We can now bound these terms individually, Finally, we can combine the previous results into one theorem.
Theorem  (Error estimate) Let A ∈ V * be a solution of the strong formulation ()-() (in the sense of distributions) and let Then there exist constants C > , C η > , both independent of h, μ, and and the discrete problem () is well-posed.The constant C η depends on ς  , k and C depends on η, ς  , k.
This theorem tells us that the total error is bounded by the best approximation error (w.r.t.suitable norms).Note that we didn't make any assumption on how the submeshes T h, and T h, meet at .In order to get rates of convergence we will have to make additional assumptions about the approximation space V h and the exact solution A. This will be the topic of Section .
Proof of Theorem  In this proof C denotes an arbitrary, positive constant that is independent of h, μ, ε, and that may take a different value every time it used.We begin by picking an arbitrary v h ∈ V h .Then, by the triangle inequality, Inserting this bound into () (which holds for arbitrary v h ) yields the assertion.Note that the bilinear form a SWIP h is coercive (Lemma ) and bounded (finite dimensions).Thus, Lax-Milgram assures the well-posedness of the discrete problem.
Remark  Observe that for ε →  the variational formulation () becomes ill-posed.To see this we observe that the • SWIP norm 'becomes' a semi-norm as ε → .In order to study the behavior as ε →  it is thus desirable to state the discrete coercivity (Lemma ) w.r.t. a norm that does not depend on ε: We use that A h  SWIP ≥ ε μ -/ A h  L  and thus Lemma () can be rewritten as (   ) We see now clearly that the coercivity constant depends linearly on ε, i.e. the discrete problem becomes ill-posed as ε → .

Rate of convergence for edge functions
In the following we will bound the best approximation error appearing in Theorem  for edge functions of the first kind.For this we assume, in addition to (), that a F is uniformly bounded from below in the sense that there exists a constant ς  such that for all h ∈ H, all T ∈ T h and all F ∈ F T we have For the remainder of this section, let us choose where R k is the space of k-th order edge functions (k =  are the lowest order, H(curl) conforming Whitney elements, cf.[], Eq. (.)).Because the sub-meshes T h, , T h, are conforming, the spaces R k (  ), R k (  ) are H(curl) conforming.We can thus use the standard projection operator r h as it is defined in [], Section ., for edge functions on  ,  to build our global projection operator π h : The following theorem then gives an upper bound for the best approximation error of Theorem : Here C depends on μ, ς  , k, ε but not on h.Moreover if ε < , C is independent of ε.
Remark  By combining Theorem  with Theorem  we see that for a sufficiently smooth exact solution A, the total error A -A h SWIP = O(h k- ) if k-th order edge functions are used.In comparison to standard FEM on conforming meshes one order of convergence is lost.Theoretically it is possible that there exists another projector πh which would give a better rate of convergence, but numerical experiments show that Theorem  is sharp for k =  (see Section ).
In order to prove the above theorem we will make use of two Lemmas to bound the face contributions.
Lemma  Let T H, be a sequence of shape regular, conforming, simplical meshes of the domain  .Suppose there exists an integer where C is independent of h T , T.
For the proof of Lemma  we refer the reader to [], Lemma . (which is proven element-wise).
Lemma  Let T H, be a sequence of shape regular, conforming, simplical meshes of  .Assume u ∈ H s (  )  for some integer  ≤ s ≤ k and u transforms such that it preserves the divergence, i.e. if F : T → T, û → u is an arbitrary mapping then u transforms as Then the following estimate holds: where w T : H  (T)  → D k is the standard (local) interpolation operator for k-th order Thomas-Raviart elements D k [], Section ..The constant C does not depend on h T , T.
Proof In order to simplify notation we will assume in this proof that C >  is an arbitrary constant independent of h, T that may take a different value every time it is used.We note that since u ∈ H s (T)  , w T u is well defined by [], Lemma ..Now split the integral over ∂T into its facet contributions, Since our mesh contains only tetrahedrons we can find for every F T ∈ F T a linear transformation T,F T : T → T which maps the reference element T onto the actual element T such that the pre-image FT of facet F T lies in the x-y plane of T, where B T,F T ∈ R × .Now using the usual change of variables together with () we obtain Here (B T,F T ) :,i denotes the i-th column of B T,F T , w T u is defined by (), and we have used that w T u = w T û [], Lemma ..Now notice that û -w T û ∈ H s ( T)  and thus we can use the trace inequality [], Theorem ., For the next step we note that ∀φ ∈ P k- ( T)  ⊂ D k ( T) we have φ = w T φ by the definition of w T (D k is the H(div; ) conforming Raviart-Thomas space, see [], Section .).Therefore, where we have used that w T : Since φ is arbitrary we can use the Deny-Lions theorem [], Theorem ., Finally we have to map | û| H s ( T)  back to the actual element T. For this observe that using (), with |α|  = s being a multi-index.Therefore, where we have used [], Lemma ., in the last step.Now combining ()-() gives.
Here we have used [], Lemma ., together with the fact that the mesh sequence is shape regular.Now summing over all facets F T ∈ F T yields the assertion.
Using these Lemmas we can finally give a bound for Aπ h A SWIP, * .
Proof of Theorem  In order to simplify notation, C denotes in this proof an arbitrary, positive constant that is independent of h.Note that the interpolation operator r h (A|  ) is well defined for s ≥  by the Sobolev Embedding Theorem and [], Lemma ..Because the sub-meshes of  ,  are conforming themselves, r h (A|  ) is tangentially continuous across all inner, conforming faces.The same holds for  and because A ∈ H(curl; ) the exact solution is also tangentially continuous across all inner faces.Therefore only jump terms across the faces

we have to bound
Since μ is piecewise constant on each i,μ ∈ P there are constants μ min , μ max such that  < μ min ≤ μ ≤ μ max .T  and T  are easily bounded using [], Theorem .: The term T  is bounded using Lemma , Here we have used that a F ≥ ς  h T .
In order to bound the term T  we first note that the global Thomas-Raviart interpolation operator w h (∇ × A| i ) i=, is well defined by [], Lemma ..Thus, ∇ × [r h (A| i )] = w h (∇ × A| i ) by [], Lemma ., and we can bound T  as follows: where we have used Lemma  and the fact that h T ≤ h.
Remark  From the proof of Theorem  it is clear that for h sufficiently small the term T  dominates the other three terms and is thus responsible for the loss of one order of convergence as pointed out in Remark .Interestingly T  sums the jump terms only over the faces F b, h .This suggests that it suffices to use (k + )-th order edge functions in elements adjacent to F b, and k-th order edge functions everywhere else to achieve O(h k ) order convergence.This can be implemented easily by using a hierarchical basis for the edge functions [].

The local length scale a F and h-convergence
So far we have assumed that the local length scale a F fulfills (), () in order to derive an a-priori error estimator, i.e.  < ς  h T ≤ a F ≤ ς  h T .We will now study the following three concrete choices for a F : • a () , where h T  , h T  are the diameters of the adjacent elements of face F and h F is the diameter of face F. It turns out that for each choice of a F we have to make additional assumptions on the mesh such that a F fulfills (), ().So once we have chosen a concrete a F we can think of ς  , ς  as mesh dependent parameters.The important point is that the constants C in Theorems  and  depend on the constants σ max , ς  , ς  but they do not depend in any other way on the shape of the underlying meshes.Hence, if we can show that σ max , ς  , ς  are independent of the way that T H, , T H, intersect at the sliding interface , then there must be an upper bound on the total error A -A h that is independent of the relative position of T H, to T H, and that tends to  as h → .
Let us now discuss the precise conditions on the mesh for each choice of a F : For a () F , a () F we require T H to be quasi-uniform at the sliding interface : Definition  A mesh-sequence T H is said to be quasi-uniform at if it is shape regular () and if there exists a constant σ  >  such that for all h ∈ H, all T ∈ T h := {T ∈ T h | ∂T ∩ = ∅}: Lemma  If the mesh is quasi-uniform at then a () F , a () F fulfill conditions (), () and the constants σ max , ς  , ς  are independent of the way T H, , T H, intersect at .

Proof a ()
F ≥   h T i follows immediately from the definition for i = , .For the other direction we use () and get a () The lemma above asserts that the choices a () F , a () F lead to a method that converges independently of the way that the two mesh-sequences T H, , T H, intersect at .In particular the faces can be very tiny 'slivers' (i.e.triangles with high aspect ratio).But note that the choice of a F determines the required minimum value of the penalty parameter (see Lemma ).
By substituting a () F into () we see that we need an estimate of the form h F ≥ ς  h T in order for Theorem  to hold.However if two meshes are sliding against each other such an estimate is not feasible since h F can become arbitrarily small in comparison to h T .In other words, the constant ς  depends on the way T H, intersects with T H, .Nevertheless using a ()  F in the variational formulation  seems to work in practice (see below).

Numerical examples
We study the behavior of the SWIP formulation for the three different choices of the local length scale a F ; As in [] we consider a D sphere with radius  that is split into two hemispheres,  and  (cf. Figure ).For each hemisphere we create a sequence of quasiuniform meshes, T H, and T H, , which approximate the boundary linearly.We impose the analytical solution A = (sin y, cos z, sin x) and choose j i , g D such that they fulfill ()-().
Figure  shows the error in the curl-semi-norm for different angles of rotation, for all three choices of a F , and for different mesh-sizes h.We can see that although the error depends slightly on the angle, it converges to zero in all three formulations as h is decreasing.Moreover we see that the choices a () F , a () F yield similar results which are slightly better than the choice a () F .In order to illustrate the estimates of Theorems  and  we plot the error for a series of quasi-uniform meshes in Figure  for a F = a () F . a As in [] no convergence is observed when first order edge functions (k = ) are employed which implies that Theorem  is sharp for k = .For k =  and k =  we observe rates of convergence O(h . ) and O(h . ), respectively, which affirms Theorems  and .
For V h = P  (T h )  we observe the rate of convergence O(h), i.e. there is no loss of one order of accuracy.This is because P  (T h )  spans the full polynomial space (see [] for a proof ).used for discretization and ε = 10 -6 , η = 50, μ ≡ 1.Note that the curve for a (2)   F is partially hidden by the one of a (1)  F .

Remark 
The observed rates for k =  and k =  are higher than the rates O(h), O(h  ) which we expect from Theorem .This is due to the better approximation properties of the edge functions in the inside of the two hemispheres (cf.Remark ).
Remark  Strictly speaking this numerical experiment does not fit the framework developed in Sections ,  because is not a polyhedron.Extending the theory to domains with curved boundary ∂ is beyond the scope of this paper but Figure  suggests that the order of convergence for k = , k =  is the same as for polyhedral domains.

The regularization parameter ε
So far we have looked at the regularized system ()-() and it was shown that the proposed method yields the expected rates of convergence for ε > .However, genuine magnetostatics amounts to choosing ε = .We will consider two approaches to solve the system ()-() with ε = : On the one hand we will try to set ε =  directly and on the other hand we will study the effect of choosing ε small enough such that the error due to regularization is negligible.

The case ε = 0
If we set ε = , the boundary value problem ()-() does not have a unique solution.Indeed the continuous curl-curl operator has an infinite-dimensional kernel and the non-zero eigenvalues are well separated from  [], Corollary ..If one uses H(curl) conforming edge functions of the first kind on a conforming mesh it can be shown that the discrete curl-curl operator has a (finite-dimensional) kernel and that the discrete eigenvalues are well separated from it [], Discrete Friedrichs inequality, Lemma ., i.e. edge functions of the first kind yield a spectrally accurate discretization of the curl-curl operator.From a theoretical point of view it remains unclear whether this property carries over to the SWIP formulation (), cf.
[].Therefore the spectrum of the a SWIP h bilinear form is investigated in a numerical experiment.The setup is very similar to the one in the previous section: The domain consists of two half-spheres which can be rotated against each other by an angle θ .However this time we only assemble the matrix of the a SWIP h bilinear form with ε = , a F = a () F b and compute its eigenvalues using the eig routine of MATLAB Ra.
Figure  shows the smallest and largest non-zero eigenvalues of the SWIP formulation for different mesh-widths h and different angles θ (dashed, blue lines) (an eigenvalue has been classified as non-zero if its absolute value is greater than  - ).For comparison we have also plotted the eigenvalues of a standard H(curl) conforming discretization using second order edge functions on the conforming grid with θ =  (green lines).
We see that the bandwidth of the SWIP eigenvalues is comparable to the bandwidth of the H(curl) conforming discretization for many angles.But we also observe that for some angles the lower end of the spectrum tends to zero.In order to better understand this phenomena we plotted the smallest/largest non-zero eigenvalues of the SWIP discretization against θ for one mesh-size (Figure ).We now see that the lower end of the spectrum deteriorates as θ → , i.e. we can expect spectral pollution for very small angles.This agrees with the observations of [].
The previous considerations indicate that the a SWIP h bilinear form is not suitable to solve the Maxwell Eigenvalue Problem.However in this work we are concerned with the curlcurl source problem ()-().Although the Galerkin matrix becomes singular for ε =  we can in principle still solve the linear system if it is consistent, i.e. if the right-hand side lies in the range of the Galerkin matrix.Then the solution A h is not unique anymore, but curlA h is.
We attempt to solve the linear system of equations using the conjugate gradient (CG) method [].In [] it is shown that the CG method converges for consistent, symmetric   This has been confirmed in a numerical experiment: We take the example from Section  with the same analytical solution and chose the right-hand side  provides the number of CG iterations required to reach the prescribed tolerance  - .We see that without a preconditioner the computational cost for the angle θ =  - is almost  times larger than for θ =  - .For comparison we also list the number of iterations needed when the multi-level ILU decomposition ILUPACK is employed c [].In this case the number of iterations also increases but the factor  is reduced to ≈ ..
Remark  Although the right-hand side j i chosen in the numerical experiment above is clearly divergence free, there is no guarantee that its discrete counterpart h is so too, i.e. it is not clear that the right-hand side vector b, that is associated with h , lies in the range of the system matrix.We have investigated this by splitting the right-hand side vector b into a part that lies in the kernel of the system matrix, b, and into it's orthogonal complement, b⊥ .It turned out that for all angles b⊥  / b  ≈  - , which seems to be sufficient for CG to converge.
We can conclude that setting ε =  is in principal possible if the right hand side vector b lies in the range of the system matrix.However checking this for non-zero right hand sides j i is a non-trivial task because we don't know a-priori the kernel of the system matrix a SWIP h .Moreover the system matrix becomes ill-conditioned as the angle θ →  which causes an increase in the number of CG iterations.
Remark  For H(curl) conforming discretizations, which fulfill the discrete sequence property, the kernel of the system matrix is known.In particular it is easily proven that div j i =  implies that h lies in the range of the system matrix.Unfortunately it is not clear whether this property carries over to the SWIP formulation () because to the best of our knowledge there exists no characterization of the kernel of a SWIP h on arbitrarily nonconforming meshes.

The case 0 < ε 1
We saw in the previous section that setting ε =  is in practice not feasible.Therefore we study a different approach: We choose ε so small that the error due to regularization becomes negligible.To make this more explicit we bound the total error between the discrete, regularized solution A ε h and the exact solution of ()-() with ε = , A  , by two contributions: herein A ε is the exact solution of the regularized system ()-().Clearly the second component is independent of the discretization and thus h, but it depends on ε for a given problem.Moreover, the first term depends on h but is independent of ε because the constant C of Theorems  and  is independent of ε.
It is thus desirable to choose ε small such that ∇ × However, as ε →  the discrete problem becomes ill-posed and solvers typically fail to converge, cf.Remark , Section ..
We try to circumvent this problem by two approaches: • For small problems we use the Sparse Cholesky Decomposition of PARDISO [] (Intel MKL Version .) and solve the linear system of equations directly.• For problems whose Cholesky Decomposition does not fit into memory we use the Conjugate Gradient Method together with ILUPACK [] as a preconditioner (using the settings of Section .).
Remark  We are only interested in the curl of the solution, i.e. the magnetic field B. If we were to look at A instead of ∇ × A then A ε h -A ε L  ( )  would not be independent of ε as can be seen from Theorem .
The following Lemma gives us a guideline for choosing ε.
Lemma  If we impose homogeneous Dirichlet data, g D ≡ , we have, where α is the smallest, non-zero eigenvalue of ∇ × (μ - ∇ × A) = αA.If μ is constant in we can choose α = √ π/l max where l max is the maximum side length of a cuboid that contains .
The first part of this lemma is proven in [], Lemma ., and for the second part we use a result of [], Section .
By comparing () with () we see that μ ) for k order edge functions.Therefore ε should be chosen such that ε ∼ h k- as the mesh is refined.

Numerical example:
We consider the same setup as in the previous section (cf. Figure ) with but we choose a different μ for the upper and lower hemisphere.The analytic solution is chosen as A  = (sin y, , μ sin x).j i is chosen such that A  fulfills ()-() with ε = .We solve the system of linear equations using PARDISO for different values of ε and μ (as in the previous section we choose a F = a () F ). Figure  shows the total relative error μ -/ ∇ × (A ε h -A  ) L  ( )  as a function of ε for various mesh-sizes.The solid lines show the error for μ upper /μ lower =   whereas the dashed lines show it for μ upper /μ lower =   .
We note that the errors are almost identical for both choices of μ.Moreover, we observe that for ε <  - the discretization error μ -/ ∇ ×(A ε h -A ε ) L  ( ) (which here includes the error due to boundary approximation, cf.Remark ) clearly dominates the regularization error μ -/ ∇ × (A ε -A  ) L  ( ) whereas for ε >  - the regularization error is dominated by the discretization error.This is what we can expect from the previous discussion.In fact, if we use Lemma  and fit the two hemispheres into a cube of side length l max = , Remark  The same results are obtained if CG together with ILUPACK is used.For brevity we omit these results here.
We would like to point out that by using the direct solver PARDISO we were able to solve the resulting system of linear equations for ε as small as  - and that the time needed to solve the problem seems to be independent of ε (see Table ).A similar result holds for preconditioned CG with ILUPACK preconditioner where the system is solvable for arbitrary small ε (cf.Section .) and the solution time seems to be independent of ε for ε small enough.
We can thus choose ε (almost) arbitrarily small without affecting the discretization error μ -/ ∇ × (A ε h -A ε ) L  ( ) and incurring rising cost for solving the resulting linear systems of equations.In other words, one should choose ε as small as possible such that the resulting linear system can still be solved.

Conclusion and outlook
We have proved a-priori error estimators for the interior penalty formulation of the regularized curl-curl source problem ()-(); If the solution is approximated by k-th order edge functions we can expect at least convergence of order O(h k- ) (provided the exact solution is sufficiently smooth).In particular, for k =  no convergence was observed in a numerical experiment [], which implies that our result is sharp.The reason for this is that R k does not span the full polynomial space P k .The bounds require the mesh to be quasi-uniform at the sliding interface but do not make any assumptions on how the sub-meshes abut at the sliding interface nor does the error estimate depend on it.This is confirmed by the numerical experiments and it is observed that the approximation is stable independent of the way the sub-meshes intersect.
Moreover the role of the regularization parameter ε has been investigated; For practical purposes one can choose ε (almost) arbitrarily small and solve the discrete problem with a direct solver or by using the preconditioned conjugate gradient method.The error due to regularization is then dominated by the discretization error of the regularized problem and is negligible.
Outlook: The proof of Theorem  suggests that it suffices to use nd order edge functions solely in elements adjacent to the non-conforming interface, respectively boundary faces, to achieve O(h) convergence.This would reduce the required number of unknowns drastically and should be pursued for practical applications.

Figure 1
Figure 1 Sliding meshes.Initially conforming sub-meshes become non-conforming when the upper sub-mesh starts moving.

Figure 2
Figure 2 Non-conforming overlap of two submeshes.Consider two cuboids stacked upon each other (similar to Figure 1) which are meshed individually.This figure shows the restriction of the two submeshes T h,1 and T h,2 to the common interface and the inner/boundary intersections are marked.

Figure 3
Figure 3The meshes for the two half spheres.The upper hemisphere is turned against the lower hemisphere by θ = 2.86 degrees to create a non-conforming mesh.

Figure 4
Figure 4 The relative H(curl) error vs. the rotation angle for three different choices of a F and 4 different mesh-sizes: h = 0.638174, 0.482025, 0.359644, 0.261798.Second order edge functions R 2 were

Figure 6
Figure 6 The smallest/largest non-zero eigenvalue for ε = 0 is plotted against the meshwidth for 50 different angles of rotation (dashed lines).For comparison the smallest/largest non-zero eigenvalue of a H(curl) conforming discretization based on second order edge functions is plotted as well.The angles are θ = 0.01n rad, n ∈ 0, . . ., 49 and R 2 edge functions were used to discretize a SWIP h , μ ≡ 1.

Figure 8 Table 2 a
Figure 8 Relative L 2 -error of curl vs. ε for multiple mesh-sizes h.The solid lines show the error for μ upper = 10, μ lower = 0.1 whereas the dashed lines show it for μ upper = 10 6 , μ lower = 0.1.The meshes have been rotated against each other by θ = 0.057 degrees and second order edge functions (k = 2) were used for discretization.