Numerical Computations of the PCD Method

abstract: The PCD (piecewise constant distributions) method is a discretization technique of the boundary value problems in which the unknown distribution and its derivatives are represented by piecewise constant distributions but on distinct meshes. It has the advantage of producing the most sparse stiffness matrix resulting from the approximate problem. In this contribution, we propose a general triangulation with the PCD method by combining rectangular elements and triangular elements. We also apply this discretization technique for the 2D elasticity problem. We conclude by presenting the numerical results of the proposed method for the 2D diffusion equation.


Introduction
We present a discretization technique of boundary value problems (BVP) in which the unknown distribution as well as its derivatives are all represented by piecewise constant distributions, each one on a specific mesh.The interest for piecewise constant approximations has been highlighted since the early days of partial differential equations as a discretization techniques, being more or less explicitly 40 A. Tahiri at the root of the control volume method widely used in engineering applications, particularly in computational fluid dynamics.This motivated the early analysis by [1], [2], [5] and [14] which raises its persistent interest as it is illustrated in many more recent works [4], [6] , [7] and [8].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 requires the introduction of distinct meshes to represent each distribution and new mathematical tools for the convergence analysis of the resulting discrete schemes.

Formulation of the problem
To keep the presentation of this discretization as simple as possible, we restrict the present contribution to the analysis of the 2D diffusion equation on a nonuniform rectangular mesh.The convergence analysis and the technical results of the PCD method can be found in [9] and [12].The extension of the presented method on polygonal or L-shaped domain does not raise any difficulties.
We consider solving the following BVP on a rectangular domain Ω : − div( p(x) ∇u(x)) + q(x)u(x) = s(x) x ∈ Ω (2.1) where n denotes the unit normal to Γ = ∂Ω and Γ = Γ 0 ∪ Γ 1 .We assume that p and q are in L ∞ (Ω), that p(x) is strictly positive on Ω, that q(x) is nonnegative on Ω and that we have a well posed problem.
The discrete version of this problem will be based on its variational formulation: where and (s, v) Ω denotes the L 2 (Ω) scalar product.

PCD spaces
The PCD discretization splits the open domain Ω under investigation into elements Ω ℓ , ℓ ∈ J , open subsets of Ω, such that

41
and defines several sub-meshes on each element Ω ℓ for the representation of v ∈ H 1 (Ω) and of 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 sub-meshes (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).Here, we consider rectangular meshes and the operators ∂ hi (i = 1, 2) are finite difference quotients taken along the element edges.
Similarly, ∂ h2 v h | Ω ℓ is the piecewise constant distribution with constant values: on the regions denoted 1, 2 on Fig. 1 (c).
In addition, v h must be continuous across element boundaries.
It is of interest to note that ∂ hi v h has the following property: for any pair of nodes {P, Q} of the mesh.
Element mesh H h0 mesh The spaces associated with this discretization are defined as follows.E denotes (L 2 (Ω)) 3 with norm ( u, v, w H h0 and H hi (i = 1, 2) 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 F and F h respectively and we let f and f h denote the bijections of H and H h into E and E h ( F = f (H) and We further denote by p h the canonical injection of E h into E.These operators are represented on Fig. 3.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 E to measure the distance between elements Figure 2 (top-left) provides an example of a rectangular element mesh.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 ; we similarly denote the cells of the other meshes by Ω ℓi , ℓ ∈ J i , i = 0, 1, 2 respectively.It is of interest to note that each node of the element mesh may be uniquely associated with a cell of H h0 ; we therefore denote them by N ℓ , ℓ ∈ J 0 .The mesh size, denoted h, is defined as:

General triangulation
We note that triangular elements may also be used in the PCD discretization.For this purpose we consider a rectangular triangular element Ω ℓ and we define the sub-meshes, on such element, used to define the representations v h , ∂ h1 v h , and and the representation ∂ h2 v h assumes 1 value (Fig. 4 (c)) where h ij is the distance between the nodes N i and N j ( 1 ≤ i , j ≤ 3 ).Also here we require the continuity of v h across the element boundaries.
In this way, the method can accommodate any shape of polygonal domains through the combined use of rectangular elements and triangular elements.
We note that the use of rectangular and triangular elements is not a restriction of the PCD discretization.Other elements and other forms of sub-meshes on such elements can be used, see [3].We note that the property (2.6) is still valid even if we have combine the rectangular elements and the rectangular triangular elements.Also this property is still valid even if we introduce a local mesh refinement, see [9], [10] and [12].

PCD equations
We now define the discrete problem to be solved in 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 u h = i∈J0 ξ i φ i and expressing the variational condition (2.8) by Whence, the linear system A ξ = b with stiffness matrix A = ( a i j ) = ( a h ( φ i , φ j ) ) , right hand side b with components b j = ( s , φ j ) Ω , j ∈ J 0 and unknown vector ξ with components ξ i , i ∈ J 0 .
It should be stated that the presented method has the advantage of producing the most compact schemes.Also, it leads to the most sparse stiffness matrix resulting from the approximate problem.

Discrete Friedrichs inequalities
Here we present the following lemmas that represent a discrete version of first and second Friedrichs inequalities and trace inequality for the proposed discretization.
Lemma 2.1.Let Ω be a connected bounded polygonal domain.We assume that Γ 0 = ∂Ω .Then, there exists a constant C > 0, independent of h such that:

10)
Proof: This lemma is a discrete version of the first Friedrichs inequality.We give its proof on a square (0, a) × (0, a) since we can include all connected bounded polygonal domain in a square.A. Tahiri Letting P and Q be two nodes such that P ∈ ∂Ω and Q ∈ Ω ⊂ (0, a) × (0, a) , there is a path of the form S = { P 0 , P 1 , P 2 , ..., P l } connecting P 0 = P and P l = Q (made up of a succession of mesh grid segments, see Fig. 5) such that where i = 1 on horizontal segments and i = 2 on vertical segments.Subdividing S into S = S 1 ∪ S 2 where S 1 is the union of horizontal segments and S 2 the union of vertical segments, one may write Then, Therefore ( by Integrating this inequality on the domain Ω This Lemma and this proof are still valid if Γ 0 has a positive measure and Γ 0 = ∂Ω. In such case, we can choose the node P such that P ∈ Γ 0 and we prove (2.10) by the same argument.
In similar way, we can prove the following lemmas which give a discrete Friedrichs second inequality and a discrete trace inequality for the proposed discretization in the case where Γ 0 has a positive measure, see [9].
Lemma 2.2.Let Ω be a connected bounded polygonal domain.Then, there exists a constant C > 0, independent of h such that:

11)
Numerical Computations of the PCD Method 47 Lemma 2.3.Let Ω be a connected bounded polygonal domain.Then, there exists a constant C > 0, independent of h such that: Using these lemmas we can prove that a h ( u h , v h ) is uniformly coercive on H h .Therefore, the associated norm is uniformly equivalent to the norm of H h , i.e. there exists positive constants α , β independent of the mesh size h such that: On the other hand, since p(x) and q(x) are in L ∞ (Ω), it is clear that a h ( u h , v h ) is uniformly continuous on H h .Therefore, the problem (2.8) has a unique solution.

Convergence analysis
As mentioned in Sect.2, the error between the solution u ∈ H of (2.1) and the discrete solution u h ∈ H h of (2.8) will be measured by the distance between their representations in E, namely As usual, the bounds that can be obtained depend on the regularity of u.Here, we assume that u ∈ H 2 (Ω).In that case, u is continuous on Ω and we can then define its interpolant u I in H h through and rely on interpolation theory in H h .More specifically, we can then use the following result proved in [9], [10] and [12].
Lemma 3.1.Under the general assumptions and notation defined above, there exists a positive constant C independent of the mesh size h such that ) By this result, it remains to bound the error between u I and the approximate solution Since the a h -norm is uniformly equivalent to the H h -norm, we may equivalently try to bound u I − u h a h defined by Since s(x) is still defined on H h and s(x) is replaced by its value in (2.1).In the term a ( u , v h ), the derivatives of v h are understood in distributions sense.By introducing some restrictions the expression a ( u , v h ) can be defined.More explicitly: To bound (3.5), we may try to bound separately the terms: The first term is bounded by using Lemma (3.1).The last terms are most easily analyzed and bounded on the cells Ω ℓi of the H hi meshes (i = 1, 2).Being similar in both cases, we consider i = 1, the contribution of an arbitrary cell Ω ℓ1 is: where (v 2 − v 1 ) is the jump of v h through the vertical median line E ℓ of the cell Ω ℓ1 .
It should be noticed that ∂ 1 v h (respectively ∂ 2 v h ) is reduced to Dirac distributions taken along the edges of the regions where v h is constant weighted by the corresponding discontinuity of v h , the vertical median line of the cell Ω ℓ1 of the H h1 meshes in this case (respectively the horizontal median line of the cell Ω ℓ2 of the H h2 meshes).However, in the expression of the error, we now see that ∂ 1 v h and ∂ 2 v h appear with Dirac behaviors across these lines.This is clearly incompatible with the coefficient 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.By introducing this restriction the expression a ( u , v h ) is well defined.Under the mentioned restrictions and the regularity of u the exact solution of (2.1), we can give the following Theorem, for the proof see [9] and [12].
Theorem 3.2.Assume that the solution u of (2.1) belongs to H 2 (Ω).Then, there exists a constant C > 0 (independent of h), such that: where u h is the solution of the problem (2.8).
The result of the previous theorem is still valid if the solution of the problem (2.8) is only locally H 2 -regular.Under higher regularity assumptions an error bound of order O(h 2 ) can be obtained.If the solution u is only in H 1 (Ω), we can still prove the convergence of u h to u, see [9] and [12].

Application: Elasticity problem
We consider the boundary value problem: find u = ( u 1 , u 2 ) in (H 1 (Ω)) 2 which satisfies the following equation in the domain Ω.
where n denotes the unit normal to Γ = ∂Ω and Γ = Γ 0 ∪ Γ 1 .We assume that µ and λ are two constants such that µ ≥ µ 0 > 0 , λ > 0 , s = ( s 1 , s 2 ) in (L 2 (Ω)) 2 .We assume that this problem is well posed and has a unique solution in H = (H 1 Γ0 (Ω)) 2 equipped with the product norm: The variational formulation of this problem can be written: where Also here we use the same discretization technique to define an approximate problem of (4.4).The space H h0 is used to define the representation for each component v i (i = 1, 2) of v.The space H h1 (respectively H h2 ) is used to define the representation for the partial derivatives components For this application (elasticity problem), the discrete space H h is the space H h0 × H h0 equipped with the norm: We now define the approximate problem to be solved in H h by: where and s(v h ) is still defined by (4.6).
Also here this discretization technique leads to the most sparse stiffness matrix from the approximate problem (4.7).
By the same analysis as in the previous section, we can bound the error between u = ( u 1 , u 2 ) the exact solution of (4.4) and u h = ( u 1 h , u 2 h ) the exact solution of (4.7).Then, we can write the following error estimate:

Numerical experiments
In this section we present experimental results concerning the PCD discretization.We consider the domain Ω as the unit square ]0, 1[×]0, 1[.We denote by A the stiffness matrix arising from the approximate problem (2.8), u h denotes its solution and N denotes the number of unknowns.u is the exact solution of (2.1) and u I its interpolant in H h .We consider the following error estimators: the relative L 2 -error ε r 0 and the relative H 1 -error ε r 1 , defined by: We consider the problem: −div(∇u(x, y)) = s(x, y) in Ω and u(x, y) = 0 on ∂ Ω .
We consider 3 examples, in the first one (respectively the second) we choose the source term s(x, y) such that the exact solution is u 1 (x, y) = x(x − 1)y(y − 1) (respectively u 2 (x, y)), the results of this example are presented in Table 1 (respectively Table 2), where In Table 3 we present the results of the third example where the exact solution is:    From Tables 1-3 we observe a monotonic improvement of the accuracy in both error estimators ε r 0 and ε r 1 , they decrease when N increases.
From Table 1-2, we observe that ε r 0 and ε r 1 are reduced by a factor of 4 when the mesh size is reduced by a factor of 2.Then, we have an O(h 2 ) convergence rate, for the both norms, since we have a very smooth solutions.That proves, under higher regularity assumption, the presented method has the standard O(h 2 ) convergence rate.In the third example (Table 3), We observe that ε r 0 and ε r 1 are reduced by a factor of 2 when the mesh size is reduced by a factor of 2. Since the solution u 3 (x, y) has a sharp peak and an important variation in Ω 1 than the remaining part of Ω, the presented method has only an O(h) convergence rate.For the three examples, our numerical results are in agreement with the theoretical results given in Sect.3. We can improve the numerical results for the example 3 by introducing a local mesh refinement without slave nodes which does not raise new complexity, see [9], [11] and [12].Furthermore, the PCD method can accommodate any shape of the domain by using triangular elements and still produces a reduced size of the approximate problem.More explicitly, in Table 3 we must have an uniform mesh with 65025 nodes to get ε r 0 = 2.817 × 10 −2 and the local L 2 -error around the point (0.5, 0.5) is 3.092 × 10 −2 .With a local mesh refinement we get for ε r 0 = 2.434 × 10 −2 only with the use of 28545 nodes and the local L 2 -error around the point (0.5, 0.5) is 7.547 × 10 −4 .For more numerical tests and for the best local refinement strategy, we refer to [9], [11] and [13].
We illustrate on Fig. 6 an example of general triangulation of a polygonal domain Ω by combining rectangular elements, triangular elements and local mesh refinement.

Concluding remarks
The main issue of the present work is the presentation of a BVP discretization method on polygonal domain.It is based on the use of piecewise constant distributions to represent the unknown distribution as well as its derivatives on distinct meshes.The PCD method has the advantage of producing the most sparse stiffness matrix resulting from the approximate problem, and has the standard first order estimate under the H 2 -regularity of the exact solution.The method is presented in the case of the diffusion equation but it can be applied in the other cases (elasticity problem for example), see [9] and [10].We note that, with the PCD method we can introduce a local mesh refinement without slave nodes and it is still producing the most sparse stiffness matrix, see [9], [11], [12] and [13].This local refinement improves the accuracy of the approximate solution without additional difficulties.Also, we note that we can combine the use of local mesh refinement and triangular elements in order to be able to follow any shape on the domain Ω and to reduce the size of the approximate problem.
a rectangular element Ω ℓ are represented on Fig. 1. v h | Ω ℓ 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:

Figure 3 :
Figure 3: Structure of the discretization analysis

Figure 5 :
Figure 5: Example of path Let the bounded open square (0, a) × (0, a) with sides parallel to the axes x 1 and x 2 .

Figure 6 :
Figure 6: General triangulation of the domain Ω Thus for example if the bottom boundary of Ω ℓ is common with the top boundary of Ω k , one must have that v h1

Table 1 :
Results of the first example

Table 2 :
Results of the second example

Table 3 :
Results of the third example