{'fraction_non_alphanumeric': 0.05714285714285714, 'fraction_numerical': 0.017217865666553017, 'mean_word_length': 3.790298873101421, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Crouzeix-Raviart Element for non-conforming dual mixed methods: A Priori Analysis

Under some regularity assumptions, we report an a priori error analysis of a dG scheme for the Poisson and Stokes flow problem in their dual mixed formulation. Both formulations satisfy a Babu\v{s}ka-Brezzi type condition within the space H(div) x L2. It is well known that the lowest order Crouzeix-Raviart element paired with piecewise constants satisfies such a condition on (broken) H1 x L2 spaces. In the present article, we use this pair. The continuity of the normal component is weakly imposed by penalizing jumps of the broken H(div) component. For the resulting methods, we prove well-posedness and convergence with constants independent of data and mesh size. We report error estimates in the methods natural norms and optimal local error estimates for the divergence error. In fact, our finite element solution shares for each triangle one DOF with the CR interpolant and the divergence is locally the best-approximation for any regularity. Numerical experiments support the findings and suggest that the other errors converge optimally even for the lowest regularity solutions and a crack-problem, as long as the crack is resolved by the mesh.


Introduction
For a simple read, we start with the Poisson's equation. Let Ω be a bounded and simply connected domain in R 2 with polygonal boundary Γ . Then, given f ∈ L 2 (Ω) and g ∈ H 1/2 (Γ) we look for u ∈ H 1 (Ω), such that (1.1) −∆u = f in Ω and u = g on Γ .
We follow [2] and introduce the gradient σ := −∇u in Ω as an additional unknown. In this way, (1.1) can be reformulated as the following problem inΩ : Find (σ, u) in appropriate spaces, such that (1.2) σ + ∇u = 0 in Ω , div σ = f in Ω and u = g on Γ .
In the classical approach [2] the scalar-valued unknown u is sought in broken H 1 and the vectorial unknown σ in L 2 . In the present work, we are interested in approximating σ in a discrete space that locally belongs to H(div) and the scalar-valued u ∈ L 2 . This allows to use piecewise constants to approximate u without slowing down convergence for σ, as it would in [2].
This kind of approach has been also applied in the previous works [3][4][5] and [7]. All of them consider the standard local Raviart-Thomas spaces, to approximate σ.
Alternatively, it is possible to use other pairs which weakly preserve the continuity of normal components of the vectorial unknown. This motivates us to approximate σ using the space that locally is the lowest order Crouzeix-Raviart (CR) element. This non-conforming element was introduced in the early work [13] in the framework of the Stokes problem analysed in the classical velocity pressure formulation. Since then, many advances have been made using this element, for an overview we refer to [9]. Sometimes, the CR-element is considered to be a special case of DG methods (on conforming meshes), where its use eliminates some terms of the stabilised bilinear form. However, the proof of the DG methods mentioned earlier (e.g. [6,7]) cannot be extended easily to this element, since the Raviart-Thomas space is not a subspace of the CR-space. Up to our best knowledge, there is no analysis developed considering this element for the Stokes problem formulated with the unknowns velocity and pseudostress, i.e., in its dual mixed formulation. However, a conforming scheme for this approach was introduced in [11].
Then, the proposal of this article is twofold: First, we explore the use of the CRelement in order to approximate the H(div) space for the mixed Poisson problem. After that, we extend the approach to approximate the solution of the Stokes equations in their velocity-pseudostress formulation. In both cases, we do require a conforming mesh and the discrete approximation of the divergence is the local best-approximation for any regularity. In numerical experiments we confirm that the approximation quality of the Raviart-Thomas space (for least regularity) can be improved within the CR space when weakly imposing the continuity of the normal component of σ.
The paper is organised as follows. We end this section introducing notation. Sections 2, 3 and 4 contain the design of the method, the existence and uniqueness of the discrete solution and stability, as well as, a priori estimates. In Section 5 we extend the applicability of the Crouzeix-Raviart element to Stokes system approximated by the unusual velocity-pseudostress formulation. Numerical examples are reported in Section 6.
In the rest of the paper we will use the following notation. Given any Hilbert space H, we denote by H 2 the space of vectors of length 2, and by H 2×2 the space of tensors, all with entries in H. Tensor notation will be used in Section 5. We also use the standard notation for Sobolev spaces and norms. In particular, let H(div; Ω) := τ ∈ [L 2 (Ω)] 2 : div(τ ) ∈ L 2 (Ω) . Finally, we use C or c, with or without subscripts, to denote generic constants, independent of the discretization parameters, which may take different values at different occurrences.

The new non-conforming formulation
In this section, we derive a discrete formulation for the linear model (1.1), applying an unusual discontinuous Galerkin method in divergence form. We begin with some definitions and notations.
2.1. Meshes. Let {T h } h>0 be a conforming family of triangulations ofΩ made up of straight-side triangles T with diameter h T and unit outward normal to ∂T denoted by ν T . As usual, the index h also denotes h := max T ∈T h h T . Then, given T h , its edges are defined as follows. An interior edge of T h is the (non-empty) interior of ∂T ∩ ∂T , where T and T are two adjacent elements of T h . We denote by E I the set of all interior edges of T h in Ω and by E Γ the list of all boundary edges, respectively. Then, E := E I ∪ E Γ denotes the skeleton inherited from the triangulation T h . Moreover, we denote by |e| the length of an edge e ∈ E and by |T | the area of T ∈ T h .

2.2.
Averages and jumps. Now, in order to define average and jump operators, let T and T be two adjacent elements of T h and x be an arbitrary point on the interior edge e = ∂T ∩ ∂T ∈ E I . In addition, let q and v be scalar-and vectorvalued functions, respectively, that are smooth inside each element T ∈ T h . We denote by (q T,e , v T,e ) the restriction of (q T , v T ) to e. Then, we define the averages at x ∈ e by: {q} := 1 2 q T,e +q T ,e and {v} := 1 2 v T,e +v T ,e . Similarly, the jumps at x ∈ e are given by [[q]] := q T,e ν T + q T ,e ν T and [[v]] := v T,e · ν T + v T ,e · ν T . On boundary edges e , we set {q} := q , {v} := v , as well as [[q]] := q ν and [[v]] := v · ν. These operators are well-defined for functions in the broken space H (T h ) with > 1/2. Since our approach is non-conforming, we also introduce the elementwise divergence operator div h and the broken Sobolev space H(div; T h ) defined in the standard way. Throughout, in order to shorten notation, we use the subspace 2.3. Discrete spaces. We define discrete spaces as follows. Let P (T ) the space of polynomials of degree on T ∈ T h . To shorten definitions we define the discontinuous space P (T h ) := v ∈ L 2 (Ω) : v| T ∈ P (T ) ∀T ∈ T h . Then, the scalarand vector-valued CR-spaces use the CR-element in each component as follows: Out method will approximate u within the space V h := P 0 (T h ) and σ in Σ h ⊂ Σ.

2.4.
The weak formulation of the Poisson problem. The global weak formulation of (1.2) is obtained as usual and reads: Find (σ, u) ∈ H(div; Ω) × L 2 (Ω) , such that In order to derive the discrete scheme, we consider the first equation in (1.2) on each T ∈ T h , multiply by a test function in H(div, T h ) and integrate by parts, to deduce Summing up the boundary terms we get where the signs denote the restrictions to the two neighbouring cells in T h . On each interior edge, we have the standard dG identity Finally, we realise that e [[τ h ]] {ũ} = 0 for all τ h ∈ Σ h and any approximatioñ u ∈ V h of u. Therefore, these terms will not appear in our discrete scheme.
2.5. The global discrete scheme. Considering (2.1b) and (2.2) we arrive at the following stabilised discrete dual mixed non-conforming Galerkin formulation: where the bilinear forms a s : where γ e = 1/ |e| on each edge. The linear functionals G : Σ → R and F : L 2 (Ω) → R are given by Remark 2.1. The jumps in the bilinear form a s will be used to control the consistency error visible when comparing (2.2) and (2.3a). The value chosen for γ e allows this error to decrease with the optimal speed (as in dG methods). Additionally, this value is later shown (Lemma 2.5) to not affect the existence of the discrete solution in a negative fashion. Even, without penalising the jumps, τ h can be thought of as weakly belonging to H(div, Ω), since e [[τ ]] = 0 for all e ∈ E I . However, the normal component of functions σ h ∈ Σ h jumps across edges and the penalty reduces this behaviour and the consistency error.
For arguments below we will use the following equivalent form of problem (2.3): In order to study the existence, uniqueness and approximation qualities of the discrete solution we introduce the following norms and properties of the CR interpolant. The space Σ is equipped with the norm · Σ : Σ → R, which is defined by where |·| a : Σ → R is defined by We note that the jumps vanish for τ ∈ Σ ∩ H(div, Ω). For V h we use the standard L 2 -norm. In addition, we define the natural norm (·, ·) nc : Σ × L 2 (Ω) → R of the method by 2.6. Properties of the Crouzeix-Raviart interpolant. Given ε > 1 2 , the CRspace and approximation operator Π cr : H ε (Ω) → CR h was first introduced in [13] and is defined as follows: ∀e ∈ E .
The vector-valued operator Π cr : [H ε (Ω)] 2 → Σ h uses the same definition for each component. A useful property of this operator is the following: We will also use the following approximation error estimate of the CR-interpolant.
) be a triangle (or tetrahedron). Then, the CR-interpolant satisfies the following estimate with a constant C independent of T : Proof. Estimate (2.10) is proved in Lemma A.2.
Furthermore, we will use the following inf-sup condition.
Proof. By [15] this result is equivalent to the inf-sup condition proved in [13].

2.7.
A modified inf-sup condition. For the existence and uniqueness of (σ h , u h ) we will show that τ in Lemma 2.3 satisfies τ Σ ≤ C τ 1,T h for appropriate parameters γ e . Lemma 2.5 proves this bound using the following lemma.
where C 1 depends on Lemma 2.4 and is independent of mesh and data. Now, a simplified version of [1, ,T with C 2 independent of T. Therefore, for γ| e satisfying the hypothesis the constant C sup ≤ (1 + 2C 1 C 2 )C Ω is independent of mesh properties which finishes the proof.

Existence and uniqueness
The well-posedness of Problem (2.3) and equivalently (2.4) is established in the next theorem. The proof is well-known and included for completeness.

Stability and a priori estimates
In this section, we focus on the stability and a priori error analysis for the scheme (2.3). The main advantage of our approach and using Σ h × V h is the inf-sup condition stated in Lemma 2.5 which allows a method without stabilisation terms for u h . This is similar to the conforming Raviart-Thomas pair. We will perform the analysis for parts of the norm (·, ·) nc . Hereafter, (σ, u) and (σ h , u h ) will be the unique solutions of (1.2) and (2.3), respectively.
Lemma 4.1. The consistency error is bounded as follows: The right term is part of |τ | a and the left term is bounded by C 1 h |u| 1,Ω because of γ −1 e = |e| and a standard trace estimate.
We now state stability estimates.
3) is stable in the following sense: Proof. Since div h σ h ∈ V h the second local estimate follows from (4.1b) using v := div h σ h . For the first bound, we consider the definition of B s (2.5) and the quasi-consistency (4.1) to deduce Therefore, from Theorem 3.1 we know there exists α > 0 , such that The last step above follows using Lemma 4.1 and the continuity of B s with respect to the norm (·, ·) nc , which is established by Cauchy's inequality.

Error estimates.
We continue with a priori estimates. The first one is a completely local best-approximation result.
T f s,T with a constant C > 0 is independent of mesh and data. In fact, if s = 1 then C = 1/π.
Proof. The consistency identity (4.1b) states |T | T div σ and then (div σ h )| T is the L 2 -projection into local constants on T. This is the best-approximation of div σ in V h which gives the first equality. The second infimum follows considering the div σ = f as stated in (1.2). Finally, applying Lemma A.1 and [8] to the second infimum gives the estimate.
where C > 0 is a constant independent of h, σ and u.
Proof. The following arguments are common and included for completeness.
As usual, we split the error in a discrete error and a projection error, i.e.
h + e u with e u h := u h − Πu and e u := Πu − u , where Π is the L 2 -projection into V h and Π σ will be chosen later. Then, the norm allows to bound the errors separately as follows: h , e u h ) nc + (e σ , e u ) nc . We first bound the discrete error (the left norm), later we concern ourselves with the projection error (the right norm). To this end, we reuse the consistency error identity (4.2). Hence, from Theorem 3.1 we know there exists α > 0 , such that Here we bounded the consistency error by Lemma 4.1. It remains to bound each part of the norm · nc defined in (2.7). We start using a Poincaré estimate (or the standard L 2 projection error, see [8]), to obtain The same Poincaré estimate allows to bound the divergence error as follows. We choose the Crouzeix-Raviart projection for Π σ . For this projection, by (2.9) we get the following equality and bound for all T ∈ T h . Moreover, by Lemma 2.2 and Lemma A.1 we get Furthermore, since [[σ]] = 0 pointwise we use Lemma A.3 to obtain Using the arguments above we get the final estimate Joining these estimates completes the proof.
Remark 4.5. It is known that the treatment of the Darcy flow is similar to Poisson's equation formulated in mixed form. In fact, the scheme presented here coincides with the scheme for the Darcy flow presented in [10]. Hence, the analysis exhibited here can be extended, in a natural way, to the porous media equations. Obviously, adapting the analysis to the boundary conditions considered in the model, which should be clear at least for Dirichlet, Neumann and Mixed type boundary conditions. In other words, the lowest order of Crouzeix-Raviart element for the velocity and piecewise constants for the pressure is an inf-sup stable but inconsistent pair for Darcy's law, too.

Stokes
In this section we extend the method to the incompressible Stokes flow for a velocity u : R 2 → R 2 and a pressure p : R 2 → R with p ∈ L 2 0 (Ω), i.e. div ν∇u − pI) = −f , div u = 0 in Ω and u = g on Γ with given ν > 0 and source terms f ∈ [L 2 (Ω)] 2 and g ∈ [H 1/2 (Γ)] 2 . Here and later on, I denotes the identity matrix and div gives a vector whose entries are the divergence of a row of the tensor. Furthermore, we use the following additional notation: the space H(div, Now, in order to use the dual mixed formulation we introduce the pseudostress σ := ν∇u − pI in Ω. Taking the trace and applying div u = 0 we realise that p = − 1 2 tr σ (independent of variations of ν ). This relation and p ∈ L 2 0 (Ω) give σ ∈ H 0 where H 0 := τ ∈ H(div, Ω) : Ω tr(τ ) = 0 ⊂ Σ 0 := τ ∈ Σ : Ω tr(τ ) = 0 . Additionally, defining the strong deviator τ d := τ − 1 2 (tr τ )I ∈ H 0 and using the previous identities we obtain σ d = σ + pI = ν∇u. Thus, we get the equivalent velocity-pseudostress formulation: Find (σ, u) ∈ H 0 × H 1 (Ω) 2 such that The pressure can be reconstructed from tr σ. A version of the previous discussion, as well as, the weak formulation of (5.1) can be found in [7]. The formulation reads: Furthermore, [7, Theorem 2.1] states that the solution (σ, u) of (5.2) is unique and stable with respect to the data, i.e., there exists a constant C > 0 independent of the solution such that 5.1. The non-conforming weak form and finite elements. Similar to Section 2.4, we approximate the tensor σ in each component by the CR h element and u by locally constant functions, that is, we consider the discrete spaces Σ h , Σ h,0 and V h defined by

Hence, integration by parts is done locally and gives
Using the same arguments as in Section 2.4 yields Using u ∈ H 1 (Ω) 2 , we obtain the weak form of the first identity in (5.1): This identity shows the consistency error which does not appear in the discrete formulation, because {u h } e ∈ R 2 and e [[τ h ]] = 0. However, this error has to be controlled by jumps just as in the other method. We will see, that the previous identities impose a significant simplification of the Galerkin scheme presented in [7].

5.2.
The theoretical and practical scheme. We now define define two numerical schemes. Both schemes approximate σ ∈ H 0 by σ h ∈ Σ h,0 , but the practical scheme imposes the average-free trace by a Lagrange multiplier and allows testfunctions with local support. The later proven equivalence of the schemes allows to reduce the theory. Identity (5.4) and the properties mentioned thereafter, give the following discrete (theoretical) scheme: where the bilinear forms a h : with the parameters γ| e = 1 |e| and G and F are defined as for (5.2).
The practical scheme reads: Given the property Σ h = Σ h,0 ⊕ RI, the following result relates the solutions of (5.5) and (5.6).
Second, we introduce a few identities to prove that every solution of (5.6) satisfies φ = 0. To this end, we write each τ ∈ Σ h as τ = τ 0 + ρI with τ 0 ∈ Σ h,0 and ρ = 1 2|Ω| Ω tr(τ ). Using this decomposition we conclude Furthermore, considering definitions we realise that These identities and (5.6a) with the test-pair (τ = I, v = 0) prove 2φ |Ω| = 0 and φ = 0. Furthermore, the previous identities show that (5.6a) is equivalent to Finally, taking a solution (σ h , u h ) of (5.5) the last equivalence shows that (σ h , u h , 0) solves (5.6a). Equation (5.6b) is satisfied by definition since σ h ∈ Σ h,0 . This finishes the proof. Lemma 5.1 proves that both schemes give the same solutions and, therefore, we can prove the existence and uniqueness for the scheme (5.5) only.

Existence and uniqueness.
We aim to prove an inf-sup condition for scheme (5.5) which confirms its existence and uniqueness. To this end, we define the following semi-norms and norms: Remark 5.2. It is not immediate that · Σ is a norm. Luckily, Lemma 3.10 in [7] proves the existence of a constant C > 0, independent of the meshsize, such that Hence, · Σ is a norm on Σ 0 Furthermore, since I Σ = 0 it is clear that this is not a norm on Σ.
Lemma 5.3. There exists a constant C sup > 0 independent of mesh and data, such that for every v ∈ V h , there existsτ ∈ Σ h,0 such that Proof. Forτ ∈ Σ h , the proof is analogous to the one of Lemma 2.5, since we may consider the separate independent cases v = (v 1 , 0) and v = (0, v 2 ). Then, considering div h τ = div h τ 0 and τ Σ = τ 0 Σ finishes the proof.

Stability, errors and convergence.
From now on (σ, u) and (σ h , u h ) will be the unique solutions of (5.1) and (5.5), respectively. As before, stability requires a bound on the consistency error. Therefore, we begin providing quasi-consistency identities. Taking the difference of (5.4) and (5.5a), as well as, (5.2b) and (5.5b) we get: Next we state a few estimates that are proven analogously to those for the Laplacian. The techniques used in Lemma 4.1 yield the following bound: The equation (5.12b) and the fact that div Therfore, just as in Lemma 4.3 we get Furthermore, (5.11), (5.12a), (5.12b) and (5.3) give the stability bound: Theorem 5.5. Let σ ∈ H 2 (Ω) 2×2 and u ∈ H 1 (Ω) 2 be the solution of (5.1) and (σ h , u h ) ∈ Σ h × V h be the solution of (5.5). Then, where C > 0 is a constant independent of data and h, σ and u.
Proof. The arguments are analoguous to Theorem 4.4, i.e., we split the error into a discrete error and a projection error: where Π is the L 2 -projection into V h and Π σ is the Crouzeix-Raviart projection. Then, the norm allows to bound the errors separately as follows: h , e u h ) nc + (e σ , e u ) nc . We bound the discrete error (the left norm), followed by the projection error (the right norm). The consistency error identity (5.4) and Theorem 5.4 yield the existence of α > 0, such that where we bounded the consistency error using (5.13). From here on, each part of the norm · nc defined just after (5.9) can be bounded as in Theorem 4.4 which ends the proof.
In summary we have shown that the same analysis works for the Poisson and the Stokes problem both formulated in a dual mixed non-conforming approach.

Numerical confirmation
In this section, we present numerical experiments that illustrate the performance of our method and confirm the convergence rates for smooth solutions and exhibit optimal convergence rates for non-smooth solutions that are not supported by our theory. The numerical experiments were performed with the finite element toolbox ALBERTA using refinement by recursive bisection [17]. The solutions of the corresponding linear systems were computed using the backslash operator of MATLAB.
In all cases, we start with an initial mesh T 0 , and in each step we solve the corresponding problem, compute the errors and refine uniformly to generate a conforming refinement T k+1 of T k by bisecting all elements twice. In what follows, DOFs k denotes the total number of degrees of freedom (unknowns) of the corresponding system on T k . Since we only use uniform meshes we have DOFs −r/2 ∼ h r . Hence, we calculate the empirical order of convergence (EOC), associated to some global error e k in step k, by EOC(e k ) := − 2 log(e k /e k−1 ) log(DOFs k /DOFs k−1 ) , k > 0.
In order to impose the continuity of the normal trace across inter-element edges, we always choose the jump parameter γ| e = |e| −1 . We present three numerical examples for each of the analyzed models (Poisson and Stokes). The first one exhibits a smooth behavior, completely supported by our theory. But, since the divergence error equals the best-approximation error for any regularity, we also analyze the performance of the methods for a singular solution on an M-shaped domain and on a crack domain (opening angle zero) where the crack is resolved by the mesh and where σ / ∈ H s (Ω), s ≥ 1/2. The used opening angle of the crack is zero, but Figure 1 shows an open crack to illustrate the "additional" DOFs on the "overlapping edges". crack Figure 1. Initial meshes for problems on M-shaped domains and the crack problem (real opening angle is zero).
Note that asymptotically 33% of the DOFs for σ h are fixed by the best approximation of div σ. This is quite a big percentage, so we analysed the perfomance of the method for the non-smooth cases to see whether the other DOFs are also close to an optimal converging approximation, as the smooth convergence result proves. It turns out, that for all cases we obtain optimal convergence (a result left open to be proven as part of future work).

Poisson problems.
6.1.1. P1. A smooth test solution. In this first test, we consider a smooth solution on an M-shaped domain with Dirichlet boundary conditions given by: The objective of this example is to confirm that the analyzed scheme provides the optimal rate of convergence. The initial mesh on the non-convex domain Ω is shown in Figure 1 (left). Additionally, we report the convergence rates in Table 1 which are in agreement with the theory. Finally, in Figure 2 we present the module of σ h and u h at iterations 2, 4 and 6, respectively. The pictures seem to indicate that our method provides a field σ h with a continuous normal trace across inter-element edges.    Figure 1 (left), and consider the following solution of (1.1): where (r, θ) denote polar coordinates. Note that u is the solution of the elliptic equation −∆u = 1 in Ω, and that u ∈ H 5/3 (Ω) (we ignore the −ε ). Figure 3 shows the discrete solution u h (piecewise constant) and the module of σ h (module of a CR field) at iteration 6, respectively. We appreciate some small oscillations in the second picture, due to the singularity, but this does not affect the rates reported in Table 2. We omit the div h (σ − σ h ) error, because it is almost zero ( f = 1 is constant). Since σ = ∇u ∈ H 2/3 (Ω) at the reentrant corner, we expect (and see) the rate h 2/3 ∼ DOFs −1/3 for the L 2 -error of σ. The jump terms converge faster and the L 2 -error of u converges with the optimal rate h ∼ DOFs −1/2 .  Figure 1 (right) with opening angle zero, and consider the following singular solution of (1.1): where (r, θ) denote polar coordinates. Note that u is the solution of the elliptic equation ∆u = 0 in Ω. The initial mesh is shown in Figure 1 (right). We emphasize that the DOFs on the line (0, 1) × {0} are different. The error behaviours for this test are reported in Table 3. We observe optimal rates for the error of variable u, and the expected rate for the error of σ (according to the singularity). Finally, in Figure 4 we show the discrete solution u h and the module of σ h at iteration 5.  6.2. Stokes problem. In order to impose the zero mean value condition for the trace of functions in the spaces Σ h,0 , we consider the equivalent problem (5.6), that was introduced in Section 5.2. We recall that from definition of pseudostress σ, the pressure can be recovered as: p = − 1 2 tr(σ). This identity allows us to analyze the pressure error, since 6.2.1. S1. Smooth solution. This example is motivated by the two-dimensional analytical solution of the Navier-Stokes equations derived by Kovasznay in [16], where we consider Ω = − 1 2 , 3 2 × (0, 2) and: with the parameter λ given by: In order to test the robustness of our scheme, we solve the problem for different values of the viscosity, ranging from 1 to 10 −5 . Figure 5 shows the decay of the errors versus the DOFs. Convergence rates predicted by the theory (Theorem 5.5) are attained in all cases, independently of the value of the viscosity. We remark that the pressure decay with rate h 2 = DOFs −1 .   Figure 5. S1. Decay of errors for ν ∈ {1, 10 −1 , . . . , 10 −5 }.
In Figure 6, we show the velocity stream lines and pressure at iteration 4, for the viscosity values analyzed. We do not appreciate any instability featured caused by a small viscosity.  Figure 1(left), and take f ≡ 0 and ν = 1. Then, if (r, θ) denote the polar coordinates, we impose an appropriate boundary condition for u so that: Here, ω = 3π 2 and the coefficient λ is the smallest positive solution of: sin(λω) + λ sin(ω) = 0 ⇒ λ ≈ 0.54448373 . . .
We emphasize that (u, p) are singular functions, and at the origin u ∈ H 1+λ and p, σ ∈ H λ . Table 4 gives the individual errors and the corresponding rates. Again, we omit the divergence error because its size is close to the rounding unit.
As in the Poisson problem, we note that the jump and the L 2 -error of u decay with rate h 1 . The L 2 error of σ decay according to the regularity of the exact solution. Figure 7 show the velocity and pressure at iteration 3 and 6. We observe some oscillations in the first steps that diminish later as the mesh is refined.  Figure 1(right) with (u, p) given by (6.5) and (6.7) Ψ(θ) = 3 sin(θ/2) − sin(3θ/2), λ = 1/2, ω = 2π . Table 5 reports the error decay. Again, we leave out the divergence error as it is close to the rounding unit. Although, this example is not covered by our theory, we obtain results equivalent to those in Section 6.1.3. We remark that the observed rates for the L 2 errors in σ and p are slightly better than the expected value (0.5). Finally, Figure 8 shows the velocity-field and pressure at iteration 2 and 5.

Final comments and conclusions
We developed an a priori error analysis for a low-order non-conforming method applied to the dual-mixed variational formulation of the Stokes and Poisson problem, and thus also for Darcy flow. Our a priori error estimates were established using Crouzeix-Raviart elements to approximate the σ = ∇u and piecewise constants to approximate u. We proved a priori estimates for smooth solutions. Low regularity and a posteriori estimates are part of a forthcoming paper. However, the divergence error was proven to converge optimally for any regularity. As a by-product we proved that the L 2 -error of the CR-interpolant (if well-defined) is equivalent to the best-approximation error in a discontinuous P 1 -space. The approximation errors of our non-conforming scheme improve those of the conforming Raviart-Thomas approach.
Finally, an advantage of this approach over the one with Raviart-Thomas elements, is that, the CR-element is defined component-wise and not vector-wise. In fact, this leaves more flexibility for a symmetry restriction when dealing with elasticity equations, which is much more difficult to achieve with the Raviart-Thomas element. In particular for the lowest order Raviart-Thomas pair, this could only be done for a tensor with constant entries. We therefore, intent to extend the CR-element to the elasticity problem in a separate paper. Recalling γ = |e| −1 ≈ h −1 T and joining the estimates and identities we arrive at the result.