The Finite Volume Formulation for 2D Second-Order Elliptic Problems with Discontinuous Diffusion/Dispersion Coefficients

We propose a finite volume method for the numerical resolution of two-dimensional steady diffusion problems with possibly discontinuous coefficients on unstructured polygonal meshes. Our numerical method is cellcentered, secondorder accurate on smooth solutions and based on a special numerical treatment of the diffusion/dispersion coefficients that makes its application possible also when such coefficients are discontinuous. Numerical experiments confirm the convergence of the numerical approximation and show a good behavior on a set of benchmark problems in two space dimensions.


Introduction
In the last decades, finite volume methods have been greatly successful in solving engineering models of flows in porous media on complex geometrical domain because the finite volume formulation works on general polygonal and polyhedral meshes.This great mesh flexibility is now combined with a strong theoretical foundation, that is, convergence analysis and error estimates are available on the simplest mathematical models 1 .In this work, we address a version of the finite volume method that is popularly known as "the diamond scheme" and was originally presented for the advection-diffusion equation in two dimensions in 2-4 and successively extended to convection-dominated problems in 5, 6 and nonlinear flow problems in partially saturated porous media 7 .
The numerical treatment of both diffusion/dispersion flux and of the partially saturated flow is based on averaging one-sided gradients and diffusion conductivity tensors.Many papers addressed the issue of conductivity averaging both in finite differences and in finite elements.A seminal work 8 started with a finite difference application to the unsaturated flow equation.Later on they, dealt with the problem in a finite volume scheme 9 .Recent reviews on this topic can be found in 10-12 .They also made performance comparisons using different averages.The latter paper shows in conclusion that their scheme, namely, an averaging schemes based on Darcian mean principle used in the framework of either vertex-centered or cell-centered approach compare favorably to other methods for a range of test cases.They compared their method with the one of 11 .Aside of the numerical comparisons, in 13 it is remarked that the usage of the harmonic average can be motivated on the basis of physical arguments.
Another reason for dealing with the tensor averaging issue is related to the scale problems.Typically, the Darcy equation is written at a scale much smaller than the scale of practical interest 14, 15 .This gives boost to new numerical techniques to able to deal with the heterogeneity of hydrodynamic parameters.One possibility is to obtain either analytically or empirically explicit equations for the scale of interest, eliminating other scales in the problem 16 .
Several methods have been recently developed, such as, the heterogeneous multiscale method HMM 16 and the multiscale finite volume element method 17 .Babuška in the 70s motivated the multiscale finite element method 18-20 .Multiscale methods have been proposed also in the case of saturated flows in heterogeneous porous media and also applied the multiscale finite volume element method to the Richards equation 21, 22 .In 14 , a multiscale method based on the framework of HMM was recently extended to solve the Richards equation with random coefficients.Another feature of their numerical scheme is the formula estimating the macroscopic flux in which the unsaturated hydraulic conductivity can be calculated as a diagonal tensor.In particular, it is worth noting that a finite volume high-order accurate approximation of the pressure head field also allows one to achieve a better resolution in the approximation of other important fields like the components of the hydraulic conductivity tensor at the mesh vertices 7 .For the sake of simplicity, in this work it was considered the steady flow in a layered porous medium in the presence of a source term.It is the Poisson equation which is suitable to test new schemes which will be applied in further work to the more general unsaturated transient cases.
We propose a novel technique that automatically adjusts when the diffusion tensor is discontinuous across a mesh interface shared by two adjacent cells.This technique is general and can be easily implemented in any finite volume scheme that has an explicit numerical flux and may result in a particularly efficient strategy for the numerical resolution of problems where both the diffusive/dispersive and convective phenomena are simultaneously significant.In fact, in such problems high-resolution finite volume methods are normally coupled with mixed finite elements following the criterion of choosing the best available technique in accordance with the nature of the equation to be discretized.In fact, the RT 0 − P 0 method is more suitable to the numerical discretization of the diffusive part of the model, and the finite volume method gives an accurate and stable discretization of the convective part of the model, even in presence of strong convective fields.This approach was proved successful in the numerical modeling of oil reservoir problems 23 and of groundwater flow and transport of contaminants in porous media, combared with 24-26 .Moreover, other different engineering areas may benefit from this new technique, such as 27-38 .
The outline of the paper is as follows.In Section 2, we introduce the mathematical model and discuss the formulation of the finite volume methods based on vertex reconstruction.In Section 3, we present the numerical treatment of the diffusion tensor.In Section 4, we confirm the theoretical results with numerical experiments.In Section 5 we offer final remarks and conclusions.

The Mathematical Model and the Finite Volume Formulation
Let Ω be a polygonal domain with boundary Γ.
We consider the steady diffusion problem for the scalar solution field u given by where Λ is the diffusion tensor describing the material properties, f is the forcing term and g the boundary function that defines nonhomogeneous Dirichlet conditions on the boundary Γ.Under suitable regularity assumptions on f, g, Λ, and Ω, it turns out that the diffusion problem is mathematically wellposed, a unique solution exists 39 and such a solution is continuously dependent on the model data.In particular, due to the elliptic nature of the model, the diffusion tensor Λ x is normally represented by a strictly positive definite matrix for every x in Ω.The components of the diffusion tensor may be discontinuous in the computational domain; in such a case and without loss of generality, we assume that the mesh is conforming with each discontinuity of Λ.
The numerical approximation to 2.1 -2.2 is performed on a sequence of polygonal partitions {Ω h } h of the domain Ω.For any mesh Ω h , the subscripted label h is the mesh size and is defined by where h K is the diameter of the polygonal cell K ∈ Ω h .We also label the numerical solution u h calculated on the mesh Ω h by h.Also, i a generic mesh vertex is denoted by V and its coordinate vector by x V .We will also find it convenient to introduce a local numbering of the vertices, for example, V i , and to ease notation to denote the vertex position of the i-th vertex by x i instead of ii a generic cell interface or a boundary edge is denoted by σ, its center i.e., its edge midpoint by x σ , and its measure edge length by |σ|; iii a generic polygonal cell is denoted by K its measure area by |K|, its center of gravity by x K , and its boundary by ∂K.
The orientation of each mesh interface σ is reflected by its unit normal vector n σ , which is fixed once and for all.For any mesh cell K and any face σ of the polygonal boundary ∂K, we define the unit normal vector n Kσ that points out of K and we also use the notation N Kσ |σ|n Kσ .
Let the algebraic vector u h u K K∈Ω h be the numerical solution, where each u K approximates the cell average of the scalar solution u over the cell K.The finite volume scheme for u h on the computational mesh Ω h reads where Λ σ is an evaluation of Λ at face σ taken from the side of K; ∇ σ is the discrete face gradient built inside the diamond cell centered at face σ; N K,σ is the geometric vector perpendicular to σ, pointing from K to L, and with lenght |N Kσ | |σ|; f K is the cell average of the right-hand side term f over the cell K: Let Λ K and Λ L be first-order approximations of the diffusion tensor Λ within the control volumes K and L. For example, we can take Λ K : Λ x K and Λ L : Λ x L , where x K and x L are internal points not necessarily the barycenters of K and L, respectively.Now, we can consider either Λ σ Λ K inside K and Λ σ Λ L inside L or a suitable mean of these two tensors.In particular, we can deal with the arithmetic mean: or with the harmonic mean: Both approaches deserve a special care when Λ is discontinuous across σ and a proper definition of the discrete gradient ∇ σ u h is also needed in order to preserve the property of flux conservation.
To define the numerical gradient, a special control volume is built around this interface, which has a quadrilateral shape in two dimensions and is named "diamond cell".The geometry of a diamond cell is shown in Figure 1, which plots the diamond cell D centered at the mesh face σ with vertices x i and x i 1 and shared by the cells with centers x K and x L .The diamond cell D can also be seen as the union of the subdiamonds D K and D L , which are the triangular cells sharing σ as common base and having, respectively, vertices x K and x L , the centers of gravity of cells K and L. The four vectors N Ki , N Li , N Li 1 , and N Ki 1 shown in Figure 1 are respectively orthogonal to the four boundary faces σ Ki , σ Li , σ Li 1 , and σ Ki 1 , and their lenght is equal to the lenght of the corresponding face.Instead, when σ is a boundary face, thus defined by σ ∂K ∩ Γ Γ being the boundary of the computational domain Ω , the diamond cell associated to σ coincides with D K , that is, it is the triangle defined by σ and the vertex x K .
The numerical diffusive flux is built by using a constant approximation of the solution gradient and an evaluation of the diffusion tensor Λ within the diamond cell, as discussed before.Let us give the formulas for an internal mesh face σ, that we suppose to be shared by two cells K and L. All these formulas can be easily adapted to the case of a boundary face by taking x L equal to the center of gravity of face σ.
Using the Green-Gauss theorem yields: where n is the unit vector orthogonal to σ ∈ ∂D and pointing out of D. If the restriction of u to the face σ of ∂D is an affine function, the boundary integral on ∂D in 2.8 only depends on the values of u at the vertexes of D and the constant vector provided by the formulas that we derive below must coincide with the gradient of u.
The Gauss-Green theorem applied to the solution gradient ∇u integrated on the 2D diamond cell D of Figure 1 yields:

2.9
Let us introduce the vectors N KL and N ii 1 , which are, respectively, orthogonal to the edge connecting x K to x L and x i to x i 1 , and whose lenght are equal to the length of these edges.
Using the geometric relations N KL − N Ki N Li and N KL − N Ki N Li into 2.9 and rearranging the terms yields:

2.10
To derive the gradient formula, we replace the function values u x K and u x L with the cell unknowns u K and u L , respectively, and the function values u x i and u x i 1 with the corresponding nodal unknowns u i and u i 1 .We obtain which is a constant approximation of the solution gradient ∇u within the quadrilateral cell D.
With any vertex x i of the mesh Ω h , we associate the reconstructed vertex value u i K∈Ω i ω Ki u K , where Ω i {K : x i ∈ K} denotes the subset of the mesh cells which share the vertex x i .The interpolation weights ω Ki are assumed to verify the consistency relations 3 : Alternative constructions of the interpolation weights are available in the literature, see, for instance, 40 .The main advantage offered by these alternative choices concerns the control of positivity of the weights, which turns out to be significant when we are aimed at developing a numerical method with a maximum or minimum principle.However, as this topic is out of the scope of the current work, we will not pursue it anymore.The interpolation weights ω Ki are obtained by solving the reconstruction problem that approximates the cell-averaged data set { x K , u K K∈Ω i } by the affine function: on the covolume V i K∈Ω i K and in a least square sense, compared to 3, 41 .The reconstructed value at vertex x i is now given by taking u i u i x i α.The coefficients α, β are the minimizers of the least squares functional:

2.14
Imposing the zero gradient condition, that is, ∇ α,β J α, β 0, yields a linear system for the coefficients α, β , whose solution returns the interpolation weights.For completeness of exposition, we briefly review the derivation of the weight formula given in 41 .Let β β x , β y ; a straightforward calculation yields:

2.15
Let m i be the number of cells of Ω i ; we reformulate the linear system written above as 17 For convenience, we consider a suitable local numbering of the cells of V i , for example, , and we introduce the two matrices: and the column vector u u 1 , u 2 , . . ., u m i T , which collects the solution averages to be imposed.Using such definitions, the linear system 2.16 -2.18 becomes since matrix A is, obviously, symmetric and positive definite and, thus, nonsingular.The coefficient α that provides the vertex value is given by

2.21
We require that such a coefficient be an average of the data in u with coefficients ω ω K 1 i , ω K 2 i , . . ., ω K m i i T , that is, α ω T u.By comparison with 2.21 , we have the final weight formula: Conditions 2.12 can be checked by a direct calculations.The functional in 2.14 can be suitably modified to take into account different kind of boundary conditions, for example, Neumann or Robin, compared to 41 .The values u i at the vertexes x i ∈ Γ on the Dirichlet boundary are constrained to the boundary data, for instance u i 0 for a homogeneous condition.We also mention that this choice of weights provides a very robust technique, compared to the study of locking effect due to the accuracy of the reconstruction 42 .

Treatment of Discontinuous Diffusion Tensors
We discuss, here, how to treat the case of a diffusion tensor that is discontinuous across an interface σ shared by the cells K and L, that is, σ K | L. Since we suppose that Λ is discontinuous across σ and the normal flux of the exact solution is continuous across σ, the normal component of the solution gradient must be discontinuous.Consequently, a numerical approximation of the diffusive flux across σ like cannot be consistent whenever Λ Kσ is some kind of average between Λ K and Λ L and the constant vector ∇ σ u h approximates ∇ u over all the diamond cell hull x K , x L , σ .To tackle this problem, we consider two distinct approximations of the solution gradient within the subcells D K and D L , denoted by ∇ Kσ u h and ∇ Lσ u h , respectively, and impose directly the flux conservation.To derive an expression for such one-sided gradients, we introduce an additional unknown u σ along face σ and we apply the Gauss-Green Theorem as for the derivation of ∇ KL u h in the previous section.Also, we recall that N Kσ N KL and |N Kσ | |σ| and we use a similar notation for the normal vector from the other side of σ, so that N Lσ −N KL and |N Lσ | |σ|.The two one-sided gradient formulas are

3.2
Now, we search for a tensor Λ Kσ , which is an average of the diffusion tensors Λ k and Λ L , and a gradient vector ∇ KL , which is an average of the gradient vectors ∇ Kσ u h and ∇ Lσ u h , that makes it possible to express directly the flux conservation as 3.3 Since we expect that the normal component of the vectors ∇ Kσ u h and ∇ Lσ u h be the same, we impose directly this condition by requiring that a scalar coefficient ϕ exists such that recall that n σ n KL .We also assume that ∇ KL u h be a weighted average of ∇ Kσ u h and ∇ Lσ u h .Let us be given two nonnegative coefficients μ K and μ L such that and such that To this purpose, a number of possible choices have been proposed and widely investigated in the literature for μ K and μ L .For example, a popular option is to take the volume of the two subdiamonds D K and D L normalized by the total volume of the diamond, for example, The diffusion tensor can be chosen as the simple average of Λ K and Λ L , for example, Λ KL μ K Λ K μ L Λ L , or as the harmonic average, for example, L .Much more attention must be paid when Λ is discontinuous across σ as none of the previous choices may provide the correct flux information across the common interface of cells K and L. Our aim in the next developments is at determining a value for ϕ and of the average diffusion tensor Λ KL in terms of Λ K and Λ L that ensures flux conservation 3.3 for any given pair of coefficients μ K and μ L .Our derivation is similar to that considered in 43 , and as pointed out therein, the resulting matrix is the same obtained for other purposes in the field of upscaling of conductivity, either by means of the large average techniques 44 or using the homogenization theory in the case of layered materials 45 .
To this purpose, we first derive an expression for ϕ.We multiply 3.4 by Λ L n σ and use flux conservation to obtain from which we have Likewise, we multiply 3.4 by Λ K n σ and use flux conservation to obtain from which we have

Mathematical Problems in Engineering
We multiply 3.8 by μ K and 3.10 by μ L , we sum the resulting relations and use 3.5 to have

3.11
Solving this last equation for ϕ gives us the formula: Now, we derive an expression for Λ KL .In view of 3.5 and 3.3 , it holds that First, in 3.13 , we substitute the expression for ∇ Lσ u h provided by 3.4 , we collect the factor ∇ Kσ u h and we obtain: Second, in 3.13 , we substitute the expression for ∇ Kσ u h provided by 3.4 , we collect the factor ∇ Lσ u h and we obtain

3.15
Third, we multiply 3.14 by μ L , 3.15 by μ K , we add the two resulting equations, we use again 3.5 and 3.6 , and we obtain Finally, in 3.16 we substitute the scalar factor ϕ given by 3.12 , we collect the vectors n F and ∇ KL u h and we obtain: 3.17 Equation 3.17 suggests us to set where

Numerical Experiments
In this section, we present and discuss a set of numerical results to show the convergence behavior of the finite volume method when we compute the diffusive flux using the technique described in the previous sections.In all the test case that we present in this section, we compare the behavior of the new numerical treatment that we consider in this paper with the weighted average, which is the standard approach in the method 5, 41, 42 .In fact, from the formulation given in equation 3.18 , it follows that the average diffusion tensor λ KL is equal to the weighted average μ K Λ K μ L Λ L plus a correction term.This correction term is specifically designed to take care of possible discontinuities in the diffusion coefficients.
The finite volume formulation based on the least squares reconstruction of vertex values leads to a linear system for the cell-average unknowns whose coefficients matrix is generally unsymmetric although displaying a symmetric nonzero pattern.The positive definiteness of such system is still an open issue and this fact is also the major difficulty for the development of a full theory of convergence of such method.Theoretical results that prove coercivity are available only for meshes of slightly deformed quadrilaterals 2-4 .
We solve such linear system using the MA41 routine of the HSL software collection, which implements an unsymmetric multifrontal sparse LU factorization technique especially designed for matrices with a symmetric nonzero pattern and unsymmetric values.Different software packages like UMFPACK and standard Krylov methods like BiCG-Stab and GMRES can be used in alternative.Efficiency issue is beyond the scope of our investigation but more details and comparison of performance when implementing different linear algebra solvers are found in the benchmark of the FVCA-6 Conference held in Prague, Czech Republic, in June 2011, 46 .
We solve 2.1 -2.2 on the domain Ω 0, 1 × 0, 1 for the data specified in the three test cases reported below.For all calculations, we measure the following relative errors: error on the solution: error on the gradient:

4.2
Test Case 1 This test case is taken from 47 and is devoted to confirm that the new treatment of tensor coefficients and the standard weighted arithmetic average show the same behavior when the diffusion coefficients are regular, for example, constant or smoothly varying functions of position, even if small anisotropies along the principal direction of diffusion are present.The exact solution that we want to approximate is given by u x, y sin 2πx sin 2πy x 3 xy 2 .

4.3
We consider two different constant diffusion tensor.The first one, called Λ iso , is isotropic, while the second one, called Λ ani , is anisotropic.The two diffusion tensors are given by: We solve this test case on mesh family M 1 .Each mesh of such family is built as follows.
First, we remap the position x, y of the nodes of an n × n uniform partition by the smooth coordinate transformation: x x 1 10 sin 2π x sin 2π y , y y 1 10 sin 2π x sin 2π y .

4.5
The corresponding mesh of M 1 is built from this "primal" mesh by splitting each quadrilateral cell into two triangles and connecting the barycenters of adjacent triangular cells by a straight segment.The mesh construction is completed at the boundary Γ by connecting the barycenters of the triangular cells close to Γ to the midpoints of the boundary edges and these latters to the boundary vertices of the "primal" mesh.The first and the second mesh of mesh family M 1 are displayed in Figure 2; mesh data are given in Table 1.Approximation errors for the isotropic diffusion tensor Λ iso are shown in Figure 6; approximation errors for the anisotropic tensor Λ ani are shown in Figure 7.In both plots, we also show the exact slopes proportional to h and h 2 .From Figures 6 and 7, we deduce that in the case of constant diffusion tensors the performance of the new diffusion average proposed in this work and the weighted arithmetic average are almost the same.This behavior is confirmed both in the case of an isotropic diffusion tensor and in the case of an anisotropic diffusion tensor.

Test Case 2
In this test case, we show the behavior of the new technique that is proposed in this paper when the diffusion coefficients are discontinuous along an internal line of the computational domain.To such purpose, we split Ω ∪ 2 i 1 Ω i where

4.6
The diffusion tensor is discontinuous across the horizontal line y 1/2 and is given by where and λ 1 1, λ 2 10.The exact solution is continuous across the line y 1/2 and is designed to ensure flux conservation, that is, continuity of the normal component of the flux field.The exact solution is given by This test case is solved on the two mesh families described below.

Mesh Family M 2 Mathematical Problems in Engineering
Mesh Family M 3 Below y 1/2, we consider a regular mesh formed by 2n × n squares, while above y 1/2 we consider a regular mesh formed by 4n × 2n squares.All the mesh nodes except those located at the boundaries and those located at the internal discontinuity line y 1/2 are then displaced by perturbing their position to a random position inside a square box centered at that original node position.The sides of this square box are aligned with the coordinate axis and their length is equal to 0.8/n note that 1/n is the distance between two consecutive nodes in each direction .Instead, the nodes lying at y 1/2 are allowed to change only the abscissa in order not to modify the shape of the interface line while the nodes at the boundary are not displaced.We treat a nonmatching mesh as a conformal polygonal mesh modifying the shape of the polygons that are immediately below the line y 1/2: these polygons are treated as degenerate pentagons with two parallel consecutive edges.The first mesh is built by taking n 10 and mesh refinement is implemented by doubling the value of n at each refinement step, thus implying that the mesh size parameter is approximately halved when passing from one mesh to the next one in the refinement process.The first mesh and the first refined mesh of this mesh suite are shown in Figure 4; mesh data are reported in Table 3.
The numerical results for the gradient and the solution approximation on mesh family M 2 are shown in the two plots of Figure 8 and on mesh family M 3 in the two plots of Figure 9.In the bottom-left corner of both plots, we show the exact slopes proportional to h 1/2 gradient errors and to h and h 2 solution errors .The solution gradient has a discontinuous normal component across such lines, and, due to this lack of regularity, the convergence rate cannot be expected to be better than that of a first-order scheme.This behavior is reflected in both left plots of Figures 8 and 9.Even if the convergence rate of the gradient is the same in the two cases, the approximation errors of the formulation using the new average are smaller than those obtained when using the standard weighted arithmetic average.Concerning the solution approximation, the results are even more spectacular because adopting the new average allows us to recover the second-order convergence rate, thus confirming the superior behavior of the new method.

Test Case 3
In this test case we aim at confirming the behavior of test case 2 for a more difficult case in which the diffusion tensor has intersecting discontinuities with an internal cross point.To such purpose, we split the computational domain as Ω ∪ 4 i 1 Ω i where

Mathematical Problems in Engineering
The diffusion tensor is discontinuous across the horizontal line y 1/2 and the vertical line x 1/2 and is given by

4.12
This test case is solved on the mesh family M 4 5 that is built by displacing randomly all the nodes except those at the subdomain interfaces x 1/2 and y 1/2 of a n × n regular partition of Ω.We consider five meshes with n 10, 20, 40, 80, 160; mesh data are shown in Table 4, and the first mesh and the first refined mesh of this mesh suite are shown in Figure 5.
The numerical results for the gradient and the solution approximation on mesh family M 4 are shown in the two plots of Figure 10.In the bottom-left corner of these plots, we show the exact slopes proportional to h 1/2 gradient errors and to h and h 2 solution errors .As in test case 2, the normal component of the solution gradient is discontinuous across the internal lines x 1/2 and y 1/2.Therefore, the convergence rate that is numerically measured in this test case is proportional to h 1/2 , cf. the left plot of Figure 10.This fact is independent of the average technique that is used for the diffusion coefficients.Nonetheless, the approximation errors given by the new average are smaller than those obtained when using the standard weighted arithmetic average.Concerning the solution approximation, such results confirms the behavior seen in test case 2. Consequently, the finite volume approximation of the scalar solution recovers the optimal second-order convergence rate when we apply the new average technique.

Conclusions
The numerical treatment of diffusion coefficients is an open problem and major problem, as, for example, the conductivity of different layers of soil is abruptly spatially variable.Moreover, it is particularly important in view of the problems of different scales involved and of random fields of the coefficients themselves.To this purpose, we proposed and tested a new technique for the numerical treatment of the conductivity tensor that interpolates the information coming from the different sides of a mesh interface shared by two adjacent cells.We point it out that this design is particularly suitable to the case of discontinuous conductivity tensors since the interpolation automatically adjusts its value by introducing an appropriate correction to the standard weighted average.Numerical experiments confirm this behavior.Furthermore, it may result in a particularly efficient strategy for the numerical resolution of problems where both the diffusive/dispersive and convective phenomena are simultaneously significant.

Figure 2 :
Figure 2: First and second mesh of mesh family M 1 .This mesh family is used in test case 1.

Figure 3 :Figure 4 :Figure 5 : 2 bFigure 6 :
Figure 3: First and second mesh of mesh family M 2 .This mesh family is used in test case 2.

Table 1 :
Data of mesh family M 1 ; l is the mesh label, N K is the number of cells, N σ is the number of mesh edges, N V is the number of mesh nodes, and h is the mesh size parameter.

Table 2 :
Data of mesh family M 2 ; l is the mesh label, N K is the number of cells, N σ is the number of mesh edges, N V is the number of mesh nodes, and h is the mesh size parameter.

Table 3 :
Data of mesh family M 3 ; l is the mesh label, N K is the number of cells, N σ is the number of mesh edges, N V is the number of mesh nodes, and h is the mesh size parameter.

Table 4 :
Data of mesh family M 4 ; l is the mesh label, N K is the number of cells, N σ is the number of mesh edges, N V is the number of mesh nodes, and h is the mesh size parameter.
α 3 1 and α 2 α 4 100.The exact solution is continuous across such line and is designed to ensure flux conservation, that is, continuity of the normal component of the flux field.The exact solution is given by i , i 1, . . ., 4.