VAGO method for the solution of elliptic second-order boundary value problems

Mathematical physics problems are often formulated using differential oprators of vector analysis - invariant operators of first order, namely, divergence, gradient and rotor operators. In approximate solution of such problems it is natural to employ similar operator formulations for grid problems, too. The VAGO (Vector Analysis Grid Operators) method is based on such a methodology. In this paper the vector analysis difference operators are constructed using the Delaunay triangulation and the Voronoi diagrams. Further the VAGO method is used to solve approximately boundary value problems for the general elliptic equation of second order. In the convection-diffusion-reaction equation the diffusion coefficient is a symmetric tensor of second order.


Introduction
In the theory of difference schemes [1] the problem of approximation of a differential problem via grid one received primary emphasis. A universal method of balance (integro-interpolation method) is used as the basic one to construct discrete analogs. It was proposed by Samarskii [2] and currently it is known as the control volume method [3,4,5]. A difference scheme is constructed using the balance method via integration of the initial equation over a control volume -a part of the computational domain adjacent to the specified grid point.
In solving mathematical physics problems on general unstructured grids there is often used the Delaunay triangulation which displays some optimal features [6]. For the Delaunay triangulation it is natural to consider as control volumes a Voronoi polygon -a set of points lying closer to the considered grid point in compare with all others. Consistent approximations of the vector analysis operators on the basis of the Delaunay triangulation were constructed in [7] for some problems of mathematical physics. A systematic investigation of approxinations for gradient and divergence operators on dual Voronoi-Delaunay partitionings has been conducted in [8,9] for model 2D and 3D problems.
The present work continues our studies on Vector Analysis Grid Operators method. In paper [10] there are discussed problems of constructing consistent approximations for basic operators of vector analysis -gradient, divergence amd rotor operators. Special representations were derived for truncation errors that allowed to obtain the corresponding estimates for the convergence rate of the approximate solution to exact one. Possibilities of the VAGO method have been demonstrated on some scalar and vector problems. Difference schemes for steady-state convection-diffusion problems were constructed on the basis of the dual Voronoi-Delaunay partitioning in [11,12]. In these problems (see, e.g., [13]) the emphasis is on properties of discrete operators for the convection and diffusion transport. The convective terms are written [14,12] in the divergent (conservative), non-divergent (characterictic) and symmetric forms. In [10] there were designed VAGO approximations as well as investigated truncation errors for the Dirichlet problem for stationary and usteady convection-diffusion problems for isotropic media.
Paricular emphasis should be placed on approximation of diffusion transport operator for an anisotropic medium where the diffusion coefficient is a symmetric tensor of second order. In this case we should control at the discrete level not only the seld-adjoint property of the grid diffusion operator but the monotonicity and conservation properties for the corresponding grid approximatioms, too [1]. The present-day view on approximation of problems with mixed derivatives on rectangular grids is given in [15]. Prop-erties of standard finite-element approximations on triangular or tetrahedral meshes are discussed in [16] in detail. To obtain the monotonicity property, a modification of standard linear Galerkin finite element discretizations of the Laplace operator was performed. Much attention is given to constructing control volume schemes for the approximate solution of convection-diffusion boundary value problems with an anisotropic diffusion coefficient (see [17] and references therein). In particular (see also [18]), the control volume schemes are developed with an anisotropic diffusion for approximating a local discrete gradient. The results of work [8] are generalized in [19] to the 2D div-curl system in anisotropic media.
The outline of this paper is the following. In Section 2, the boundaryvalue problem for the general elliptic equation of second order is formulated in an invariant form as the problem for a convection-diffusion equation based on using vector analysis operators -divergence and gradient. The diffusion coefficient is a symmetric positive definite tensor of second order. On dual Voronoi-Delaunay partitionings (Section 3) there are introduced the corresponding spaces of grid functions. Using the VAGO method the difference operators of divergence and gradient are constracted and investigated in Section 4. Difference schemes for the considered convection-diffusion-reaction problem are studied in Section 5. Approximation of the diffusion transport operator with an anisotropic diffusion coefficient has reseived much consideration.

The boundary value problem
The boundary value problem is considered for an elliptic equation of second order. In n-dimensional bounded domain Ω with smooth enough boundary ∂Ω function u(x), x = (x 1 , ..., x n ) satisfies the following equation The coefficients at the higher order derivatives meet the ellipticity condition: (2) For other coefficients of equation (1) we have Assume that the right hand side of the equation can be represented as follows The simplest boundary value problem is considered for (1) with the following boundary conditions In Hilbert space L 2 (Ω) the scalar product and norm are introduced like this In H(grad) we introduce For the solution of problem (1) -(5) the following apriori estimate takes place (see, e.g., [20]) Unstructured computational grids are often used to solve approximately boundary problems of mathematical physics. In some cases it is not suitable to use the coordinate-wise form of the considering equations and boundary conditions in one or another coordinate system. It seems more reasonable to employ the invarinat formilation of the problem. In this case the considered problem is represented via operators of vector and tensor analysis.
Each node x D i , i = 1, 2, . . . , M D , connect a certain part of the computational domain, namely, the Voronoi polyhedron or its part belonging to Ω. The Voronoi polyhedron (polygon) for a separate node is a set of points lying closer to this node than to all the other ones: where | · | is the Euclidean distance. Each vertex x V k , k = 1, 2, . . . , M V of the Voronoi polyhedron is associated with the tetrahedron constructed by the appropriate nodes contacting the Voronoi polyhedrons. We will assume that all vertices of the Voronoi polyhedrons lie either inside the computational domain Ω or on its boundary ∂Ω. These tetrahedrons determine the Delaunay triangulation -a dual triangulation to the Voronoi diagram. The D-grid in the domain Ω is determined by the set of nodes (vertices of tetrahedrons of the Delaunay triangulation) x D i , i = 1, 2, . . . , M D , the V -grid is defined by the set of nodes (vertices of polyhedron of the Voronoi diagram) We mark a separate tetrahedron D k of the Delaunay triangulation. This tetrahedron is identified by the number k of the Voronoi polyhedronm vertex, For common planes of the tetrahedron we use the notations The boundary of the computational domain ∂Ω consists of the planes of Delaunay tetrahedrons. Let We associate the Voronoi polyhedron V i , i = 1, 2, . . . , M D with the node of the main grid i. Thus, we have In this case, m = 0 means that the tetrahedron D k contacts the boundary. We define also the neighbors for each Voronoi polyhedron V i , i = 1, 2, . . . , M D : We assume that the introduced Delaunay triangulation and the Voronoi diagram are regular [21]. For the notations the regularity condition of the Delaunay triangulation is Likewise, for the Voronoi diagram we have We will approximate the scalar functions of the continuous argument by the scalar grid functions that are defined in the nodes of the D-grid or in the nodes of the V -grid. We denote by H D the set of grid functions defined on the D-grid For the functions y(x) ∈ H D , vanishing on the boundary ∂ω, we define We consider the scalar product and the norm for the scalar grid functions from H D by This scalar product and the norm are grid analogs of the scalar product and the L 2 (Ω)-norm for the scalar functions of the continuous argument.
To determine the vector field in the control volume, it is natural to use the components of the sought function normal to the corresponding planes of the control volume. Choosing the initial and the final node, we connect with each tetrahedron edge D k , k = 1, 2, . . . , M V or polyhedron edge V i , i = 1, 2, . . . , M D , a vector -the directed edge. For Delaunay triangulation the normals to the planes are the directed edges of Voronoi diagram and vice versa. For the approximation of the vector functions thereby we can use projections of the vectors on the directed edges. We will further use exactly this variant with the description of the vector field in the control volume by means of vector projections on the edges of the control volume.
We will orient the Delaunay triangulation edges by the unit vector directed from the node with a smaller number to the node of a larger number. The vector function y(x) on the Delaunay triangulation is defined by the components that are given in the middle of the edges For the Delaunay triangulation used and the Voronoi diagram we define the length of the edges in the following way: We denote by H D the set of grid vector functions determined by the components y D ij , i = 1, 2, . . . , M D , j ∈ W V (i) that are given in the middle of the edges. In a similar way we denote by H V the set of grid vector functions defined by the components y V km , k = 1, 2, . . . , M V , m ∈ W D (k). If the tangential components of the grid vector functions y ∈ H D vanish on the boundary, we define Consider the scalar product and norm in H D .

Gradient and divergence approximations
Let us present the main results (see [10] for deatails) connected with approximation of gradient and divergence operators on the basis of the Delaunay triangulation and Voronoi diagram.
The set of grid functions H D can be the domain of definition of the grid gradient operator. We denote the grid gradient operator by grad D . Taking into account the chosen edge orientation, at the points x D ij we set i.e., the range of values of the operator grad D : H D → H D is the set of vector grid functions H D . In (8), we use the following notation: For the truncation error of the grid operator grad D we have provided that u(x) is a sufficiently smooth functions. For the Voronoi polyhedron this equality is written in the following form: where n V ij is the normal to the edge ∂V ij outside with respect to V i . To construct the grid operator div D : H D → H D we use the elementary formulas of integration for the left-and right-hand sides (10). This leads to the grid analogs of the operator div in the form of for the truncation error. The truncation error for the grid divergence operators div D equals to O(1). However, there does exist a special divergence expression of the truncation error saving the situation in the case of approximation of problems of mathematical physics. Let us highlight the adjoint property of the considered grid operators of gradient and divergence. Taking into account the above notations we have From (13) it follows that div * D = − grad D on the set of grid functions v ∈ H 0 D and y ∈ H D .

Approximation of the convection-diffusion-reaction problem
In domain Ω we consider the regular D-grid and V -grid. Taking into account boundary condition (5), we will find the approximate solution of problem (5), (7) as grid function y(x) ∈ H 0 D . A particular discussion should be given to approximation of convective, and moreover, diffisive terms in equation (7).
We consider approximations of the convective transport operator both in the divergent and in non-divergent (characteristic) form: at specified vector field v(x), x ∈ Ω.
The tangential components are approximated at the edge midpoints of the Delaunay triangulation as follows For scalar functions it is natural to introduce The direct approximation of the convective transport operator in the divergent form results in (15) As for the convective transport operator in the non-divergent form, its approximation is performed [12] using the grid analog of the following equality For the trancation error similarly to (12) we have Let us construct a grid analog of the diffusion transport operator Du = div(D grad u), x ∈ Ω.
We put formally D D y = div D (D D grad D y). Above (see (8)) we presented the grid approximation for the projection of the gradient operator onto edges of the Delaunay triangulation. It is enough for approximation of the diffusion operator in isotropic media (with a scalar diffusion coefficient). For general elliptic equations with mixed derivatives (1) (with a tensor diffusion coefficient) it is necessary to approximate the gradient not only along the edges of triangulation. To investigate this problem, let us start from the 2D case (n = 2, Fig.1).
Here i, j, k, m -nodes of the Dalaunay triangulation, i, j, k i, j, mtwo adjacent triangles of this triangulation with common edge (i, j), p qnodes of the Voronoi diagram. For nodes of the Voronoi diagram there are used notations x V p = x D ijk and x V q = x D ijm . Highlighted in Fig.1 edges of the Delaunay triangulation and Voronoi diagram are orthogonal to one another. Introduce at the point of their intersection (at point x = x D ij ) local coordinate system (s 1 , s 2 ) with direction unit vectors e 1 , e 2 ).
Divergence of a vector field at the discrete level (see (11)) is defined using the tangential component. To approximate the diffusion transport, it is necessary to evaluate the diffusion tensor (18). Assume for the introduced local coordinate system Coefficients (d D αβ ) ij , α, β = 1, 2 are calculated via specified smooth enough coefficients d αβ (x), α, β = 1, 2 of equation (1). If local coordinate system (s 1 , s 2 ) coincides with (x 1 , x 2 ) we have, for instance, In the general case the recalculation rule for a tensor of second order should be taken into account when going to a new orthogonal coordinate system. Using the above notations in local coordinate system (s 1 , s 2 ) the tangential component of the flow can be written as corresponding nodes of the Voronoi diagram (vortices of the Voronoi polygon) we again use notations x V p and x V q . The line connecting these vortices is perpendicular to triangle i, j, k and intersects it at point r (Fig.2). Point s is the midpoint of edge (i, j). Stright line (r, s) is perpendicular to both edge (i, j) of the Delaunay triangulation and edge (p, q) of the Voronoi diagram. On the basis of these three lines (i, j), (p, q) and (r, s) highlighted in Fig.2 we construct local coordinate system (s 1 , s 2 , s 3 ) with direction unit vectors e 1 , e 2 , e 3 ).
Similarly to the 2D case we can evaluate now the diffusion tensor in the local coordinate system For the smooth enough coefficients d αβ (x), α, β = 1, 2, 3 of equation (1) and in case of coincidence of local coordinate system (s 1 , s 2 , s 3 ) with (x 1 , x 2 , x 3 ) it is possible to take (d D αβ ) ij = d αβ (x D ij ), α, β = 1, 2, 3. In local coordinate system (s 1 , s 2 , s 3 ) the flow component along edge (i, j) of the Delaunay triangulation is equal to (D grad D y) D ij = (d D 11 ) ij (grad D y) 1 + (d D 12 ) ij (grad D y) 2 + (d D 13 ) ij (grad D y) 3 .
The tangential component of gradient (grad D y) 1 is calculated in accordance with (20) whereas for component (grad D y) 2 (along the edge of the Voronoi diagram) approximation (21) is used. For the third term in (26) we employ (grad D y) 3 =ȳ r −ȳ s l rs , where l rs -distance between points r and s. Values of scalar functionȳ r ,ȳ s in (27) are approximated with the second order via (22)-type formulas using values defined at nodes of the Delaunay triangulation. So, representation (23) for the truncation error of the flow is valid in the 3D case, too. Hence we aso obtain grid equation (24) with representation (25) for the truncation error. On the basis of this divergent representation for the error we can prove convergence of the discrete solution to the exact one [10] in the grid analog of space H(grad) with the first order with respect to h.