{'fraction_non_alphanumeric': 0.04184190330007675, 'fraction_numerical': 0.020598618572524943, 'mean_word_length': 4.233933161953727, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mimetic Spectral Element Method for Anisotropic Diffusion

This paper addresses the topological structure of steady, anisotropic, inhomogeneous diffusion problems. Two discrete formulations: a) mixed and b) direct formulations are discussed. Differential operators are represented by sparse incidence matrices, while weighted mass matrices play the role of metric-dependent Hodge matrices. The resulting mixed formulations are point-wise divergence-free if the right hand side function f = 0. The method is inf-sup stable and displays optimal convergence on orthogonal and non-affine grids.


Introduction
Anisotropic and inhomogeneous diffusion appears in many applications such as heat transfer [15], flow through porous media [87], turbulent fluid flow [116], image processing [98] or plasma physics [112]. In 2D, steady, anisotropic diffusion is governed by the following elliptic partial differential equation Here, p is the flow potential, f the source term, with p =p along Γ p and (K∇p, n n n) = u n along Γ u . Here, for all x x x, K(x x x) is a symmetric, positive definite tensor.
In the presence of strong anisotropy, i.e. large ratio between the smallest and largest eigenvalues of the diffusion tensor, the construction of robust and efficient discretizations becomes particularly challenging. Under these conditions, the convergence rates of the discretization error can be considerably reduced; this effect is commonly referred in the literature as locking effect, see for example [4,5,12,84]. For sufficiently refined discretizations, the deterioration of the convergence rates eventually disappears. Unfortunately, this may occur only when the grid cell size is prohibitively small.
Another important aspect is mesh flexibility. In many applications of diffusion equations, particularly in porous media flow, typical grids are highly irregular. In many of these situations the results obtained are strongly dependent on the grid type, see [11] for a discussion of the use and properties of different grids in reservoir modelling.

Overview of standard discretizations
In order to overcome these limitations and improve the efficiency and robustness of the discretization of the anisotropic diffusion equations, several approaches have been proposed.
The discretization of the anisotropic diffusion equations in complex media in many situations is still a trade-off between, e.g. [89]: • Accuracy in the representation of the medium (complex grids).
• Accuracy in the discretization of the equations.
The need for such a choice is rooted in the use of numerical schemes based on two-point flux approximations (TPFA), see for example, [3,89,120]. These methods produce good approximations on orthogonal grids when the diffusion tensor K is diagonal, but are known to introduce significant discretization errors in the presence of a non-diagonal diffusion tensor. This introduces severe limitations into the possible grid choices. Under these conditions, the geometric flexibility introduced by perpendicular bisector (PEBI) grids, [11,67,90], is considerably limited, for example.
It has been known that the discretization error is related to the misalignment between the grid and the principal directions of the diffusion tensor K. In fact, Aavatsmark showed in [3] that for TPFA this misalignment leads to the discretization of the wrong diffusion tensor.
These ideas initially led to the construction of grids aligned with the principal axis of the diffusion tensor, so called K-orthogonal grids, see for example [65,67]. This approach significantly improves the performance of the numerical method but substantially limits the geometric flexibility.
More recently, multipoint flux-approximation (MPFA) schemes have been introduced specifically to address these limitations, see e.g. the initial works by Aavatsmark [4,5] or a more recent presentation [2], and by Edwards and Rogers [57]. This method is based on a cell-centred finite volume formulation and introduces a dual grid in order to generate shared sub-cells and sub-interfaces. This in turn produces a discretization of the flux between two cells that involves a linear combination of several adjacent cells. This method is robust and locally conservative but does not guarantee a resulting symmetric discrete diffusion operator. More recently, this work has been connected to the mixed finite element method, [56].
Alternative approaches based on the finite element formulation have also been proposed by several authors. We briefly mention the work on the control-volume finite element discretization by Forsyth [60] and Durlofsky [54], on nodal Galerkin finite elements by Young [122], and on mixed finite elements by Durlofsky [53].

Overview of mimetic discretizations
Over the years, the development of numerical schemes that preserve some of the structures of the differential models they approximate has been identified as an important ingredient of numerical analysis. One of the contributions of the formalism of mimetic methods is to identify differential geometry as the proper language in which to encode these structures/symmetries. Another novel aspect of mimetic discretizations is the identification and separation of physical field laws into two sets: (i) topological relations (metric-free), and (ii) constitutive relations (metric dependent). Topological relations are intimately related to conservation laws and can (and should) be exactly represented on the computational grid. Constitutive relations include all material properties and therefore are approximate relations. For this reason, all numerical discretization error should be included in these equations. A general introduction and overview of spatial and temporal mimetic/geometric methods can be found in [39,42,66,100].
The relation between differential geometry and algebraic topology in physical theories was first established by Tonti [117]. Around the same time Dodziuk [52] set up a finite difference framework for harmonic functions based on Hodge theory. Both Tonti and Dodziuk introduce differential forms and cochain spaces as the building blocks for their theory. The relation between differential forms and cochains is established by the Whitney map (k-cochains → k-forms) and the de Rham map (k-forms → k-cochains). The interpolation of cochains to differential forms on a triangular grid was already established by Whitney, [119]. These generalized interpolatory forms are now known as Whitney forms.
Hyman and Scovel [75] set up the discrete framework in terms of cochains, which are the natural building blocks of finite volume methods. Later, Bochev and Hyman [18] extended this work and derived discrete operators such as the discrete wedge product, the discrete codifferential, and the discrete inner products. Robidoux,Hyman,Steinberg and Shashkov,[74,[76][77][78]107,108,111,113,114] used symmetry considerations to construct discretizations on rough grids, within the finite difference/volume setting . In a more recent paper by Robidoux and Steinberg [110] a finite difference discrete vector calculus is presented. In that work, the differential operators grad, curl and div are exactly represented at the discrete level and the numerical approximations are all contained in the constitutive relations, which are already polluted by modeling and experimental error. For mimetic finite differences, see also the work of Brezzi et al. [32,33] and Beirão da Veiga et al. [47].
The application of mimetic ideas to unstructured triangular staggered grids has been extensively studied by Perot, [99,[101][102][103]123], specially in [100] where the rationale of preserving symmetries in numerical algorithms is well described. The most geometric approach is presented in the work by Desbrun et al. [49,58,86,97] and the thesis by Hirani [72].
The Japanese papers by Bossavit, [26][27][28][29][30], serve as an excellent introduction and motivation for the use of differential forms in the description of physics and the use in numerical modeling. The field of application is electromagnetism, but these papers are sufficiently general to extend to other physical theories.
Another approach to develop a discretization of the physical field laws is based on a discrete variational principle for the discrete Lagrangian action. This approach has been used in the past to construct variational integrators for Lagrangian systems, e.g. [79,85]. Kraus and Maj [80] have used the method of formal Lagrangians to derive generalized Lagrangians for non-Lagrangian systems of equations. This allows to apply variational techniques to construct structure preserving discretizations on a much wider range of systems. Recently, Bauer and Gay-Balmaz presented variational integrators for elastic and pseudo-incompressible flows [14].

Outline of chapter
In Section 2 the topological structure of anisotropic diffusion problems is discussed. In Section 3 spectral basis functions are introduced which are compatible with the topological structure introduced in Section 2. In Section 4 transformation to curvilinear elements is discussed. Results of the proposed method are presented in Section 5.

Anisotropic diffusion / Darcy problem
Let Ω ⊂ R d be a contractible domain with Lipschitz continuous boundary ∂ Ω = Γ p ∪ Γ u , Γ p ∩ Γ u = ∅. The steady anisotropic diffusion problem is given by with p =p along Γ p and (K∇p, n n n) =ū n along Γ u . Here, for all x x x, K(x x x) is a symmetric, positive definite tensor, i.e. there exist constants α,C > 0 such that If Γ p = ∅, then (2) has a unique solution. If Γ p = ∅ then (2) only possesses solutions if ∂ Ωū in which case the solution, p, is determined up to a constant. An equivalent first order system is obtained by introducing u u u = −K∇p in which case (2) can be written as An alternative first-order formulation is given by Formulation (3) is generally referred to as the Darcy problem, while the relation u u u = −K∇p is called Darcy's law, [87]. The Darcy problem plays an important role in reservoir engineering. In this case u u u is the flow velocity in a porous medium and p denotes the pressure.
While the formulations (2), (3) and (4) are equivalent, (2) only has 1 unknown, p, (3) has (d + 1) unknowns, p and the d components of u u u, and (4) has (2d + 1) unknowns. Formulation (4) is of special interest, because it decomposes the anisotropic diffusion problem into two topological conservation laws and one constitutive law 1 . By making a suitable choice where and how to represent the unknowns on a grid, the topological relations, v v v − ∇p = 0 and ∇ · u u u = f reduce to extremely simple algebraic relations which depend only on the topology of the mesh and are independent of the mesh size, independent of the shape of the mesh, and independent of the order of the numerical scheme. We will refer to such discretizations as exact discrete representations.

Gradient relation
Consider two points A, B ∈ Ω and a curve C which connects these two points, then where dl is a small increment along the curve C .
Suppose that we take another curveC which connects the two points A and B then we also havev v vC : The integral along C is equal to the integral alongC . We will refer tov v v as an integral value, since it denotes an integral and not a point-wise evaluation of v v v. The advantages of integral values are: 1 An even more extended system is, see for instance [16] This seems an unnecessarily complicated system. If we eliminate ψ from (5) we obtain (4). The usefulness of this system lies in the fact that by introducing ψ, the conservation ∇ · u u u = f becomes independent of the data of the PDE, in this case the right hand side function. A similar situation occurs when K = I, the identity tensor, then the equation u u u+Kv v v = 0 in (4) seems redundant, but we have good reason to keep this seemingly redundant equation as we will show in this paper. 1. The velocity-gradient relation is exact. It is not obtained by truncated Taylorseries expansions or does not depend on the choice of basis functions/interpolations. 2. Does not depend on mesh parameters. The mesh size h does not appear in (6).
Whether the curve which connects two points is straight or curved is irrelevant in this relation, therefore this relation is directly applicable on curved domains. 3. Integral quantities are additive.
Consider the points and lines segments as shown in Figure 1. In this figure the arrow along the curves indicates the direction in which v v v is integrated 2 . Application of (6) shows, for instance, thatv v v 14 = P 6 − P 2 .
The additivity property implies that and even more paths can be constructed that connect P 2 and P 7 . The independence of the path depends critically on the assumption that the space is contractible, i.e. there are no holes in the domain (Poincaré's Lemma). A special case is the curve from a point to itself, say P 2 → P 2 in Figure 1. This integral is zero and if the integral is independent of the path this implies that, for instance, where we once again use the additivity property. We see that the circulation vanishes if v v v is a potential flow, which in turn implies that the circulation of the velocity field over the boundary of any surface vanishes. Or, using Stokes' theorem, the integrated vorticity w vanishes. Here the vorticity w is represented as the integral over a surface. We can collect all the integrated velocity fields and pressures in Figure 1 in the following form If we store allv v v i in a vector v and all P j in a vector P and denote the matrix by E 1,0 , we have v = E 1,0 P . If we now also collect all the integrated vorticities, w i , we can relate them to the integrated velocities in the following way If we store all vorticity integrals, w i in the vector w, then we can write this as The matrices E 1,0 and E 2,1 are called incidence matrices. We have E 2,1 · E 1,0 ≡ 0. This identity holds for this particular case, but is generally true; it holds when we would have used triangles or polyhedra instead of quadrilaterals and it holds in any space dimension d. If E 1,0 represents the gradient operation and E 2,1 represents the curl operation, then E 2,1 · E 1,0 ≡ 0 is the discrete analogue of the vector identity ∇ × ∇ ≡ 0, [24-27, 49, 88, 110]. If boundary conditions for p are prescribed along ∂ Ω , then these degrees of freedom can be removed from the grid in Figure 1.
If p is known along the boundary then the integral of v v v is also known along the boundary, so the degrees of freedom for v can also be removed. Relabeling the remaining unknowns gives the geometric degrees of freedom as shown in Figure 2.

Divergence relation
Consider a bounded, contractible volume V ⊂ Ω then we have If the boundary ∂ V can be partitioned into n sub-boundaries, Collecting all fluxes and source terms in vectors u and f, respectively, we can write this equation The matrixẼ 2,1 is the incidence matrix which represents the divergence operator, not to be confused with E 2,1 in (8) which represents the curl operator. If, in the 2D case, the flow field is divergence-free, i.e. f = 0, we know that a stream function ψ ψ ψ exists which is connected to u u u by If we represent the stream function in the nodes of the grid shown in Figure 3, then we have the exact topological equation We can write this in terms of incidence matrices as 3 If the flux u u u is prescribed along the Γ u the associated edges (2D) or surfaces (3D) can be eliminated from the system E d,d−1 u = f and transferred to the right hand side.
For the discretization of (4) the first and last equation in that system can be represented on the mesh by 3 Note that if we performed the same steps in 3D, then the divergence relation (10) would bẽ E 3,2 u = f , and the 2D stream function becomes the 3D stream vector field and we would have So clearly the incidence matricesẼ depend on the dimension of the space d in which the problem is posed. Note that this is not the case for the incidence matrices E. Alternatively, we could refer to the dimension-dependent incidence matrices as in which case it is immediately clear that these matrices depend on the d. From now on we will use the incidence matrices with the d, because then the results are valid for any space dimension d. v Prescription of boundary conditions p along Γ u and u u u along Γ u can be done strongly.
The degrees of freedom can be eliminated and transferred to the right hand side. The equation between p and v v v is exact on any grid and the discrete divergence relation between u u u and f is exact on any grid. Note the (v v v, p)-grid is not necessarily the (u u u, f )-grid, so in principle we can use different grids for both equations. Unfortunately, neither of the two problems, v v v = ∇p and ∇·u u u = f has a unique solution on their respective grids. It is the final equation in (4), u u u = −Kv v v, that couples the solution on the two grids and renders a unique solution. It is also in this equation that the numerical approximation is made; the more accurate we approximate this algebraic equation, the more accurate the solution to the first order system (4) will be.
For many numerical methods 4 well-posedness requires that the number of discrete degrees of freedomv v v i is equal to the discrete number of degrees of freedom u u u j , or more geometrically, that the number of k-dimensional geometric objects on one grid is equal to the number of (d − k)-dimensional geometric objects on the other grid. Here k = 0 refers to points in the grid, k = 1 to edges in the grid, k = 2 the faces in the grid, and k = 3 the volumes in the grid.
The requirement #k = #(d − k) cannot be accomplished on a single grid, so this requires two different grids which are constructed in such a way that #k = #(d − k) holds, [24,25,49,82,88,110]. A dual grid complex is shown in Figure 4. The integral quantities (v v v, p) can be represented on the gray grid. If p is prescribed along the entire boundary, then those degrees of freedom are eliminated (including the gray edges along the boundary for which the integral value v v v is then known also), see for instance Figure 2. In that case flux u u u along the boundary cannot be prescribed. In Figure 4, the number of points in the gray grid, 9, equals the number of surfaces in the black grid, the number of edges in the grey grid is equal to the number of edges on the black grid, 24, and the number of surfaces on the gray grid equals the number of points in the black grid, 16, therefore, we have #k = #(d − k) for d = 2.
Alternatively, we could have represented (u u u, f ) on the gray grid with u u u and the stream function ψ prescribed and (v v v, p) on the black grid.In this case Γ u = ∂ Ω and Γ p = ∅.

Dual grids
If dual grids, such as described above, are employed then we have two properties: 1. There exists a square, invertible matrix H d−1,1 2. The incidence matrices on the primal and dual grid satisfy 5 If we use dual grids and these properties hold, we can write (4) as where the vectors p, v, u and f contain the integral quantities in the mesh as discussed in the previous sections.
In the diagram below, we place the various integral values in appropriate 'spaces' This relation is true if the orientations on primal and dual grid agree. This is not always the case and then the relation reads E d−k,d−k−1 = −E k,k−1 T . A well known example is the duality between grad and div.
Here H k denotes the space of values assigned to k-dimensional objects in the H-grid for k = 0, 1, 2. IfH denotes the dual grid, thenH l is the space of values assigned to l-dimensional objects in theH-grid. For dual grids the number of points in the H-grid is equal to the number of ddimensional volumes in the dual gridH. Let H d,0 and H 0,d be square, invertible matrices which map between H 0 andH d as shown in the diagram above.
If we eliminate v and u from (12) we have This discretization corresponds to (2). We will refer to this formulation as the direct formulation.
If p ∈H d we can set up the diffusion problem as This formulation, where we solve for p and u u u simultaneously, resembles (3), and will be called the mixed formulation, [35].

Mimetic spectral element method
The incidence matrices introduced in the previous section are generic and only depend on the grid topology. The matrices H which switch between the primal and the dual grid representation explicitly depend on the numerical method that is used. In this section we will introduce spectral element functions which interpolate the integral values in a grid. With these functions we can construct the H-matrices, which turn out to be (weighted) finite element mass matrices. The derivation in this section will be on an orthogonal grid. The extension to curvilinear grids will be discussed in the next section.

One dimensional spectral basis functions
The explicit form of the Lagrange polynomials in terms of the Legendre polynomials is given by .
Let f (ξ ) be a function defined for ξ ∈ [−1, 1] by Using property (15) we see that f (ξ j ) = a j , so the expansion coefficients in (17) coincide with the value of f in the GLL nodes. We will refer to this expansion as a nodal expansion, because the expansion coefficients, a i in (17) The functions e i (ξ ) are polynomials of degree (N − 1). These polynomials satisfy, [61,82,96] ξ j Let a function g(ξ ) be expanded in these functions then using (19) ξ j So the expansion coefficients b i in (20) coincide with the integral of g over the edge [ξ i−1 , ξ i ]. We will call these basis functions edge functions and refer to the expansion (20) as an edge expansion, see for instance [16,82,96] for examples of nodal and edge expansions.
Let f (ξ ) be expanded in terms Lagrange polynomials as in (17), then the derivative 6 of f is given by, [61,82,96] If we collect all the expansion coefficients in a column vector and all the basis functions in a row vector we have then the derivative is given by 7 (23) So taking the derivative essentially consists of two step: Apply the matrix E 1,0 to the expansion coefficients and expand in a new basis.

Expanding p (Direct formulation)
In finite element methods the direct finite element formulation for the anisotropic diffusion problem is given by: For (K∇p, n n n) = 0 along Γ u and f ∈ H −1 (Ω ), find p ∈ H 1 0,Γ p (Ω ) such that 6 Note that the set of polynomials {h i }, i = 0, . . . , N is linearly dependent and therefore does not form a basis, while the set {e i }, i = 1, . . . , N is linearly independent and therefore forms a basis for the derivatives of the nodal expansion (17). 7 The matrix E 1,0 is the incidence matrix as was discussed in Sections 2.1 and 2.2. It takes the nodal expansion coefficients and maps them to the edge expansion coefficients. The incidence matrix is the topological part of the derivative. It is independent of the order of the method (the polynomial degree N) and the size or the shape of the mesh. The incidence matrix only depends on the topology and orientation of the grid, see [18,81,82].
The vector f in (28) is given by If we compare (28) with (13), we see that the H d−1,1 K -matrix from (13) is represented in the finite element formulation by the weighted mass matrix M (1) K given by (29), see also [18,115].
In contrast to the pressure expansion in Section 3.2.1 in the direct formulation, (25), in the mixed formulation the pressure is expanded in terms of edge functions The velocity u u u is expanded as Application of the divergence operator to (32) and using (21) we obtain Note that E d,d−1 is the incidence matrix which also appeared in (10) and footnote 3. If we insert the expansion (32) in (ũ u u, K −1 u u u) we obtain with M (d−1) Note that pressure is expanded in the same basis as the divergence of the velocity field, (31) and (33), therefore we can write with With (34) and (37) we can write (30) as Comparison of (38) with (14) shows that the topological incidence matrices also appear in the finite element formulation and that the (weighted) mass matrices M (d−1) and M (d) once again play the role of the H-matrices which connect solutions on dual grids. In this section only the discretization on a single spectral element is discussed. Transformation of the domain [−1, 1] 2 to more general domains will be discussed in Section 4. The use of multiple elements follows the general assembly procedure from finite element methods. Results of this approach are presented in Section 5.

Transformation rules
The basis functions used in the discretization of the different physical field quantities have only been introduced for the reference domainΩ = [−1, 1] 2 . For these basis functions to be applicable in a different domain Ω , it is fundamental to discuss how they transform under a mapping Φ : (ξ , η) ∈Ω → (x, y) ∈ Ω ⊂ R 2 . Within a finite element formulation this is particularly useful because the basis functions in the reference domainΩ can then be transformed to each of the elements Ω e , given a mapping Φ e :Ω → Ω e .
Consider a smooth bijective map Φ : (ξ , η) ∈Ω → (x, y) ∈ Ω such that and the associated rank two Jacobian tensor J The transformation of a scalar function ϕ discretized by nodal values is given bỹ and of a scalar function ρ discretized by surface integrals is given bỹ The transformation of vector fields v v v discretized by line integrals is and of vector fields u u u discretized by flux integrals is These transformations affect only the mass matrices and not the incidence matrices. This is fundamental to ensure the topological nature of the incidence matrices.

Numerical results
In this section three test cases are presented to illustrate the accuracy of the discretization scheme developed in this work. The first test case, 5.1, is an analytical solution taken from [68] to assess the convergence rates of the method. The second test case, 5.2, is the flow through a system of sand and shale blocks with highly heterogeneous permeability in the domain, see for more details [54]. The third test case, 5.3, is a highly anisotropic and heterogeneous permeability tensor in the domain, see for more details, [53].

Manufactured solution
We first test the method using the exact solution with the permeability tensor given by The mixed formulation (3) in the form of (38) is then solved in the domain (x, y) ∈ Ω = [0, 1] 2 with the source term f = −∇ · (K∇p exact ) and the Dirichlet boundary condition p| ∂ Ω = 0. A benchmark of this test case for α = 0 using multiple numerical schemes can be found in [68]. When α = 0, K is multi-valued at the origin which makes this test case a challenging one. To see this, we can first convert the Cartesian coordinates (x, y) to polar coordinates (r, θ ) by x = r cos θ , y = r sin θ . Then we have K| α=0 = 10 −3 cos 2 θ + sin 2 θ 10 −3 − 1 cos θ sin θ 10 −3 − 1 cos θ sin θ cos 2 θ + 10 −3 sin 2 θ .
It can be seen that we get different K| α=0 when we approach the origin along different angles, θ . It must be noted that inverse of K does not exist at the origin. The inverse of the tensor term appears in 3.2.2. We use Gauss integration and thus the inverse term is not evaluated at the origin. The meshes we use here are obtained by deforming the GLL meshes in the reference domain (ξ , η) ∈ Ω ref = [−1, 1] 2 with the mapping, Φ, given as where c is the deformation coefficient. The two meshes, for c = 0.0 and c = 0.3, are shown in Figure 5. The method is tested for α ∈ {0, 0.01} and c ∈ {0, 0.3}. In Figure 6, the results for ||∇ · u u u h − f h || L 2 are presented. They show that the relation ∇ · u u u h = f h is conserved to machine precision even on a highly deformed and coarse mesh i.e. of 2 × 2 elements with N = 2 and c = 0.3.
When α = 0.01, K is no longer multi-valued at the origin. In this case the source term f is smooth over the domain, see Figure 7 (bottom). For this smooth case, the method displays optimal convergence rates on both the orthogonal mesh and the deformed mesh, i.e. see Figure 8 (bottom) and Figure 9 (bottom).
When α = 0, both the h-convergence rate and p-convergence rates are suboptimal, see Figure 8 (top) and Figure 9 (top). This is because K is multi-valued and therefore f becomes singular at the origin when α = 0, see Figure 7 (top left).

The Sand-Shale system
This example is taken from [54,76,78]. The domain is a 2D unit square, Ω = [0, 1] 2 , with 80 shale blocks, Ω s , placed in the domain such that the total area fraction of shale blocks is A shale = 20%, as shown in Figure 10. We solve the mixed formulation (38) with f = 0 in this domain. The flux across the top and the bottom boundaries is u u u · n n n = 0. The flow is pressure driven with the pressure at the left boundary, p = 1, and the pressure at the right boundary, p = 0. The permeability in the domain is defined as K = kI, where k is given by: For this test case an orthogonal uniform grid of 20 × 20 elements is used. The polynomial degree is varied to achieve convergence. Streamlines through the domain for 20 × 20 elements and polynomial degree N = 15 are shown in Figure 11. It can be seen that the streamlines do not pass through, but pass around the shale blocks of low permeability.
The ||∇ · u u u h || L 2 over the entire domain as a function of polynomial degree is shown in Figure 12. We observe that ∇ · u u u h = 0 is satisfied up to machine precision.
The net flux entering the domain (the same as the net flux leaving the domain) is given in Table 1 for varying polynomial degree. A reference value for this solution is given in [54] as 0.5205, and in [78] as 0.519269. In this work the maximum resolution corresponds to 20 × 20 elements and a polynomial degree N = 19, for which the net flux entering the domain is obtained as 0.52010. In Figure 13 we compare the net flux entering the sand-shale domain, calculated using the mixed and the direct formulation of equations, as a function of polynomial degree for different values of k in the shale blocks. The data for these figures is given in Table 2. Note that the direct formulation converges from above towards the correct inflow flux, whereas the mixed formulation converges from below.

The Impermeable-Streak system
The next example is from [53,76,78]. The physical domain is a 2D unit square, Ω = [0, 1] 2 . The domain is divided into three different regions, Ω 1 , Ω 2 , and Ω 3 , as shown in Figure 14 (left). For calculations, each region is further divided into K × K elements. Therefore, the total number of elements in the domain is given by K × K × 3. In Figure 14 (right) we show the domain with each region divided into 2 × 2 elements. The mixed formulation (38) is solved, with f = 0 and mixed boundary conditions, such that at the top and the bottom boundaries the net flux u u u · n n n = 0, and at the left and the right boundaries, p = 1 and p = 0, respectively. Permeability in Ω 1 and Ω 3 is given by K = I. Ω 2 has a low permeability and defined such that the component parallel to the local streak orientation is k = 10 −1 , and the component perpendicular to the local streak orientation is k ⊥ = 10 −3 . The analytical expression for the permeability in terms of Cartesian coordinates is given in [76] as, The flow field in the domain is shown in Figure 15. The magnitude of velocity in Ω 2 is small due to low values of the permeability tensor in this region. The velocity vectors bend in the direction of the permeability streak Ω 2 . The L 2 -norm of ∇·u u u over the entire domain as a function of polynomial degree, N, is shown in Figure 16. We can see that the flow field is divergence free up to machine precision because f = 0.  The net flux through the system for varying number of elements and polynomial degree is given in Table 3. In this work the finest resolution corresponds to 12×12× 3 elements and N = 15. For this case the net influx at the left boundary is 0.75668.
The net influx and outflux from the region Ω 1 , Ω 2 and Ω 3 is given in Tables 4, 5, and 6, respectively. The net influx for Ω 1 is larger than the net outflux. And the net outflux for Ω 2 and Ω 3 is larger than the net influx.  Table 6: Net flux through the left boundary of the region Ω 3 for K × K elements, K = 4, 6, 8, 10, 12 and N = 1, ..., 15.

Future Work
In the above sections, mixed and direct formulations of mimetic spectral element method are discussed. The next step is to explore this framework in the direction of hybrid formulations [20,35,43]. Additionally, the focus will be on developing multiscale methods [118], using these formulations, for reservoir modelling applications.