Primal and mixed finite element formulations for the relaxed micromorphic model

The classical Cauchy continuum theory is suitable to model highly homogeneous materials. However, many materials, such as porous media or metamaterials, exhibit a pronounced microstructure. As a result, the classical continuum theory cannot capture their mechanical behaviour without fully resolving the underlying microstructure. In terms of finite element computations, this can be done by modelling the entire body, including every interior cell. The relaxed micromorphic continuum offers an alternative method by instead enriching the kinematics of the mathematical model. The theory introduces a microdistortion field, encompassing nine extra degrees of freedom for each material point. The corresponding elastic energy functional contains the gradient of the displacement field, the microdistortion field and its Curl (the micro-dislocation). Therefore, the natural spaces of the fields are $[\mathit{H}^1]^3$ for the displacement and $[\mathit{H}(\mathrm{curl})]^3$ for the microdistortion, leading to unusual finite element formulations. In this work we describe the construction of appropriate finite elements using N\'ed\'elec and Raviart-Thomas subspaces, encompassing solutions to the orientation problem and the discrete consistent coupling condition. Further, we explore the numerical behaviour of the relaxed micromorphic model for both a primal and a mixed formulation. The focus of our benchmarks lies in the influence of the characteristic length $L_\mathrm{c}$ and the correlation to the classical Cauchy continuum theory.


Introduction
A common problem in the computation of materials with a pronounced micro-structure is the internal complexity of the geometry. In order to fully capture the kinematics one might resolve the underlying micro-structure. This can be done either with multi-scale finite element methods [21] or by modelling the finite element mesh to fully incorporate the microstructure. In both cases, the computational cost increases, leading to longer computation times and decreasing applicability of such approaches. Alternatively, one may enrich the mathematical model in order to account for the increase in the kinematical complexity. This approach gives rise to generalized continuum theories such as higher order gradient methods [29,36,48] or micromorphic continua [28,44,67]. Micromorphic continuum theories extend the kinematics of the material point with additional degrees of freedom, the choice of which defines the specific theory. Common examples are micropolar Cosserat [27,39], microstretch [60], and microstrain models [23,26]. The latter represent sub-types of the micromorphic continuum. In its most general case, micromorphic continua assume an affine deformable microbody for each material point. As such, this deformation, called here the microdistortion P , is fully captured by three-by-three matrices and introduces nine extra degrees of freedom. The full micromorphic continuum, as introduced by Eringen and Mindlin [22,37], incorporates the gradient of the microdistortion DP into the free energy functional. The resulting hyperstress term is a third order tensor. As such, it is unclear how this term is to be interpreted or applied. Typically, micromorphic continuum theories introduce a characteristic length scale parameter L c , which abstractly relates the dimension of the micro-body to that of the macro-body. In the full micromorphic model, when the characteristic length L c becomes very large, the microdistortion P must become constant in order for the theory to generate finite energies, which may lead to boundary layer problems [66].
The relaxed micromorphic continuum theory [46] takes a different approach by instead incorporating only the Curl of the microdistortion Curl P into the free energy function. The latter term, known as the micro-dislocation, remains a second order tensor and as such, induces a matrix-valued right-hand-side, known as the micro-moment M . Further, large characteristic lengths L c maintain finite energies [47,66]. The theory aims to capture the mechanical behaviour of both highly homogeneous materials and materials with a pronounced microstructure by governing the relation to the classical Cauchy continuum using the characteristic length L c [6] and shows great promise with respect to applications utilizing metamaterials, such as band-gap materials [7,16,34,35] and shielding against elastic waves [54,59]. Furthermore, analytical solutions have been derived for bending [56], torsion [55], shear [57], and extension [58] kinematics. The inclusion of the Curl of microdistortion in the free energy function implies the existence of unique solutions [24,45] in the space X = [H 1 ] 3 ×[H (curl)] 3 , as shown in Section 3.1.1. While the construction of finite elements for the space H 1 is well-known in the field of mechanics, H (curl)-finite elements are commonly used for the Maxwell equations, for example in magnetostatics [64]. For the construction of finite elements for H (curl) one may use Nédélec subspaces [40,41]. The formulation of higher order elements is detailed in [68]. For large characteristic length values L c the computation with the primal formulation (X h ⊂ X = [H 1 ] 3 × [H (curl)] 3 ) becomes unstable [66]. However, it can be re-stabilised by using a mixed formulation. The latter requires the employment of Raviart-Thomas- [53] or Brezzi-Douglas-Marini elements [13] and fully discontinuous finite elements as per the de Rham diagram [5,19,20].
In this work we demonstrate the existence and uniqueness of the primal formulation using the Lax-Milgram theorem, Section 3.1.1. Further, we introduce a mixed formulation, which is stable for large characteristic length L c values, Section 3.1.3. In Section 3.2 we derive the corresponding convergence rates for the discrete spaces. The construction of lower order finite elements for both formulations is explained in Section 4, with focus on a solution to the orientation problem and application of the discrete consistent coupling condition. Section 5 is devoted to numerical benchmarks of the finite element formulations but also features the convergence characteristics of higher order elements using NETGEN/NGSolve [61,63]. Finally, we present our conclusions and outlook in Section 6.

The linear relaxed micromorphic continuum
The linear relaxed micromorphic continuum [43,45,46] is described by its free energy functional, incorporating the gradient of the displacement field, the microdistortion and its Curl I(u, P ) = 1 2 Ω C e sym(Du − P ), sym(Du − P ) + C micro sym P , sym P (2.1) + C c skew(Du − P ), skew(Du − P ) + µ macro L 2 c CurlP 2 dX − Ω u, f + P , M dX , with u : Ω ⊂ R 3 → R 3 and P : Ω ⊂ R 3 → R 3×3 representing the displacement and the non-symmetric microdistortion, respectively. Here, C e and C micro are standard fourth order elasticity tensors and C c is a positive semi-definite coupling tensor for (infinitesimal) rotations. The macroscopic shear modulus is denoted by µ macro and the parameter L c ≥ 0 represents the characteristic length scale motivated by the geometry of the microstructure. The body forces and micro-moments are denoted with f and M , respectively. The differential operators are defined as Du =   u 1,1 u 1,2 u 1,3 u 2,1 u 2,2 u 2,3 u 3,1 u 3,2 u 3,3   , CurlP =    (curl P 11 P 12 P 13 ) T (curl P 21 P 22 P 23 ) T (curl P 31 P 32 For isotropic materials the material tensors have the following structure C e = 2µ e S +λ e 1 ⊗ 1 , C micro = 2µ micro S +λ micro 1 ⊗ 1 , where S : R 3×3 → Sym(3) and A : R 3×3 → so(3) are the fourth order symmetry and anti-symmetry tensors, respectively. Taking variations with respect to the displacement u δ u I = Ω C e sym Dδu, sym(Du − P ) + C c skew Dδu, skew(Du − P ) − δu, f dX , (2.4) and the microdistortion P δ P I = Ω − C e sym δ P , sym(Du − P ) − C c skew δ P , skew(Du − P ) + C micro sym δ P , sym P + µ macro L 2 c Curl δ P , Curl P − δ P , M dX , (2.5) Figure 1: Depiction of the reference domain Ω with an underlying micro-structure, the Dirichlet boundary Γ D according to the consistent coupling condition and the Neumann boundaries Γ u N and Γ P N for the displacement and the microdistortion, respectively. The vector ν is the unit normal vector on the boundary. The unit tangent vectors are illustrated with τ . The body force is given by f and the micro-moment by M . yields the symmetric bilinear form and the linear form for the load The corresponding strong form follows from partial integration, see Appendix A where ν denotes the outer unit normal vector, such that P × ν is the projection to the tangent surface on the boundary. The terms u and P are the prescribed displacement and microdistortion fields on Γ u D and Γ P D , respectively, see Fig. 1.
From a physical point of view, it is impossible to control the micro-movements of the material point on the Dirichlet boundary without also controlling the displacement. Consequently, the relaxed micromorphic theory introduces the so called consistent coupling condition [17] where the prescribed displacement u on the boundary automatically prescribes the tangential component of the microdistortion P on the same boundary, effectively inducing the definition Γ D = Γ P D = Γ u D .
3 Solvability and limit problems

Continuous case
In this section existence and uniqueness of the weak formulation of the relaxed micromorphic continuum for both a primal and a mixed method is discussed. In order to simplify the proof, homogeneous boundary conditions are assumed on the entire boundary, Γ u D = ∂Ω. The proof can be easily adjusted for mixed or inhomogeneous boundary conditions as long as the Dirichlet boundary does not vanish |Γ u D | > 0.
For the following formulations, we define corresponding Hilbert spaces and their particular norms from which we derive the corresponding spaces for our higher dimensional problem where both spaces are to be understood as row-wise matrices of the vectorial spaces.
Remark 3.1. In the following sections we assume a contractible domain Ω for all proofs.

Primal form
The proof of existence and uniqueness has already been given in [46] and subsequently generalized to the dynamic setting in [24]. For the reader's convenience and for later use we present a proof based on the Lax-Milgram theorem and Korn's inequalities for incompatible fields [8, 9, 31-33, 49, 50], similar to [45], where only the case µ c = 0 together with Γ P D = ∂Ω has been considered. Therefore, we define the product space X = [H 1 0 (Ω)] 3 × H (Curl, Ω) with its standard product norm {u, P } 2 X = u 2 H 1 + P 2 H (Curl) and examine the problem: Find {u, P } ∈ X such that where a(·, ·) and l(·) are defined as in (2.6) and (2.7), respectively.
, and C e , C micro be positive definite on Sym(3). Further, assume that C c is positive definite on so(3), or positive semi-definite and Γ P D = ∅. Then Problem 3.3 is uniquely solvable and there holds the stability estimate Proof. We show continuity and coercivity of (2.6). During the proof we denote with c > 0 a generic constant which may change from line to line. Continuity follows immediately with Cauchy-Schwarz inequality For the coercivity we first consider the case of positive definiteness of the tensors First, we consider the symmetric terms. Analogously to [66], with Young's inequality 1 and Korn's inequality 2 [42] we obtain We can choose ε = 1 2 (1 + c1 c1+cmicro ) such that both terms are positive. Next, we estimate the skew-symmetric part With 0 < δ < 1 only the second term is positive. By combining both estimates we conclude by choosing Thus, with Poincarè-Friedrich's inequality 3 we obtain the coercivity In the case of a positive semi-definite C c (even C c = 0 is allowed: absence of rotational coupling) together with Γ P D = ∅ we must use the generalized Korn's inequality for incompatible fields, cf. [9,9,24,32,33,49,50] sym P 2 and estimate Continuity of the right-hand side l is obvious and thus we can apply the Lax-Milgram theorem finishing the proof.
Remark 3.2. The proof above fails if we consider P ∈ [H 1 (Ω)] 3×3 instead of H (Curl, Ω). In fact, a(·, ·) is then no longer coercive as we cannot bound P 2 L 2 + Curl P 2 L 2 ≥ c P 2 H 1 uniformly in P . In [66] and [65] it is demonstrated that already in two spatial dimensions a loss of optimal convergence rates is obtained when using P ∈ [H 1 (Ω)] 2×2 . Remark 3.3. Note that in [51,52] it is shown that for certain sufficiently smooth data, the regularity of P can yet be improved to P ∈ [H 1 (Ω)] 3×3 .
3.1.2 Limit of vanishing characteristic length L c → 0 In the limit L c → 0 the bilinear form (2.6) reduces to (3.10) Therefore, we lose control over the Curl of P yielding a loss of regularity for P from H (Curl, Ω) to [L 2 (Ω)] 3×3 . We emphasise that the proof of Theorem 3.1 can be directly applied with the adapted product space together with the requirement of positive definite C c on the set of skew-symmetric matrices. Otherwise, control over the skew-symmetric part of P is completely lost in the limit leading to unstable results. Note that in this case no prescription of boundary conditions for P is possible.
Eq. (2.8b) can be reformulated as and used to express P algebraically P = (C e + C micro ) −1 C e sym Du + skew Du +M . Applying the latter to Eq. (2.8a) yields − Div(C macro sym Du) = f , The upper definition is derived in [6] and relates the meso-and micro-elasticity tensors to the classical macroelasticity tensor C macro of the Cauchy continuum, allowing to extract the macro material constants in the isotropic case .
In fact, C macro contains the material constants that arise from classical homogenization for large periodic structures. From Eq. (3.12) we observe that in general P is not a gradient field, even if M = 0. Therefore, by setting λ e = λ micro = 0 we obtain that skew P = skew Du and sym P = µe µe+µmicro sym Du. As µ micro > 0 we deduce that P is not a gradient field. This is a significant deviation from the two dimensional relaxed micromorphic model of antiplane shear analyzed in [66], where a gradient field as a right-hand side leads to a gradient field for the microdistortion.

Mixed form
A major aspect of the relaxed micromorphic continuum is the relation to the classical Cauchy continuum theory. This relation is governed by the material constants, where the characteristic length L c plays a significant role.
We are therefore interested in robust computations with respect to L c . To that end, we reformulate the problem as a mixed formulation. The first step consists in introducing the new unknown D = µ macro L 2 c Curl P ∈ H (Div, Ω) , (3.16) reminiscent of the micro-dislocation, and examining its distribution with a test function In fact, D must be solenoidal as and thus the appropriate space for D is We again assume that Dirichlet boundary conditions are prescribed on the whole boundary and thus, from P ∈ H 0 (Curl, Ω) there follows D ∈ H 0 (Div, Ω) as Dν = Curl P ν = 0.
The (bi-)linear forms are now given by where the Lagrange multiplier D has the physical meaning of a hyperstress. Notice we now approximate the hyperstress directly with its own variable, therefore recovering lost precision due to differentiation.
Therefore, we have to prove that a(·, ·) is coercive on the kernel of b(·, ·) The last equality follows from the exact property Curl H 0 (Curl, Ω) = • H 0 (Div, Ω). However, we already know that a({u, P }, {δu, δP }) + Curl P , Curl δP L 2 is coercive, compare the proof of Theorem 3.1. This leaves us with the Ladyzhenskaya-Babuška-Brezzi (LBB) condition to be satisfied We choose u = 0 and P such that The existence of such a potential P for D ∈ • H 0 (Div, Ω) follows directly from the theory of stable Helmholtz decompositions (3.23).
Thus, with Brezzi's theorem [11], there exists a unique solution independent of L c for 1 µmacro L 2 c ≤ 1 fulfilling the stability estimate where the constant c does not depend on L c .

Remark 3.4.
In the proof we made use of the exactness property of the de Rham complex by restricting the space of D to be divergence free. Otherwise we cannot prove the LBB condition in the limit L c → ∞. Numerically we observed that this restriction is necessary, otherwise the resulting matrix becomes singular.
The limit case lim L c → ∞ of Eq. (3.21) is well-defined, resulting in the problem: Consequently, at the limit lim L c → ∞ we have Curl P = 0, i.e. the rows of P are gradient fields. This can be observed by considering that δD only tests functions in range(Curl), where the range is fully given by range(Curl) = • H (Div). As such, by the Helmholtz decomposition only the gradient part of the microdistortion P remains. The latter is also clearly observable through the minimization of the energy function, as P = DΨ is needed for finite energies. Therefore, the consistent coupling condition is crucial, prescribing only gradient fields as Dirichlet boundary conditions for P . From the theory of mixed methods we obtain, analogously to the 2D case in [66], quadratic convergence in L c towards the limit case Remark 3.5. Note that for the limit L c → ∞ one finds P → DΨ, inducing finite energies due to Curl DΨ = 0.
Since L c defines a zoom into the microstructure, the latter can be interpreted as the entire domain being the micro-body. Consequently, setting

30)
and as such, taking the divergence of Eq. (2.8b) results in The divergence of the micro-moment Div M can be interpreted as the micro body-force. The latter implies that the limit L c → ∞ defines again a classical Cauchy continuum theory with a finite stiffness C micro , representing the upper limit of the stiffness for the relaxed micromorphic continuum. Further, due to the consistent coupling condition and Eq. (3.30) there holds DΨ = Du (for a thorough derivation see [6]).

Reformulation of the mixed divergence free constraint
The construction of finite elements which are exactly divergence free, H (Div, Ω), is possible for higher polynomials, [68]. However, at least for the lowest order shape functions, one has to use an additional Lagrange multiplier q ∈ [L 2 0 (Ω)] 3 , compare Remark 3.6 below, forcing Div D = 0. We now show how the mixed formulation has to be adapted leading to a method which can be directly implemented. By defining X = X ×[L 2 0 (Ω)] 3 we introduce the adapted bilinear forms and the problem reads: find ({u, P , q}, D) ∈ X × H 0 (Div, Ω) such that Theorem 3.3. Problem 3.33 is uniquely solvable and the solution solves also Problem 3.21.
Proof. The kernel of the bilinear form b is now given by The last equality follows by choosing on the one hand D = Curl P ∈ H 0 (Div, Ω) yielding Curl P = 0 and on the other hand from the de Rham sequence we find that Div (H 0 (Div, Ω)) = L 2 0 (Ω) such that we can choose D ∈ H 0 (Div, Ω) such that Div D = q. With the same argument as in the proof of Theorem 3.2 the kernel coercivity follows. For the LBB condition we use the Helmholtz decomposition

Now, with
Brezzi's theorem we conclude that Problem 3.33 is uniquely solvable with a stability constant independent of L c . Due to the exactness property of the de Rham complex the additional Lagrange multiplier q enforces that D ∈ • H 0 (Div, Ω). Therefore, taking ({u, P }, D) ∈ X × • H 0 (Div, Ω) from the solution solves also Problem 3.21.
Remark 3.6. In case of prescribed Dirichlet boundary conditions on Γ P D for P ∈ H (Curl, Ω), due to the de Rham complex, compatible Dirichlet conditions have to be used for D ∈ H (Div, Ω). If Γ P D = ∂Ω then from P ∈ H 0 (Curl, Ω) there follows D ∈ H 0 (Div, Ω). Thus, with Div(H 0 (Div, Ω)) = [L 2 0 (Ω)] 3 the space for q has a zero mean. Consequently, another vector valued variable in R 3 is needed for numerics to ensure the zero average over the domain of q, since elements of L 2 cannot be prescribed on the boundary.

Discrete case
We note that existence and uniqueness in the primal formulation is given by the Lax-Milgram theorem. As a result, the use of a commuting diagram for the relation between the displacement u and the microdistortion P is to also be satisfied in the discrete case. Consequently, we rely on the commuting de Rham diagram for the construction of our finite elements. Specifically, for the lower order elements we make use of linear and quadratic tetrahedral Lagrangian elements, linear elements from the first and second Nédélec spaces and lowest-order Raviart-Thomas elements, see Fig. 3. We denote in the following On each tetrahedral element T we define P p (T ) as the set of all polynomials up to order p, P p (T ) = {x i y j z l i + j + l ≤ p, i, j, l ≥ 0} and denote with P p (Ω) the set of piece-wise polynomials of order p on the triangulation of Ω.
Remark 3.7. The Raviart-Thomas element from the space RT 0 can only produce constant normal projections on an element's outer surface. Further, the dimension of the space is smaller than the full linear space, The lowest-order Brezzi-Douglas-Marini elements BDM 1 have linear normal components as the full linear polynomial space is considered. On the one hand better L 2 -estimates are therefore possible, however, on the other hand the range of the divergence operator is not increased. There holds the sequence From Section 3.1.3 we know that the hyperstress field D ∈ H (Div, Ω) is solenoidal, Div D = 0. Therefore, aside from the lowest order case RT 0 , it is desirable to use a basis for RT p (or BDM p ) where all higher-order shape functions are divergence-free. We denote the adapted space by RT p 0 = {d ∈ RT p | d ∈ RT 0 or div d = 0} and the corresponding interpolation operator as Π p d . With such a construction, we can on the one hand save several redundant degrees of freedom, reducing the computational costs for assembly and solution steps without any loss of accuracy, and on the other hand, only a constant correction term q ∈ [P 0 (Ω)] 3 is required to enforce Div D = 0 for the remaining lowest-order RT 0 shape functions in the mixed formulation since In the numerical examples we use NGSolve for the high-order elements where such high-order divergence-free shape functions exist following the construction of [68].
Remark 3.8. The polynomial order of the sequence can be increased by using the Nédélec elements of the first type N 1 I instead of N 1 II , increasing the range of the curl differential operator. We note that both types yield a tangential projection of the same polynomial power on the outer boundaries of the element. Similar to the Raviart-Thomas and Brezzi-Douglas-Marini elements there holds In the Lax-Milgram setting we define X h = V h ×U h . Due to X h ⊂ X the solvability of the discrete problem follows directly from the continuous case [10,15]. We use Cea's lemma for the quasi-best approximation as a basis for the convergence estimate Figure 3: The de Rham complex for linear (upper) and quadratic (lower) sequences. The differential and interpolation operators commute between the continuous and the discrete spaces.
Remark 3.9. When using the canonical interpolation operators in the de Rham complex depicted in Figure 3 additional smoothness of the function spaces has to be assumed, as e.g., point evaluation is not well-defined for functions u ∈ H 1 (Ω), Ω ⊂ R d in dimension d > 1. First commuting, but non-local projections without additional regularity assumptions were constructed in [62] and [14]. Very recently, local L 2 -bounded and commuting projection operators from the function spaces into the finite element spaces have been established in [4]. Therefore, in the theoretical proofs we will make use of these novel projection operators, whereas for the finite element implementations the usual canonical operators are considered.
Proof. During the proof Π I,p−1 c denotes, with an abuse of notation, the projection operator without extra regularity assumptions, see Remark 3.9. With Cea's lemma and the approximation properties of the interpolation operators [11,15] we obtain where | · | H p denotes the standard Sobolev semi-norm.
When using the Nédélec elements of second type N p−1 II instead of N p−1 I , we lose one order of convergence for the microdistortion tensor due to the decreased range of the curl. This should also lead to sub-optimal convergence of the displacement. For example, the quadratic sequence depicted in Figure 3 would only have linear convergence for u although quadratic elements are used. In the extreme case of vanishing characteristic length L c = 0 the optimal convergence rates are achieved since only the L 2 -norm of P is considered, compare Section 3.1.2. On the subspace U 0 = {v ∈ H (curl, Ω) | curl v = 0} the Nédélec elements of first and second kind coincide, N p I = N p II , leading to optimal rates, which corresponds to the limit case L c = ∞. However, in the numerical experiments we continued to observe improved convergence rates for the L 2 -norm of P and H 1 -norm of u for L c ∈ (0, ∞). Under the assumption of an s-regular problem we can prove improved estimates for these norms with an Aubin-Nitsche technique. We call Problem 3.3 s-regular if there holds for the solution {u, P } ∈ X with s ∈ (0, 1] where Proof. We solve the following dual problem, where the differences u − u h and P − P h are used as right-hand side: find {w, Q} ∈ X such that By using the test-function Using Young's inequality on the left side, dividing through u − u h L 2 + P − P h L 2 , and using that {u − u h , P − P h } X ≤ ch p−1 for P h ∈ [N p−1 II ] 3 yields the claim for the L 2 -norm of u and P . For the H 1 -norm we first note that with the Galerkin-orthogonality, adding and subtracting Π p g u, and the specific structure of the bilinear form we obtain Thus, we can estimate ) where we used Cauchy-Schwarz and the explicit structure of the bilinear form. With the triangle inequality and positive definiteness of C e + C c (in combination with the generalized Korn's inequality for incompatible tensor fields if C c = 0) the claim follows.
Remark 3.10. We do not discuss the s-regularity of Problem 3.3 in this work, but mention that for the Laplace as well as for Maxwell's equations there holds s = 1 for convex domains Ω.
Corollary 3.1 also shows that we cannot expect improved convergence for u in the L 2 -norm compared to the H 1 -norm, which is confirmed by the numerics.
Unfortunately, the approximation constant depends a priori on L c . As we are interested in the limit L c → ∞ robust estimates with respect to L c are desirable. We show such an estimate with the help of the equivalent mixed formulation of the problem.
Normally, in contrast to the primal formulation, the kernel coercivity as well as the LBB-condition for the discrete case are not inherited from the continuous proof. However, thanks to the de Rham complex and a discrete Helmholtz decomposition [3], these properties hold also in the discrete case as long as the appropriate conforming finite element spaces are employed [11]. Consequently, Cea's lemma yields Additionally, with {u ∞ , P ∞ , D ∞ , q ∞ } for the (smooth) solution of the limit problem one obtains Proof. The proof of the first claim starts by following the same lines as in the primal case. The estimates for D and q are shown in detail where we used that the exact hyperstress D is divergence free. This estimate, however, would only yield linear convergence for all fields. By noting that on the subspace of divergence free hyperstresses the solution is independent of q we obtain the optimal rates for u h , P h , and D h . For the second inequality we add and subtract the exact solution, use the triangle inequality, and (3.29) and (3.40).
Remark 3.11. Note that only linear convergence is obtained for the correction term q. If one is interested in this quantity, a cheap post-processing step can be applied element-wise to obtain a high-order approximation [30]. For U h = [N p−1 II ] 3 the results of Lemma 3.2 can directly be adapted by following the same lines as in Corollary 3.1. Inequalities (3.41) show that besides the model error in L c , also the approximation error has to be considered for convergence studies in the limit L c → ∞.

Finite element formulations
In this work we employ tetrahedral finite elements. The construction of H 1 , H (div), H (curl), and L 2 -conforming finite elements to obtain the lowest order combination is presented in the following in detail. Higher order elements are used from the open source finite element software NETGEN/NGSolve 4 [61,63] and we refer to [68] for the construction of higher order finite elements. The elements are mapped from the reference element to the physical element using the barycentric base functions

Nédélec base
For the interpolation of the microdistortion tensor P we make use of Nédélec base functions of the first and second types [40,41,68]. The edge degrees of freedom are defined by the functionals where s i is the curve of an edge on the element and q j is a test function. In accordance with the de Rham complex, the appropriate polynomial space for the linear sequence is where P is the space of homogeneous polynomials. The corresponding test function q j = 1 yields the base functions on the reference element (see Fig. 6) representing the base of the lowest order Nédélec element of the first type N 0 I . For the sequence of the quadratic element we employ the linear Nédélec elements of the second type N 1 II . The corresponding polynomial space reads In other words, the space is given by 12 Nédélec functions. Using the ansatz 15) and the degrees of freedom from Eq. (4.11) for each edge on the reference element where µ i is the parameter for the corresponding edge, we find the base functions (see Fig. 7) Remark 4.1. The resulting Nédélec base functions are not unique in the sense that they are determined by the choice of the test functions q j . In our case we chose the linear Lagrangian base on µ, leading to a Nédélec-Lagrangian base.
Remark 4.2. Note that the decrease in polynomial order of the second discrete sequence (see Fig. 2) can be alleviated by employing first order Nédélec elements of the first type N 1 I . This increases the interpolation power of the solenoidal part of the microdistortion field P , but not of the irrotational part.
The base functions are defined on the reference element. In order to preserve their tangential projection on edges in the physical domain to find the covariant Piola transformation [38] θ, ds = θ, J dµ = ϑ, dµ ⇐⇒ θ = J −T ϑ . Further, for the transformation of the curl we find which is the so called contravariant Piola transformation. The formula is derived using the identity 5 ∇ x × J −T = 0. The microdistortion can now be interpolated by Θ P e , (4.22) In other words, the microdistortion is interpolated by 18 Nédélec base functions in the linear case and 36 in the quadratic case on each element.

Raviart-Thomas base
For the mixed formulation we also require Raviart-Thomas [53,68] elements. The lowest order face degrees of freedom of Raviart-Thomas elements are defined as The polynomial space for the construction of base functions is given by where P is the space of homogeneous polynomials. In the lowest order one finds Using the ansatz (4.28) and the degrees of freedom from Eq. (4.25) for each face on the reference tetrahedron where Λ denotes a face on the reference tetrahedron, we find the base functions (see Fig. 8) we consider the transformation of surfaces Consequently, the base functions are mapped using the contravariant Piola transformation Remark 4.3. Note that this the same transformation as in Eq. (4.21) since curl θ ∈ H (div) , θ ∈ H (curl).
Considering the divergence of functions undergoing a contravariant Piola transformation we observe As a result, one finds Thus, the hyperstress field is interpolated by such that, on each element 12 Raviart-Thomas base functions define the hyperstress.

Discontinuous basis
Finally, for the discontinuous elements in L 2 we employ piece-wise constants The entire space has the dimension dim[P 0 ] 3 = 3 and the interpolation reads

The orientation problem
The co-and contravariant Piola transformations do not suffice to assert the consistent orientation of the tangential or normal projections of the Nédélec and Raviart-Thomas base functions, respectively. The transformations control the size of the projections, but not whether these are parallel or anti-parallel with respect to neighbouring elements. Consistent projections is a key requirement in ensuring no jumps occur in the trace of the respective space and as such, there exist various methods for dealing with this so called orientation problem [1, 2,66,68]. In this work we present a solution based on the sequencing of vertices and the separation of orientational data. We define the following rule for the orientation of edges This means each edge starts at the lower vertex index and ends at the higher vertex index. This definition determines the orientation of the edge tangent vector (see Fig. 4) and consequently, the tangential projection of the Nédélec base functions. Analogously, for surfaces we define such that each surface is given by a sequence of increasing vertex indices. The orientation of the surface normal is given according to the left-hand rule. In other words, the direction of the normal is determined by the cross product of the vectors arising from the edges The orientation of the Nédélec-and Raviart-Thomas base functions is according to these rules. Consequently, in order to map each tetrahedron in the mesh to this orientation, we define each element as an increasing vertex-index sequence (see Fig. 9) The latter ensures the consistent orientation of the base functions, since they are all mapped from the same reference domain. However, integration in the reference element is determined by the determinant of the Jacobi matrix which may be negative due to a reflection of the element in the mapping from the reference to the physical domain. We correct for the error by taking only the absolute value of the determinant Consequently, consistency is guaranteed by mapping from a single reference element and the use of correction functions or considerations of neighbouring elements are circumvented.
Remark 4.5. The absolute value of det J is only used for the integration over the element. In all other use-cases, the information of the sign is necessary.

The discrete consistent coupling condition
In order to exactly satisfy the consistent coupling, one may use the degrees of freedom from [20] to set the boundary conditions, ensuring a commuting projection. However, this requires the solution of a Dirichlet type problem for every edge on Γ D . Alternatively, one can first find the Lagrangian interpolation of the boundary condition of the displacement on each edge and derive the boundary condition for the microdistortion directly from it. We observe that on each element's edge, the displacement is given by where n i are the Lagrangian base functions on the reference edge and m is the mid-point of the edge. The consistent coupling condition on a curve transforms according to the covariant Piola transformation where p represents a row of the microdistortion P and Π c p, ς results from the transformations of the tangent vector Eq. (4.19) and the Nédélec base functions Eq. (4.20). The gradient of the Lagrangian interpolation of the displacement Π g u can be reformulated in dependence of the curve parameter µ ∈ [0, 1] We note that in order to maintain an invariant transformation, the tangent vectors on the reference domain are given by and are not necessarily unit vectors, consider edges Fig. 4. Consequently, one finds (Π g u) ,µ D ξ µ, ς = (Π g u) ,µ D ξ µ, ξ ,µ = (Π g u) ,µ . (4.49) As a result, we can simplify the consistent coupling condition to For the linear element, each edge has one Nédélec base function of the first type with a constant tangential projection allowing to easily embed the consistent coupling condition Consequently, on each boundary edge the interpolation of one row of the microdistortion is given by In the quadratic sequence, each edge is associated with two Nédélec base functions of the second type ϑ v1 , ϑ v2 . Their projections on the tangent vector at their respective vertices is one and zero on the other vertex Consequently, it suffices to control the tangential projection of each base function at its respective vertex in order to impose the boundary condition (4.56) The resulting interpolation of one row of the microdistortion at the boundary edge reads For both sequences the above demonstrated methodology satisfies the discrete consistent coupling condition exactly by construction.

Element stiffness matrices
The matrices N , Θ, Φ, and 1 define the interpolation for each element. The resulting finite element stiffness matrices for the primal formulation read The linear element's stiffness matrix has the dimensions K e ∈ R 30×30 and the matrix of the quadratic element the dimensions K e ∈ R 66×66 . The load vector of the primal formulation is given by with the dimensions f e ∈ R 30 for the linear sequence and f e ∈ R 66 for the quadratic sequence, respectively. The mixed formulation induces the following additional stiffness matrices where K uu , K uP , and K P u remain unchanged with respect to the primal formulation. The linear sequence yields the matrix dimensions K e ∈ R 45×45 , whereas the quadratic sequence results in K e ∈ R 81×81 . The load vector changes its dimensions in the mixed formulation, but not its content with dimensions corresponding to the element stiffness matrix.
Remark 4.6. Note that the symmetry and anti-symmetry operators have been dropped. This is because the material tensors are defined as C e , C micro : R 3×3 → Sym(3) and C c : R 3×3 → so (3), implying the former. For materials where this is not guaranteed, the operators can be recovered by modifying the material tensors C e sym P = C e S P = C e P , C micro sym P = C micro S P = C micro P , C c skew P = C c A P = C c P . (4.62)

Numerical examples
In the following we test the convergence rates of the model. The behaviour with respect to the characteristic length L c is crucial, as this parameter governs the impact of the micro-structure on the displacement field u and the relation to the Cauchy continuum.
In order to study convergence we construct artificial analytical solutions for the displacement field u and the microdistortion field P (see Appendix A) and compare them with the numerical approximation of the model. Errors are measured in the L 2 -norm

Convergence for L c → 0
In this benchmark the domain is defined as the axis-symmetric cube Ω = [−1, 1] 3 . For simplicity the material constants in the isotropic case are set to λ e , µ e , µ c , λ micro , µ micro , µ macro = 1 .
The characteristic length is varied between L c ∈ (0, 1] in order to test the stability of the formulations. The prescribed fields are defined as being of higher polynomial order than the quadratic element. The entire boundary is set to ∂Ω = Γ D . By applying the strong form one finds the force and micro-moment Reducing the characteristic length L c or setting it to zero does not disturb the stability of the computation, even though control of the Curl of P is lost, as proven by the convergence rates in Fig. 10. This is a direct result of the derivation in Section 3.1.2, asserting stable computations for P ∈ [L 2 (Ω)] 3×3 .
In an alternative example we take the Cosserat couple modulus µ c = 0 (i.e., C c = 0), resulting in the force 6) and the micro-moment The primal formulation maintains its stability up to small L c -values coupled with small elements, see Fig. 11. For L c = 10 −10 the formulation becomes completely unstable and the matrix singular. The deterioration in the convergence is clearly depicted in Fig. 12. Since µ c = 0 and L c → 0, no term controls the skew-symmetric part of the microdistortion and the formulation only converges for the displacement u, up to L c = 10 −10 , where the stiffness matrix becomes singular.

Robustness in L c
To test for robustness of the formulation in the characteristic length L c we again consider Ω = [−1, 1] 3 , the parameters λ e = µ e = µ c = λ micro = µ micro = µ macro = 1 and µ c = 0, (5.8) and the following exact solution fields  where the first part of P is curl-free. The second term is scaled by L −2 c to obtain a well-defined limit for L c → ∞. The corresponding right-hand sides read We use NGSolve in combination with the direct solver UMFPACK [18] for inverting the arising stiffness matrix for the primal and mixed method. As can be clearly observed in Figure 13 for larger values of L c ≥ 10 6 the primal method becomes unstable, whereas the mixed method still converges optimally.

Convergence for L c → ∞
Next, we consider the convergence to the limit solution L c → ∞. Rewriting the exact solution (5.9) in the mixed formulation (3.28) and taking the limit L c → ∞ yields the following displacement, microdistortion, and hyperstress fields (5.10) Note, that Curl P = 0 and Div D = 0. The hyperstress field D acts as Lagrange multiplier forcing the microdistortion tensor to be curl-free. By taking the limit L c → ∞ we obtain directly the forces We consider five different structured tetrahedral grids consisting of 6, 48, 384, 3072, and 24576 elements and solve the primal and mixed formulation with NGSolve and UMFPACK using different values for L c . As depicted in Figure 14 the expected quadratic convergence rate in L c (3.29) is observed until the discretization error is reached (3.41). For finer grids we can reduce this error. However, the primal method again becomes unstable for large values of L c . The mixed method, where also the limit problem with L c = ∞ is well defined, is again completely stable with respect to L c .
Assuming micro to be always stiffer than macro, we set the macro and micro parameters λ macro = 115.4 , µ macro = 76.9 , λ micro = 10 λ macro = 1154 , µ micro = 10 µ macro = 769 ,  The behaviour of the deviation in both the displacement and the energy in Fig. 15 and Fig. 16 demonstrates the derivation in Eq. (3.14). Further, it is apparent that this characteristic strongly depends on the approximation capacity of the finite subspace. In other words, the use of p-refinement greatly influences the result. Decreasing the element size (h-refinement) can also be used to alleviate the error.

Bounded stiffness
In the example in Section 5.4 one can observe the influence of the characteristic length L c on the reaction of the system. In fact, the characteristic length acts as an interpolation parameter between the micro-stiffness and the macro-stiffness, as shown in Remark 3.5 and Section 3.1.2. In order to demonstrate this phenomenon, we employ the axis-symmetric cubic domain Ω = [−1, 1] 3 and the material parameters from the previous example. The characteristic length L c is varied in [10 −3 , 10 3 ]. The Dirichlet boundary is set to Note that while the problem is reminiscent of simple shear [57], it is different since we are using a finite volume and only the upper and lower parts of the boundary are constrained. The problem is defined entirely as a boundary value problem and no body forces or micro-moments are applied. In order to extract an indicator for the stiffness of the system we measure the Lebesgue integral of the displacement reaction forces in x-direction on the upper Dirichlet boundary where t x represents the equivalent traction in x-direction (see Appendix B). Considering Fig. 17, it is clear that the linear finite element formulation largely overestimates the stiffness of the system by about a factor of 28% in the lower limit for the finer discretization. Further, the effects of mesh refinement are considerable. From the depiction of the displacement field we recognize a transition from a nearly linear displacement field for large L c values to a higher order displacement function for low L c values. Lastly, we observe that the intensity of the microdistortion field (first row of P ) shifts from the centre to the sides of the domain for a decreasing characteristic length L c .
Remark 5.1. The lower limit C macro was estimated using 6000 finite elements of the quadratic sequence in the primal formulation while setting the characteristic length to L c = 0. In order to approximate the upper limit C micro , we make use of the quadratic sequence in the mixed formulation with L c = 10 9 and 6000 elements, since the mixed formulation is stable for very large values of the characteristic length.

Conclusions and outlook
Existence and uniqueness in the primal formulation of the relaxed micromorphic model requires the employment of Nédélec elements, since N p I , N p II ⊂ H (curl). The construction of the lowest order Nédélec elements was demonstrated along with a solution to the arising orientation problem, such that no correction functions are necessary and a single reference element suffices.
While not being a mixed formulation, the use of the exact de Rham sequence is still recommended in order to exactly satisfy the discrete consistent coupling condition, as shown in Section 4.6. This is otherwise not possible for the general case if the exact sequence is not respected.
The example in Section 5.4 depicted the behaviour of the relaxed micromorphic theory with L c → 0 as given by the derivation in Eq. (3.14). Clearly, the characteristic length L c is a crucial component of the theory as it governs the relation to the classical Cauchy continuum. Further, it is demonstrated that the latter comparison is strongly dependent on the approximation capacity of the finite subspace. Due to the relatively large error with respect to the expected result using the linear finite element sequence, u Cauchy − u relaxed L 2 / u Cauchy L 2 ≥ 9.5%, it is recommended to use at least quadratic or higher order finite elements. This conclusion is further supported by the example in Section 5.5, where the linear formulation largely overestimated the expected stiffness of the system for L c → ∞.
While the primal formulation is relatively robust in L c , it may become unstable as L c → ∞. As demonstrated in the example in Section 5.2, the computation can be re-stabilised by using a mixed formulation. As per the de Rham complex, this requires the usage of H (div) elements, such as Raviart-Thomas RT or Brezzi-Douglas-Marini BDM, and completely discontinuous finite elements. Further, for efficient simulations with higher order   Figure 17: (a) Reaction forces with varying characteristic lengths L c , (b) initial geometry with 3072 elements, (c) cutout view of the displacement field of the quadratic element with L c = 10 3 , (d) cutout view of the displacement field of the quadratic element with L c = 1, (d) cutout view of the displacement field of the quadratic element with L c = 10 −3 , (f) microdistortion field of the quadratic element with L c = 10 3 , (g) microdistortion field of the quadratic element with L c = 10 −3 .
elements it is advantageous to split the H (div) elements into source and solenoidal parts as the hyperstress field D is completely solenoidal, meaning the source part of the interpolant should be dropped or must be otherwise compensated. The construction of the lowest order Raviart-Thomas elements along with the lowest order discontinuous elements was demonstrated in Section 4.3 and Section 4.4, respectively. The orientation problem for the latter elements is circumvented by the same methodology as for the Nédélec elements. For further development we note the computation time increases rapidly with more elements and higher order polynomials. The development of an efficient preconditioner is to be considered for future works.

A Derivation of the strong form and analytical solutions
The variation of the free energy functional with respect to the displacement field u reads δ u I = Ω C e sym Dδu, sym(Du − P ) + C c skew Dδu, skew(Du − P ) − δu, f dX ∀ δu ∈ C 1 Γ D (Ω, R 3 ) , (A.1) where C 1 Γ D (Ω, R 3 ) denotes the set of differentiable functions which are zero at the Dirichlet boundary. By using the Green-type identity where P × ν on Γ P D is directly embedded into the space. Consequently, the strong form reads

B Equivalent traction
The traction on the Neumann boundary is defined using Eq. Remark B.1. Note that using the latter, no new computation of the stiffness matrix K, the nodal values vector u d or the force vector f d is required.