The PCD Method on Composite Grid

abstract: We introduce a discretization method of boundary value problems (BVP) in the case of the two dimensional diffusion equation on a rectangular mesh with refined zones. The method consists in representing the unknown distribution and its derivatives by piecewise constant distributions (PCD) on distinct meshes together with an appropriate approximate variational formulation of the exact BVP on this piecewise constant distributions space. This method, named the PCD method, has the advantage of producing the most compact possible discrete stencil. Here, we analyze and prove the convergence of the PCD method and determine upper bounds on its convergence rate.


Introduction
We propose a discretization technique of boundary value problems (BVP) in which the unknown distribution and its derivatives are all represented by piecewise constant distributions (PCD) but on distinct meshes.The only difficulty of the method lies in the appropriate choice of these meshes.Once done, it becomes 2000 Mathematics Subject Classification: 65N12, 65N15, 65N50

122
A. Tahiri rather straightforward to introduce an appropriate approximate variational formulation of the exact BVP on this piecewise constant distribution space.We apply and analyze the method, named the PCD method, in the case of the diffusion equation on a rectangular mesh with refined zones.It has the advantage, compared with other discretization methods, of producing the most compact discrete stencil.Essentially, the PCD method does not make use of the so-called slave nodes that appear in some finite element discretizations with local mesh refinement.Particularly, the graph of the discrete matrix turns out to be the grid itself of the mesh used for the unknown distribution.
Our approach cannot however rely on the results obtained so far because we use piecewise constant approximations not only for the unknown distribution itself but also for its derivatives.This feature makes its mathematical analysis more difficult and requires further developments.To keep them as simple as possible, we restrict the present contribution to the analysis of the two dimensional diffusion equation on rectangular mesh with local mesh refinement.
More specifically, we consider solving the following partial differential equation on a rectangular domain Ω : − div ( p(x) ∇u(x)) + q(x)u(x) = s(x) x ∈ Ω (1.1) where n denotes the unit normal to Γ = ∂Ω and Γ = Γ 0 ∪ Γ 1 .By a rectangular domain we understand any connected subset of the plane with orthogonal sides.It may be L-shaped and need only to be simply connected.We assume that p(x) is bounded and strictly positive on Ω and that q(x) and ω(x) are bounded and nonnegative on Ω and Γ 1 respectively, and we assume that s(x) is in L 2 (Ω).In the case where q(x) and ω(x) would be identically zero, we require that Γ 0 has positive measure.Otherwise, one would extend the results developed below by requiring that s(x) and u(x) be orthogonal to unity in the L 2 (Ω) scalar product, as usually done in such analysis.We use homogeneous boundary conditions for simplicity.The extension of the theory to general boundary conditions does not raise new difficulties.
The PCD Method on Composite Grid

123
The discrete version of this problem will be based on its variational formulation: find u ∈ H such that ∀ v ∈ H a( u, v ) = (s, v) Ω (1.4) where and (s, v) Ω denotes the L 2 (Ω) scalar product.We use here the standard notation for Sobolev spaces as defined in Adams [1]: where A is a linear operator with domain D(A) ⊂ L 2 (Ω) and range R(A) = L 2 (Ω) (distributions belonging to L 2 (Ω) space acting on L 2 (Ω) distributions).So that Equation (1.6) is actually equivalent to stating that: where we have used the bilinear notation < , > rather than the usual scalar product notation in order to stress that A u, s and v are distributions.This holds in particular when v ∈ H and integration by parts is then generally used to get Equation (1.4) which is at the root of the study of the properties of the operator A.
Integration by parts can however be done more generally, for example when v has only piecewise square summable derivatives, with possible jumps at the inner boundaries of the subdomains of Ω on which it has square summable first derivatives (and provided that p(x) is not discontinuous across these boundaries).See for example [7, pp. 94-95].It turns out that the issue may again be written: although the roles of u and v are no more symmetrical in that case, since ∇v may present Dirac behaviors across some inner boundaries on which p ∇u is required to be continuous.
We denote by h the mesh size defined by: h = max( h ℓ ) (ℓ ∈ J) where h ℓ = diam( Ω ℓ ) (ℓ ∈ J) and we denote by h ℓ1 and h ℓ2 , the width and the height of the element Ω ℓ .We define several submeshes on each element Ω ℓ for the representation of v ∈ H 1 (Ω) and its derivatives ∂ i v (i = 1, 2).These representations, denoted v h and ∂ hi v h (i = 1, 2) respectively, are piecewise constant on each of these submeshes (a specific one for each) with the additional requirement for v h that it must be continuous across the element boundaries (i.e.along the normal to the element boundary).In the case of the rectangular meshes considered in the present work, the operators ∂ hi (i = 1, 2) will be finite difference quotients taken along the element edges.They would need to be appropriately adapted for other elements.Figure 1 gives an example of submeshes used to define The PCD Method on Composite Grid on the regions denoted 1, 2 on Fig. 1 (c).In addition, v h must be continuous across element boundaries.Thus for example if the bottom boundary of Ω ℓ is common with the top boundary of where h 1 = h 11 + h 12 , and Remember that v h must be continuous across element boundaries.Thus if the top boundary of Ω ℓ is common with the bottom boundaries of the 2 cells Ω k1 and Ω k2 of widths h 11 and h 12 we have that: For simplicity, we consider a mesh refinement by a factor 2. Higher order ratio can of course be handled in a similar way, see Tahiri [15].

A. Tahiri
Figure 3 provides a particular case of submeshes used to define  The further handling and analysis of our discretization method require appropriate notation for the various spaces involved, which we introduce now.
By X we denote (L 2 (Ω)) 3 with norm ( u, v, w ) 2 X = u 2 + v 2 + w 2 .By Y we denote the subspace of X of the elements of the form By H h0 and H hi (i = 1, 2) we denote the spaces of piecewise constant distributions used to define v h and ∂ hi v h (i = 1, 2), equipped with the L 2 (Ω) scalar product.
We further denote by H h the space H h0 equipped with the inner product: and its associate norm is denoted .h .Clearly H and H h are isomorphic to Y and Y h respectively and we let f and f h denote the bijections of H and H h into X and The PCD Method on Composite Grid 127 By r hi we denote the L 2 -orthogonal projection from L 2 (Ω) onto H hi (i = 0, 1, 2) and we let r h v, for any v ∈ H, denotes the element of H h determined by r h0 v ∈ H h0 .We precise, On Figure 4 we represent the structure of this discretization and the used operators.The motivation for using this space structure is that, while we cannot directly compare the elements of H and H h , we can use the norm of X to measure the distance between elements Figure 5 (top-left) provides an example of a rectangular element mesh with a refined zone in the right upper corner.On the same figure we also represent the meshes H h0 and H hi (i = 1, 2) used to define the piecewise constant distributions v h and ∂ hi v h , i = 1, 2. Each of these meshes defines cells which are useful for distinct purposes.The elements are denoted by Ω ℓ , ℓ ∈ J with boundaries ∂Ω ℓ and closures V ℓ = Ω ℓ .We similarly denote the cells of the other meshes by Ω ℓi , ℓ ∈ J i , i = 0, 1, 2 respectively, with boundaries ∂Ω ℓi and closures V ℓi = Ω ℓi .The measures of these cells will be denoted by It is of interest to note that each node of the mesh may be uniquely associated with a cell of H h0 ; we shall therefore denote them by N ℓ , ℓ ∈ J 0 .Also we note that ∂ hi v h has the following property: for any pair of nodes {P, Q} of the mesh.
At this stage, we do not introduce restrictions on the choice of the discrete mesh.It will be seen however in the convergence analysis that the grid lines of the element mesh should include all material discontinuities, i.e all lines of discontinuity of p(x).
Before closing this section we note that triangular elements may also be introduced.In this way, the PCD method can accommodate any shape of the domain through the combined use of local mesh refinement and triangular elements, see Tahiri [15].We note that the use of rectangular and triangular elements is not a restriction of the PCD discretization.Other elements and other forms of submeshes on such elements can be used, see Beauwens [6].

A. Tahiri
Element mesh H h0 mesh We now define the discrete problem to be solved in H h by: where The discrete matrix is obtained as usual by introducing a basis (φ i ) i∈J0 of the space H h .Expanding the unknown u h in this basis and expressing the variational condition (2.3) by: ) Ω for all i ∈ J 0 whence the linear system A ξ = b with stiffness matrix: A = ( a i j ) = ( a h ( φ j , φ i ) ) , right-hand side b with components b i = ( s , φ i ) Ω , i ∈ J 0 , and unknown vector ξ with components ξ j , j ∈ J 0 .The basis (φ i ) i∈J0 of H h is defined as usual through the conditions: The stiffness matrix is also built as usual by assembling element stiffness matrices.
To this end a few reference elements Ω may be used.

A. Tahiri
On Figure 6 we represent a reference elements Ω for which it is readily seen that the element matrix graph is the element grid itself.On the same figure, we have indicated the offdiagonal element matrix entries a (ℓ) ij along the edges in the case where p(x) = 1.The diagonal entries follow from the formula: The values represented on Figure 6 are obtained by using equation (2.3) and by introducing the local basis of each reference element.We recall that this local basis is reduced to the characteristic function of each volume defined on the Figure 1 (a) , Figure 2 (a) and Figure 3 (a) respectively (in the case where p(x) = 1 and q(x) = 0).

Notation :
We investigate in this section general properties of the discrete space H h and of the possible discrete representations of v ∈ H in H h .We sometimes need to assume that the discretization is regular.We hereby mean that there exist positive constants C 1 and C 2 independent of h such that: By x ℓ = (x ℓ1 , x ℓ2 ), ℓ ∈ J 0 , we denote the grid points of the element mesh, by x ℓE we denote the right neighbor of x ℓ if it exists.By x ℓ + (respectively x ℓ − ) we denote the nearest neighbor of x ℓE and both are located on the same vertical grid line (x ℓ + the top neighbor of x ℓE and x ℓ − the bottom neighbor of x ℓE ).
We denote by J R the subset of J containing the fine element subscripts (elements located in the refined zone).J i Ir , (i = 1, 2) (Ir for irregular) denotes the subset of J containing the subscripts of the irregular element in x i -direction (i = 1, 2).We denote by (r h v) ℓ the value of r h v on Ω ℓ 0 .We split the domain Ω into two subdomains Ω C (the coarse zone) and Ω R (the refined zone) with Ω = Ω C ∪ Ω R .Finally we denote by Ω Ir the union of all irregular elements, ; Ω Ir is a strip in Ω with a width O(h) and has the interface boundary as part of its boundary.The notation C is used throughout the paper to denote a generic positive constant independent of the mesh size.

Discrete Friedrichs inequalities :
The PCD discretization has the following properties which represent a discrete version of the first and the second Friedrichs inequalities and the trace inequality, for the proof we refer to Tahiri [15].
The PCD Method on Composite Grid 131 Lemma 3.1.Let Ω be a bounded polygonal domain.We assume that Γ 0 = ∂Ω .Then, there exists a constant C > 0, independent of h such that: When Γ 0 has a positive measure and Γ 0 = ∂Ω, one may prove the following lemmas.
Lemma 3.2.Let Ω be a bounded polygonal domain.Then, there exists a constant C > 0, independent of h such that: Let Ω be a bounded polygonal domain.Then, there exists a constant C > 0, independent of h such that: The results given in the previous lemmas are independent of the presence or not of the local mesh refinement.We note that, with the PCD discretization, for any pair of nodes of the mesh we can find a path connecting these nodes (succession of mesh grid segments).The proofs of the previous lemmas are based on this property and the property (2.2).

Discrete representation of v ∈ H 1 (Ω) :
Considering now discrete representation of elements of H, we first investigate 2) converges to zero, see for example Brezzi and Fortin [7] and Douglas et al. [11].
We shall prove that Then we can write for all v ∈ C 1 (Ω): for all ε > 0 there exists a h c such that for all h ≤ h c : Then, where n is the unit outward normal vector on ∂Ω ℓ1 and e 1 is the unit vector in x 1 direction.Then, using (3.1), (3.2), (3.7) and (3.8) we get: We next consider the interpolant v I ∈ H h0 of v ∈ H 1 (Ω) ∩ C 0 (Ω), defined by: We try to obtain similar results for v I as those obtained for r h0 v with v ∈ H 1 (Ω).
Proof: By a density argument it is sufficient to prove this in C 1 (Ω).For all ℓ ∈ J 0 and for all x ∈ Ω ℓ0 , we have: where S 1 is an horizontal segment, S 2 is a vertical segment and S 1 ∪ S 2 is a path connecting x and x ℓ .Then, The PCD Method on Composite Grid 133 Then, Integrating this inequality on the cell Then, (3.10) follows immediately.

Case of higher regularity :
We now assume that: , where Ω j is a subdomain of Ω (for example subdomains where p(x) is continuous).For all v ∈ H 2 L (Ω), the interpolant v I ∈ H h0 is defined by (3.9).For all v ∈ H 2 L (Ω), we have (see Brezzi and Fortin [7] and Douglas et al. [11].) r h v and r hi (∂ i v) are defined by: Let Ω be a rectangular bounded open set of R 2 , there exists a constant C > 0 (independent of h) such that: for all v in H 2 L (Ω) (3.16)

134
A. Tahiri Proof: Using (3.10) we have: For the partial derivatives, by a density argument, it is sufficient to prove these in For all ℓ ∈ J 1 and for all v ∈ C 2 (V ℓ1 ) , there exists θ ∈ ] 0, 1 [ such that: Furthermore, we have for all x ∈ V ℓ1 : where S 1 is an horizontal segment, S 2 is a vertical segment and S 1 ∪ S 2 is a path connecting x and (x ℓ + θh ℓ1 ).Then, Integrating this inequality on the cell Ω ℓ1 where Ω j is the subdomain of Ω containing the cell Ω ℓ1 .In the case where V ℓ1 is a subset of Ω 1 ∪ Ω 2 , where Ω 1 and Ω 2 are two subdomains of Ω, we consider the proof in V ℓ1 ∩ Ω 1 and in V ℓ1 ∩ Ω 2 .We get: Then, The same argument can be used for ∂ 2 v(x) − ∂ h2 v I .Then (3.16) follows easily.
For other proofs see Tahiri [15].✷ We note that the results presented in this section are independent of the presence or not of the local mesh refinement.
The PCD Method on Composite Grid In this section we analyze the convergence of the solution u h of the approximate problem (2.3) to the solution u of the continuous problem (1.1).The approximate bilinear form (2.4) is symmetric over H h × H h .By using the Lemma 3.1 or the Lemmas 3.2 and 3.3, we prove the uniform coercivity of the approximate bilinear form.The coercivity of (2.4) implies that the stiffness matrix is positive definite.Moreover, it is also symmetric since (2.4) is a symmetric form.The convergence is analyzed in two steps.First we give the error bounds induced by some specific approximation u d ∈ H h of u ∈ H (u d = r h u or u d = u I ).Such bounds have already been obtained for r h u (Lemma 3.4) and for u I (Lemma 3.6).It remains to give an error bound for u d − u h h between u d and the approximate solution u h .Since H h norm and a h -norm are equivalent, we may equivalently try to bound: Since s(x) is still defined on H h and s(x) is replaced by its value in (1.1).By introducing the following restrictions the expression a ( u , v h ) is well defined.It should be noticed that ∂ 1 v h is reduced to Dirac distributions taken along the edges of the regions where v h is constant (the vertical median line E ℓ of the cell Ω ℓ1 in this case), weighted by the corresponding discontinuity of v h .We required, in the definition of our BVP, that p(x) be a bounded distribution, i.e. p ∈ L ∞ (Ω).The reason was that with ∂ i u ∈ L 2 (Ω), we then also have p(x) ∂ i u ∈ L 2 (Ω) and all operators are clearly defined.However, in the expression of the error, we now see that ∂ 1 v h appear, with Dirac behaviors across specific lines and this is clearly incompatible with coefficients p(x) that would be discontinuous across the same lines.To avoid such situations, we must introduce restrictions on the choice of the mesh, namely that material discontinuities (i.e.discontinuities of p(x)) should never match grid lines of the H h0 mesh.The best practical way to ensure this restriction is to require that material discontinuities be always grid lines of the element mesh.Then we have: where u h is the solution of the problem (2.3) with local mesh refinement and Ω Ir is a strip in Ω.
Proof: We assume that u ∈ H 1 (Ω), by density argument, the proof is only considered for u ∈ C 1 (Ω).We have for all v h ∈ H h and for r h u : Using lemma 3.4, one may write: The other terms A i (i = 1, 2) are most easily analyzed on the cells Ω ℓi (i = 1, 2) of the H hi meshes (i = 1, 2).Being similar in both cases, we consider i = 1.Suppose first that we have a regular case (Fig. 7 (a)).Then, the contribution of an arbitrary cell Ω ℓ1 is : The PCD Method on Composite Grid 137 Then, where Q denotes average on Q, here the H h1 mesh cell Ω ℓ1 or its vertical median line E ℓ , defined by: Taylor expansion and (3.7) give that: It follows that: Now we consider an irregular case (Fig. 7 (b)), the contribution of this irregular cell Ω ℓ1 is: Then, ) we denote the value of v h at the grid point x ℓ + (respectively x ℓ − ) and by E ℓ = E − ∪ E 0 ∪ E + we denote the vertical median of the cell Ω ℓ1 .We add together all terms A ℓ 1 , regular and irregular, and using (3.1) and (3.2) we get: Since we assume that u belongs to H 2 (Ω Ir ), we use the Cauchy − Schwarz inequality and the trace inequality on the strip Ω Ir .Then, Similarly, the same results are derived for A 2 .
Finally we have The PCD Method on Composite Grid

139
Triangle inequality completes the proof.✷ It is of interest to consider here the case where u d = u h , the exact solution of the approximate problem, in the particular case where q(x) = 0.In this case indeed the error must be zero showing that u h is such that: Beyond giving us a relation between the exact and approximate solution, it shows that in this case, the method reduces to the control volume method under its corner mesh version.In more general cases, we may still consider our approach as an interpretation of this box method.The presented method is not a control volume method and cannot in general be correctly framed in that way.However, in the case of a rectangular mesh (whether uniform or not) it can be recovered via the corner mesh version of the control volume method and it does conserve the mass in that case.This does not seem feasible in more general cases as for example along a refined zone.
Remark 4.1.For an element with two refined sides, as on Fig. 3, the presented analysis is still valid.Indeed, the Dirac distribution ∂ 1 v h is also weighted by the inner product of the unit outward normal vector on the boundary cell with the unit vector in x 1 -direction.

Convergence rate :
Under additional regularity for u the exact solution of (1.1), we derive a bound on the convergence rate.Now, we assume that u belongs to H 2 L (Ω).The nodal values of u are well defined for each point of Ω.We assign the unknowns of the approximate solution u h to the nodes x ℓ .We denote by u ℓ the nodal value of u at the grid point x ℓ .Also here, the error analysis is done in two steps.The first one is f (u) − f h (u I ) X that is given in Lemma 3.6 and the second one is the error u I − u h h .As previously we try to bound u I − u h a h defined in (4.1) with u d = u I .Theorem 4.2.Let Ω be a rectangular bounded open set of R 2 .Assume that the unique variational solution u of (1.1) belongs to H 2 L (Ω), there exists a constant C > 0 ( independent of h), such that: where u h is the solution of the problem (2.3) without local refinement.C > 0 ( independent of h), such that: where u h is the solution of the problem (2.3) with local mesh refinement and Ω Ir is a strip in Ω.
Proof: The proofs of Theorems 4.2 and 4.3 are similar, we give only the proof of Theorem 4.3.We have for all v h ∈ H h and for u d = u I : where A i (i = 0, 1, 2) are defined in (4.2), (4.3) and (4.4) with u I instead of r h u.By Lemma 3.6 one may write: As previously, we consider the contribution of an arbitrary regular cell Ω ℓ1 : Note that x ℓ and x ℓE are in V ℓ1 , and E ℓ is a subset of V ℓ1 .
Taylor expansion gives on a closed domain K , for all x , x 0 ∈ K: where H(v)(x) denotes the Hessian matrix of v at point x.Using (4.7) for (x = x ℓ , x 0 = x) and for (x = x ℓE , x 0 = x), ∀ x ∈ E ℓ , subtracting one from the other: where: The PCD Method on Composite Grid 141 Integrating over E ℓ and since ∇u(x) • (x ℓE − x ℓ ) = (h ℓ1 ∂ 1 u(x) , 0), we get: Then, where E t ℓ is the transformation of E ℓ by the first change of variable.Since E t ℓ ⊂ E ℓ and using the second change of variable, we obtain: We denote by Ω 1 ℓ1 the left half part of Ω ℓ1 , and using Cauchy-Schwarz inequality: In similar way, we can get the same error bound for E ℓ |φ ℓE | ds.Therefore, Now we consider an irregular case (Fig. 7 (b)), the contribution of this irregular cell Ω ℓ1 is: Then, where Ω 1 and Ω 2 are two subdomains of Ω, we consider the proof in We add together all terms A ℓ 1 , regular and irregular, and using (3.1) and (3.2) we get: The PCD Method on Composite Grid 143 Using the Cauchy − Schwarz inequality and the trace inequality, we get: Similarly, one may write for A 2 : Therefore, Triangle inequality and Lemma 3.6 complete the proof.✷ Remark 4.2.The presented method is well adapted for multilevel local refinement.Successive local mesh refinement can be handled in similar way.The results established here are still valid.

Concluding remarks
The first and main issue of the present work is the introduction of the PCD discretization method that relies on the systematic use of piecewise constant distributions to represent the unknown distribution as well as its derivatives, each one on a specific mesh.This method was only formulated for rectangular grids because the choice of the appropriate meshes is rather simple in that case and this choice led us to a new Ritz-Galerkin formulation of the corner mesh box scheme.Its extension to more general (irregular) meshes clearly raises geometrical difficulties but a priori no basic principle objection.
Staying still with rectangular meshes we have investigated the question of feasibility of introducing local mesh refinement without slave nodes, with the issue of getting the most compact possible discrete stencil.The main question to solve here was to prove the convergence of the resulting scheme.
Our conclusion is that, provided that the solution belongs to H 1 (Ω) ∩ H 2 (Ω Ir ), there is always convergence in energy norm.Further, if the solution is locally H 2regular, we then showed an O(h 1/2 ) convergence rate.
One can also use slave node techniques to improve the convergence at the expense of increased complexity, see Tahiri [15].In our study of this method, we have used both theoretical and numerical investigations.The present work summarizes our theoretical findings.For the numerical results, we refer to Tahiri [15,17].Our numerical results are in agreement with the theoretical results presented here.Furthermore, the numerical results show that the L 2 -error u − u h has an O(h)-convergence rate independently of the presence or not of the local mesh refinement.

125 on the regions denoted 1 ,
2 on Fig. 1 (b) and ∂ h2 v h | Ω ℓ is similarly the piecewise constant distribution with constant values:

Figure 4 :
Figure 4: Structure of the discretization analysis

Theorem 4 . 3 .
Let Ω be a rectangular bounded open set of R 2 .Assume that the unique variational solution u of (1.1) belongs to H 2 L (Ω), there exists a constant A. Tahiri regular rectangular element Ω ℓ , i.e. an element from a regular rectangular mesh.vh| Ω ℓ is the piecewise constant distribution with 4 values v hi on the regions denoted i with i = 1, ..., 4 on Fig.1 (a).∂ h1 v h | Ω ℓ is the piecewise constant distribution with constant values: when the boundary of a refined zone covers two different sides ofΩ ℓ .In this case v h | Ω ℓ assumes 6 values v hi , i = 1, ..., 6 , ∂ h1 v h | Ω ℓ 3 values, where h 1 = h 11 + h 12