Numerical solution of 2D boundary value problems on merged Voronoi-Delaunay meshes

TL;DR

This paper presents a numerical method for solving 2D boundary value problems using merged Voronoi-Delaunay meshes, enabling efficient approximation of differential operators on irregular grids.

Contribution

It introduces a novel mesh structure combining Voronoi and Delaunay elements with orthodiagonal quadrilaterals for improved numerical solutions.

Findings

Effective approximation of gradient and divergence operators.

Successful application to steady-state diffusion-reaction problems.

Mesh structure facilitates handling anisotropic media.

Abstract

Computational technologies for the approximate solution of multidimensional boundary value problems often rely on irregular computational meshes and finite-volume approximations. In this framework, the discrete problem represents the corresponding conservation law for control volumes associated with the nodes of the mesh. This approach is most naturally and consistently implemented using Delaunay triangulations together with Voronoi diagrams as control volumes. In this paper, we employ meshes with nodes located both at the vertices of Delaunay triangulations and at the generators of Voronoi partitions. The cells of the merged Voronoi-Delaunay mesh are orthodiagonal quadrilaterals. On such meshes, scalar and vector functions, as well as invariant gradient and divergence operators of vector calculus, can be conveniently approximated. We illustrate the capabilities of this approach by solving a steady-state diffusion-reaction problem in an anisotropic medium.

Tables

Table 1. Table 1: Mesh characteristics.

Mesh	D-mesh ($M_D$)	min/max angle	V-mesh ($M_V$)	MVD-mesh ($M$)
1	16	42.7 / 81.2	30	35
4	36	42.9 / 83.3	70	87
7	129	42.2 / 85.3	256	346
10	427	43.9 / 87.2	852	1204
13	1266	42.9 / 85.5	2530	3665
16	4432	40.3 / 85.4	8862	13047
19	16879	40.8 / 85.4	33756	50150
22	64601	41.3 / 85.3	129200	192848

Equations

$$-\mathrm{div}\left(K(\mathbf{x})\mathrm{grad}u(\mathbf{x})\right)+r(\mathbf{x})u(\mathbf{x})=f(\mathbf{x}),\mathbf{x}\in \Omega .$$

$$u(\mathbf{x})=0,\mathbf{x}\in \partial \Omega .$$

$$k_{\alpha \beta }(\mathbf{x})=k_{\beta \alpha }(\mathbf{x}),\alpha ,\beta =1,2,\sum _{\alpha ,\beta =1}^2{k}_{\alpha \beta }(\mathbf{x}){\zeta }_{\alpha }{\zeta }_{\beta }\ge \kappa \sum _{\alpha =1}^2{\zeta }_{\alpha }^2,\kappa >0,\mathbf{x}\in \Omega ,$$

$$(v,w)={\int }_{\Omega }v(\mathbf{x})w(\mathbf{x})d\mathbf{x},\Vert v\Vert =(v,v{)}^{1/2}.$$

$$\Vert \mathrm{grad}u\Vert \le \frac{c}{\kappa }\Vert f\Vert$$

$$\mathbf{q}(\mathbf{x})=-K(\mathbf{x})\mathrm{grad}u(\mathbf{x}).$$

$$\mathrm{div}\mathbf{q}+r(\mathbf{x})u(\mathbf{x})=f(\mathbf{x}),\mathbf{x}\in \Omega .$$

$${\overline{\omega }}^D=\{\mathbf{x}|\mathbf{x}={\mathbf{x}}_i^D,{\mathbf{x}}_i^D\in \overline{\Omega },i=1,2,\dots ,M_D\}.$$

$${\Omega }_i^D=\{\mathbf{x}\in \Omega ||\mathbf{x}-{\mathbf{x}}_i^D|<|\mathbf{x}-{\mathbf{x}}_{i^{\pm }}^D|,i^{\pm }=1,2,\dots ,M_D,i^{\pm }\ne i\},i=1,2,\dots ,M_D.$$

$${\overline{\omega }}^V=\{\mathbf{x}|\mathbf{x}={\mathbf{x}}_j^V,{\mathbf{x}}_j^V\in \overline{\Omega },j=1,2,\dots ,M_V\},$$

$$\overline{\Omega }=\bigcup _{i=1}^{M_D}{\overline{\Omega }}_i^D,{\overline{\Omega }}_i^D={\Omega }_i^D\cup \partial {\Omega }_i^D,{\Omega }_i^D\cap {\Omega }_{i^{\pm }}^D=\varnothing ,i^{\pm }\ne i.$$

$$\overline{\Omega }=\bigcup _{j=1}^{M_V}{\overline{\Omega }}_j^V,{\overline{\Omega }}_j^V={\Omega }_j^V\cup \partial {\Omega }_j^V,{\Omega }_j^V\cap {\Omega }_{j^{\pm }}^V=\varnothing ,j^{\pm }\ne j.$$

$$(y,v{)}_D=\sum _{\mathbf{x}\in {\overline{\omega }}^D}y(\mathbf{x})v(\mathbf{x})S^D(\mathbf{x}),\Vert y{\Vert }_D=(y,y{)}_D^{1/2}.$$

$$S^D(\mathbf{x})=\mathrm{meas}({\Omega }_i^D),\mathbf{x}={\mathbf{x}}_i^D\in {\overline{\omega }}^D.$$

$$H_0({\omega }^D)=\{u\in H({\omega }^D)|u(\mathbf{x})=0,\mathbf{x}\in \partial {\omega }^D\},$$

$$(y,v{)}_D=\sum _{\mathbf{x}\in {\omega }^D}y(\mathbf{x})v(\mathbf{x})S^D(\mathbf{x}).$$

$$(y,v{)}_V=\sum _{\mathbf{x}\in {\omega }^V}y(\mathbf{x})v(\mathbf{x})S^V(\mathbf{x}),\Vert y{\Vert }_V=(y,y{)}_V^{1/2},$$

$$S^V(\mathbf{x})=\mathrm{meas}({\Omega }_j^V),\mathbf{x}={\mathbf{x}}_j^V\in {\omega }^V.$$

$$\mathbf{v}(\mathbf{x})={\mathbf{v}}_D(\mathbf{x})+{\mathbf{v}}_V(\mathbf{x}),{\mathbf{v}}_D(\mathbf{x})=v_D(\mathbf{x}){\mathbf{e}}_D(\mathbf{x}),{\mathbf{v}}_V(\mathbf{x})=v_V(\mathbf{x}){\mathbf{e}}_V(\mathbf{x}),{\mathbf{x}}_m^{\ast }\in {\omega }^{\ast }.$$

$$(\bm{v},\bm{w})_{*}=\sum_{\bm{x}\in\omega^{*}}\big{(}v_{D}(\bm{x})\,w_{D}(\bm{x})+v_{V}(\bm{x})\,w_{V}(\bm{x})\big{)}\,S^{*}(\bm{x}),\quad\|\bm{v}\|_{*}=(\bm{v},\bm{v})_{*}^{1/2}.$$

$$S^{\ast }(\mathbf{x})=\mathrm{meas}({\Omega }_m),\mathbf{x}={\mathbf{x}}_m^{\ast }\in {\omega }^{\ast }.$$

$$y_D(\mathbf{x})=y(\mathbf{x}),\mathbf{x}\in {\overline{\omega }}^D,$$

$$y_V(\mathbf{x})=y(\mathbf{x}),\mathbf{x}\in {\overline{\omega }}^V.$$

$$({\mathrm{grad}}_{h}y)(\mathbf{x})=({\mathrm{grad}}_D{y}_D)(\mathbf{x})+({\mathrm{grad}}_V{y}_V)(\mathbf{x}),\mathbf{x}={\mathbf{x}}_m^{\ast }\in {\omega }^{\ast }.$$

$$({\mathrm{grad}}_D{y}_D)(\mathbf{x})=({\partial }_D{y}_D)(\mathbf{x}){\mathbf{e}}_D(\mathbf{x}),({\mathrm{grad}}_V{y}_V)(\mathbf{x})=({\partial }_V{y}_V)(\mathbf{x}){\mathbf{e}}_V(\mathbf{x}),\mathbf{x}\in {\omega }^{\ast }.$$

$${\mathrm{grad}}_D{y}_D({\mathbf{x}}_m^{\ast })=\frac{y({\mathbf{x}}_{i^{+}}^D)-y({\mathbf{x}}_i^D)}{|{\mathbf{x}}_{i^{+}}^D-{\mathbf{x}}_i^D|}{\mathbf{e}}_D({\mathbf{x}}_m^{\ast }),$$

$${\mathrm{grad}}_V{y}_V({\mathbf{x}}_m^{\ast })=\frac{y({\mathbf{x}}_{j^{+}}^V)-y({\mathbf{x}}_j^V)}{|{\mathbf{x}}_{j^{+}}^V-{\mathbf{x}}_j^V|}{\mathbf{e}}_V({\mathbf{x}}_m^{\ast }).$$

$$({\partial }_D{y}_D)({\mathbf{x}}_m^{\ast })=\frac{y({\mathbf{x}}_{i^{+}}^D)-y({\mathbf{x}}_i^D)}{|{\mathbf{x}}_{i^{+}}^D-{\mathbf{x}}_i^D|},({\partial }_V{y}_V)({\mathbf{x}}_m^{\ast })=\frac{y({\mathbf{x}}_{j^{+}}^V)-y({\mathbf{x}}_j^V)}{|{\mathbf{x}}_{j^{+}}^V-{\mathbf{x}}_j^V|}.$$

$$\mathbf{H}({\omega }^{\ast })={\mathbf{H}}_D({\omega }^{\ast })+{\mathbf{H}}_V({\omega }^{\ast }),$$

$${\mathbf{H}}_D({\omega }^{\ast })=\{\mathbf{v}\in \mathbf{H}({\omega }^{\ast })\mid {\mathbf{v}}_V(\mathbf{x})=0,\mathbf{x}\in {\omega }^{\ast }\},$$

Full text

††journal: arXiv

1 Introduction

Applied mathematical models provide quantitative tools for analyzing and predicting physical, chemical, and biological processes [1]. Their theoretical description is based on systems of coupled partial differential equations, which include elliptic, parabolic, and hyperbolic types [2]. The corresponding initial–boundary value problems are formulated in multidimensional domains with complex geometry, and numerical methods are the primary tool for their solution [3, 4].

Modern computational practice makes extensive use of unstructured meshes for irregular geometries [5]. For two-dimensional problems, triangular meshes are most common, while in three dimensions tetrahedral meshes are used [6]. Spatial approximations are most often constructed using the finite element method (FEM) [7, 8]. The widely known FEniCS computational platform [9], for example, is based on FEM technology.

Finite difference methods for boundary value problems are well developed on regular meshes [10, 11]. On irregular meshes, approximate solutions are constructed using the balance (integro–interpolation) method [12]. In this case, the discrete problem is regarded as a conservation law for individual cells of the computational domain (control volumes). This approach is now widely referred to as the finite volume method (FVM) [13, 14]. A finite-element version of FVM is presented in [15].

A key element of finite-volume technology is the use of irregular meshes, where the approximate solution is defined at the nodes. For each mesh node, a control volume is constructed, and a balance relation is written, which formulates the conservation law at the discrete level. This requires the introduction of flux quantities that provide interface relations at the boundaries of control volumes. Such discrete problems can be related to the simplest approximations using mixed finite elements [16].

Finite-volume approximations are most straightforwardly built on optimal meshes: Delaunay triangulations in two dimensions and their tetrahedral analogues in three dimensions [17, 18]. In this case, the control volumes are Voronoi polygons. It is crucial that the mesh operators of vector calculus, invariant under a change of coordinates, inherit the main properties of the corresponding differential operators by satisfying the basic integral relations [19]. Such consistent approximations form the basis of mimetic discretization (MD) technology [20, 21].

In standard mesh generation, the nodes of the computational mesh are placed at the vertices of the cells. Another possibility is to define the solution at the cell centers (e.g., at their centroids). In this case, the cell itself serves as a natural control volume for finite-volume approximation. In a Delaunay triangulation, Voronoi polyhedra uniquely define control volumes if the solution is assigned to mesh nodes. Conversely, if the solution is defined at the vertices of the Voronoi diagram, then the control volumes coincide with the cells of the Delaunay mesh. This duality motivates the use of both Delaunay and Voronoi meshes in finite-volume approximations, thereby extending the class of solvable problems.

The idea of a two-mesh technology employing both Delaunay and Voronoi meshes was proposed earlier in [22, 23]. There are still relatively few works in this direction. The recent study [24] constructs operator-difference approximations on a merged rectangular mesh for scalar and vector boundary value problems. A general treatment of vector operators on merged two-dimensional Voronoi–Delaunay meshes is presented in [25]. Here the computational domain is partitioned into orthodiagonal quadrilaterals, and consistent approximations for gradient, divergence, and curl operators are developed. The applicability of this approach to mesh formulations of scalar and vector boundary value problems is briefly discussed.

In this work, we construct finite-volume approximations for a two-dimensional diffusion–reaction problem in an anisotropic medium. In this setting, the use of a single mesh (Delaunay or Voronoi) is problematic, while the merged Voronoi–Delaunay mesh offers an effective alternative. We emphasize the main features of the approach and present numerical results for a model problem, which demonstrate the efficiency of FVM on such meshes.

The structure of the paper is as follows. Section 2 formulates the two-dimensional boundary value problem for the diffusion–reaction equation with anisotropic coefficients. Section 3 describes Delaunay and Voronoi meshes, scalar and vector mesh functions, and introduces the merged quadrilateral Voronoi–Delaunay mesh. Section 4 develops approximations of the gradient and divergence operators. Section 5 presents the mesh boundary value problem and establishes an a priori stability estimate. Section 6 provides numerical results for test problems with different diffusion tensors. Section 7 concludes with a summary of the main findings.

2 Problem formulation

We consider the linear steady–state diffusion–reaction problem in a bounded domain $\Omega\subset\mathbb{R}^{2}$. The concentration of the substance $u(\bm{x})$, $\bm{x}=(x^{(1)},x^{(2)})$, satisfies the equation

[TABLE]

Equation (2.1) is supplemented with homogeneous Dirichlet boundary conditions:

[TABLE]

The medium is assumed to be anisotropic. The diffusion coefficient $K(\bm{x})=\{k_{\alpha\beta}(\bm{x})\}$ is a symmetric, positive-definite tensor of rank two:

[TABLE]

for any vector $\bm{\zeta}=(\zeta_{1},\zeta_{2})$. The reaction coefficient $r(\bm{x})$ is assumed nonnegative: $r(\bm{x})\geq 0$.

Let $\mathcal{H}=L_{2}(\Omega)$ with inner product and norm

[TABLE]

For solutions of problem (2.1), (2.2), one can derive a number of a priori estimates. In particular, the fundamental estimate

[TABLE]

holds, where $c$ is the constant from Friedrichs’ inequality. This relation will serve as a reference when analyzing discrete problems.

To construct finite-volume approximations of problem (2.1), (2.2), we introduce auxiliary vector variables, following the mixed finite element methodology [16]. Define the flux vector as

[TABLE]

With this notation, equation (2.1) can be written as

[TABLE]

In the finite-volume method, we approximate the system (2.4)–(2.5) together with the boundary condition (2.2), using discrete analogues of vector operators.

The finite-volume construction involves the introduction and subsequent elimination of flux variables. The main stages are as follows:

Assign the primary variable $u(\bm{x})$ on a mesh, seeking its approximation at the mesh nodes. These nodes define a partition of the computational domain into primary cells.

• 2.

Introduce a separate mesh for the auxiliary flux variable $\bm{q}(\bm{x})$. Each primary mesh node is associated with a control volume, and fluxes are evaluated at the boundaries of these volumes.

• 3.

For each control volume, write a discrete analogue of equation (2.5) by integrating over the control volume and incorporating primary and auxiliary mesh quantities.

• 4.

Approximate the flux variables using differential relations such as (2.4).

• 5.

Eliminate the auxiliary fluxes to obtain a discrete problem solely for the primary unknowns.

In this way, the finite volume method provides a constructive framework for developing mesh-based formulations of boundary value problems on irregular meshes. Below we consider the key elements of this approach for problem (2.1), (2.2).

3 Meshes and mesh functions

Meshes are constructed for the approximation of scalar and vector functions. The primary mesh is a Delaunay triangulation without obtuse angles (D-mesh). Based on it, we build a Voronoi diagram (V-mesh). From both meshes we introduce a merged Voronoi–Delaunay mesh (MVD-mesh) with quadrilateral cells, in which vector function values are defined at the intersections of the diagonals.

3.1 Delaunay mesh

As the primary mesh, we consider a triangular mesh in the computational domain $\Omega$, obtained by Delaunay triangulation. We assume that $\Omega$ is a convex polygon with boundary $\partial\Omega$.

The Delaunay triangulation is assumed to contain no obtuse angles. Such triangulations are often called non-obtuse meshes. A simple example is shown in Figure 1. The usual mesh quality requirement [18] is related to the minimum interior angles of the triangles. This characteristic is typically used in mesh generation algorithms and in software design. Computational algorithms for controlling the maximum angle of triangulations are described in [26].

In the closure $\overline{\Omega}=\Omega\cup\partial\Omega$, we denote the D-mesh nodes as $\bm{x}_{i}^{D}$, $i=1,2,\ldots,M_{D}$:

[TABLE]

We distinguish interior nodes $\omega^{D}$ and boundary nodes $\partial\omega^{D}$, which include the vertices of $\Omega$: $\overline{\omega}^{D}=\omega^{D}\cup\partial\omega^{D}$.

3.2 Voronoi mesh

Each node $\bm{x}_{i}^{D}$, $i=1,2,\ldots,M_{D}$, is associated with a Voronoi polygon $\Omega_{i}^{D}$, or its intersection with $\Omega$. The Voronoi polygon for node $i$ is the set of points in $\Omega$ closer to $\bm{x}_{i}^{D}$ than to any other node $\bm{x}_{i^{\pm}}^{D}$:

[TABLE]

If the primary mesh is acute, then the Voronoi vertices inside $\Omega$ always lie within the corresponding triangles. Vertices of Voronoi polygons lying on $\partial\Omega$ are the midpoints of Delaunay edges belonging to the boundary.

Let $\bm{x}_{j}^{V}$, $j=1,2,\ldots,M_{V}$, be the Voronoi vertices. For the triangulation in Figure 1, the Voronoi mesh is shown in Figure 2.

The set of all Voronoi vertices defines the V-mesh:

[TABLE]

with $\overline{\omega}^{V}=\omega^{V}\cup\partial\omega^{V}$.

In discretization, mesh cells are defined by their vertices, while in the finite volume method the domain is additionally partitioned into control volumes associated with mesh nodes.

The control volumes for the D-mesh are the regions $\Omega_{i}^{D}$:

[TABLE]

Similarly, each interior Voronoi node $\bm{x}_{j}^{V}$, $j=1,2,\ldots,M_{V}$, corresponds to a control volume $\Omega_{j}^{V}$, with

[TABLE]

For boundary nodes $\bm{x}_{j}^{V}\in\partial\Omega$, no control volumes are assigned: $\Omega_{j}^{V}=\emptyset$.

3.3 Merged Voronoi–Delaunay mesh

From the nodes of the D- and V-meshes we construct a partition of $\Omega$ into quadrilateral cells. Near the boundary these cells degenerate into triangles. The quadrilaterals are orthodiagonal, i.e., their diagonals are orthogonal. The center of each MVD-cell $\Omega_{m}$, $m=1,2,\ldots,M$, is the intersection of its diagonals. An example of the MVD-mesh, corresponding to Figures 1 and 2, is shown in Figure 3.

Figure 4 shows the three types of cells: Delaunay triangles, Voronoi polygons, and MVD orthodiagonal quadrilaterals.

3.4 Scalar mesh functions

For scalar functions $u(\bm{x})$, $\bm{x}\in\Omega$, we introduced the Hilbert space $\mathcal{H}=L_{2}(\Omega)$ with its inner product and norm. Analogously, we define Hilbert spaces of mesh functions on the D- and V-meshes. Mesh functions are defined at mesh nodes, each associated with a control volume.

For D-mesh functions at nodes $\bm{x}_{i}^{D}$, $i=1,2,\ldots,M_{D}$, we introduce the space $H(\omega^{D})$ with inner product

[TABLE]

Here $S^{D}(\bm{x})$ is the measure of the corresponding control volume:

[TABLE]

If homogeneous Dirichlet conditions (2.2) are imposed, then

[TABLE]

and the inner product reduces to

[TABLE]

Similarly, for V-mesh functions at interior nodes $\bm{x}_{j}^{V}$, $j=1,2,\ldots,M_{V}$, we define the space $H(\omega^{V})$:

[TABLE]

where

[TABLE]

3.5 Vector mesh functions

The centers of MVD-cells are denoted $\bm{x}_{m}^{*}$, $m=1,2,\ldots,M$, and the mesh is $\omega^{*}$. Vector mesh functions are defined at these nodes.

In each MVD-cell we introduce a local orthogonal coordinate system with unit vectors $\bm{e}_{D}(\bm{x})$ and $\bm{e}_{V}(\bm{x})$. The components of a vector $\bm{v}(\bm{x})=\{v_{D}(\bm{x}),v_{V}(\bm{x})\}$ are interpreted as tangential components on the faces of the Delaunay and Voronoi cells, respectively:

[TABLE]

The MVD-mesh with local coordinate directions is shown in Figure 5.

For vector mesh functions at $\bm{x}_{m}^{*}$ we define the space $\bm{H}(\omega^{*})$ with inner product

[TABLE]

Here

[TABLE]

4 Operator-difference approximations

When studying boundary value problems for second-order elliptic and parabolic equations, the main focus is on approximations of the fundamental differential operators of vector calculus: the gradient of a scalar function and the divergence of a vector function.

4.1 Mesh gradient operator

A scalar mesh function $y(\bm{x})$ is defined by its values at the nodes of the Delaunay mesh,

[TABLE]

and at the nodes of the Voronoi mesh,

[TABLE]

In the local coordinate system (see Figure 6), the mesh gradient is written as

[TABLE]

The components of the gradient vector are associated with finite differences along the faces of the D- and V-meshes. We set

[TABLE]

Thus, for the components we obtain

[TABLE]

[TABLE]

Hence,

[TABLE]

We distinguish subspaces of vector mesh functions $\bm{H}(\omega^{*})$ by

[TABLE]

where

[TABLE]

[TABLE]

Accordingly, the domains and ranges of the gradient operators are

[TABLE]

4.2 Flux vector

The flux vector is defined according to (2.4). Within each MVD-cell, the diffusion tensor must be expressed in the corresponding local coordinate system.

Let $\theta(\bm{x}_{m}^{*})$ denote the rotation angle from Cartesian to the local coordinates of cell $\Omega_{m}$. The clockwise rotation matrix is

[TABLE]

In local coordinates, the diffusion tensor is

[TABLE]

with components

[TABLE]

Taking into account (4.1), the mesh flux vector is

[TABLE]

where

[TABLE]

Thus, in the local coordinate system we obtain convenient formulas for the approximate computation not only of the gradient but also of the flux vector in anisotropic media.

4.3 Mesh divergence operator

The divergence approximation is constructed separately at Delaunay and Voronoi mesh nodes when vector functions are given at $\omega^{*}$ (the centers of MVD-cells). The approach relies on the integral equality for the control volume.

For $\bm{v}(\bm{x})$ defined at $\bm{x}\in\omega^{*}$, the approximation

[TABLE]

is evaluated at $\bm{x}\in\omega$. At Delaunay nodes we write $\operatorname{div}_{D}$, and at Voronoi nodes — $\operatorname{div}_{V}$.

Consider the case where a node of $\omega$ belongs to the Voronoi diagram, so the control volume is a Delaunay cell (Figure 7). By the divergence theorem,

[TABLE]

Thus, the operator $\operatorname{div}_{V}:\bm{H}(\omega^{*})\to H(\omega^{V})$ is defined as

[TABLE]

where $\widetilde{\bm{v}}(\bm{x})$ is the piecewise constant interpolation of $\bm{v}$ on $\partial\Omega_{j}^{V}$.

Let $\bm{n}_{j}^{V}(\bm{x})$ be the outward normal to $\partial\Omega_{j}^{V}$. Using the notations of Figure 7, from (4.5) we obtain

[TABLE]

Here $\delta_{j}^{V}(\bm{x})=\bm{e}_{V}(\bm{x})\cdot\bm{n}_{j}^{V}(\bm{x})$ equals $1$ if $\bm{e}_{V}(\bm{x})$ coincides with $\bm{n}_{j}^{V}(\bm{x})$, and $-1$ otherwise.

Similarly, at Delaunay nodes the operator $\operatorname{div}_{D}$ is constructed using expressions like

[TABLE]

where $\delta_{i}^{D}(\bm{x})=\bm{e}_{D}(\bm{x})\cdot\bm{n}_{i}^{D}(\bm{x})$, and $\bm{n}_{i}^{D}(\bm{x})$ is the outward normal to the Voronoi control volume $\Omega_{i}^{D}$.

Finally, taking into account representation (3.1), we write

[TABLE]

with

[TABLE]

Thus, $\operatorname{div}_{D}$ is defined on $\bm{H}_{D}(\omega^{*})$ with values in $H(\omega^{D})$, while $\operatorname{div}_{V}$ is defined on $\bm{H}_{V}(\omega^{*})$ with values in $H(\omega^{V})$.

5 Mesh boundary value problem

A discrete analogue of the boundary value problem (2.1), (2.2) with sufficiently smooth coefficients and right-hand side is constructed using the finite volume method. To this end, the mixed formulation (2.2), (2.4), (2.5) is employed, together with the mesh gradient and divergence operators defined on the considered grids.

5.1 Finite-volume approximation

The discrete system is derived by approximating the gradient and divergence operators in equations (2.4), (2.5). Equation (2.5) leads to separate relations for the interior nodes of the Delaunay and Voronoi meshes:

[TABLE]

[TABLE]

The approximation of (2.4) is performed according to (4.3). Substituting the mesh flux vector (4.3) into (5.1), (5.2) yields a system of equations for the approximate solution:

[TABLE]

[TABLE]

For boundary nodes, the mesh functions satisfy homogeneous Dirichlet conditions:

[TABLE]

The coupled nature of the system (5.3), (5.4) appears whenever $k_{DV}(\bm{x})=k_{VD}(\bm{x})\neq 0$ for at least some nodes of the MVD mesh, $\bm{x}\in\omega^{*}$. In the particular case of a scalar diffusion coefficient, $K(\bm{x})=k(\bm{x})$, the system reduces to

[TABLE]

[TABLE]

In this case, the approximate solutions are computed independently on the Delaunay and Voronoi meshes. If no strategy is used to combine these two solutions (e.g., to improve the overall accuracy), the use of two meshes becomes of limited practical value.

5.2 Consistency of gradient and divergence approximations

For sufficiently smooth continuous functions, the following integral identity holds:

[TABLE]

In the case $u(\bm{x})=0$ for $\bm{x}\in\partial\Omega$, this reduces to

[TABLE]

which expresses that the operator adjoint to the divergence is equal to the gradient with a negative sign, i.e. $\operatorname{div}^{*}=-\operatorname{grad}$, in the corresponding Hilbert spaces. The discrete analogue of (5.8) is tied to the construction of consistent mesh approximations of the gradient and divergence, as realized in the framework of mimetic discretization methods [19, 20, 21].

At the mesh level, the first term in (5.8) corresponds to scalar products in $H_{0}(\omega^{D})$ and $H(\omega^{V})$. To obtain a discrete analogue, we consider individual cells $\Omega_{m}$, $m=1,2,\ldots,M$, of the MVD mesh. For the Delaunay nodes, one has

[TABLE]

When computing $I_{m}$ over the cell $\Omega_{m}$ (Figure 6), two Delaunay nodes $\bm{x}_{i}^{D}$ and $\bm{x}_{i^{+}}^{D}$ are involved. In local coordinates, this gives

[TABLE]

For functions satisfying $y_{D}(\bm{x})=0$ on $\partial\omega^{D}$, this yields a mesh analogue of (5.8):

[TABLE]

A similar argument applies for the Voronoi mesh, leading to

[TABLE]

Thus, for the introduced function spaces,

[TABLE]

5.3 A priori estimate

For the mesh problem (5.3)–(5.5), the standard issues of approximation, stability, and convergence of the discrete solution to the exact one are considered. Here we restrict ourselves to establishing a mesh analogue of the a priori estimate (2.3) for the boundary value problem (2.1), (2.2).

Multiply (5.3) scalarly in $H(\omega^{D})$ by $y^{D}$, and (5.4) scalarly in $H(\omega^{V})$ by $y^{V}$. Using (5.9) and the introduced notation, we obtain

[TABLE]

[TABLE]

Adding these relations and using the properties of the coefficients in (2.1), we arrive at

[TABLE]

For mesh functions on Delaunay and Voronoi grids, discrete Friedrichs-type inequalities hold. In [27], it is shown that

[TABLE]

where the constant $c_{h}$ depends only on the computational domain $\Omega$. An analogous inequality holds for functions on the Voronoi mesh:

[TABLE]

Introducing the combined norm

[TABLE]

from (5.10) we obtain

[TABLE]

Inequality (5.11) is the mesh analogue of the continuous a priori estimate (2.3).

6 Numerical experiments

We present the results of the approximate solution of a diffusion–reaction problem in a rectangle using the finite volume method on a merged Voronoi–Delaunay mesh. The main focus is on the accuracy of the approximate solution in various norms for a problem with a known solution under changing sizes of the computational meshes.

6.1 Test problem

Below are the results of calculations for the boundary value problem (2.1), (2.2) in the rectangle

[TABLE]

We limited ourselves to the case of a homogeneous medium, where the diffusion tensor $K$ is constant. The main attention was on problems in anisotropic media, for which the developed computational technology of joint use of Delaunay and Voronoi meshes is most interesting.

In the calculations, the constant diffusion tensor $K$ was set to

[TABLE]

Thus, problems with strong anisotropy of the medium were considered. The right-hand side of (2.1) corresponded to the exact solution

[TABLE]

The accuracy of the approximate solution $y(\bm{x})$ was evaluated separately on the Delaunay mesh and on the Voronoi mesh. In the $L_{2}$ norm, we have

[TABLE]

and in the $L_{\infty}$ norm,

[TABLE]

The accuracy of gradient computation at the nodes of the MVD-mesh was also evaluated:

[TABLE]

6.2 Meshes

A sequence of successively refined meshes was used for numerical experiments. Triangular meshes were generated using the gmsh mesh generator [28], and based on them, Voronoi and MVD meshes were constructed. No special efforts were made to generate meshes with guaranteed non-obtuse angles.

We used 22 meshes in total. The starting (coarsest) meshes are shown in Figs. 1–3. The number of nodes for three meshes is given in Table 1. The quality of the initial Delaunay triangulation can be assessed from the data on the minimum and maximum angles of the triangular mesh cells.

6.3 Accuracy of the approximate solution

The error norms of the approximate solution for different diffusion tensors are shown in Figs. 8–12. The main conclusions can be formulated as follows:

The error of the approximate solution in the $L_{2}$ norm decreases approximately linearly with the number of nodes on both the Delaunay mesh and the Voronoi mesh (Figs. 8, 9).

• 2.

In the uniform $L_{\infty}$ norm, the error decreases, preserving the same asymptotic dependence on the number of mesh nodes (Figs. 10, 11).

• 3.

The error of gradient computation remains sufficiently small (Fig. 12).

• 4.

The dependence on medium anisotropy is pronounced. We observe a significant drop in accuracy when solving the diffusion problem in an anisotropic medium compared to the isotropic case. This is most evident in the approximate computation of the gradient.

7 Conclusions

For a Delaunay triangulation with acute-angled cells, a Voronoi diagram is constructed. In computational practice, the mesh solution is usually determined at Delaunay nodes in two variants of domain partitioning: (a) node-centered approximation using Delaunay triangles, and (b) cell-centered approximation using Voronoi polygons. Merged Voronoi–Delaunay (MVD) meshes are introduced, whose nodes include both Delaunay vertices and Voronoi vertices. The cells of such a mesh are orthodiagonal quadrilaterals. 2. 2.

In the finite volume method, both the solution and fluxes are approximated. The approximate solution is obtained by integrating the equation over a control volume. In the MVD approach, for Delaunay nodes the control volumes are Voronoi cells, while for Voronoi nodes the control volumes are Delaunay cells. Fluxes are approximated at the centers of MVD cells, i.e. at the intersections of the orthogonal diagonals. 3. 3.

On irregular computational meshes, accurate discretization requires approximations of invariant vector (tensor) analysis operators. Discrete analogues of gradient and divergence operators have been constructed on MVD meshes in a form consistent with the Mimetic Discretization methodology. 4. 4.

The developed computational technology has been applied to a two-dimensional diffusion–reaction problem in anisotropic media. The diffusion coefficient is represented by a symmetric tensor of rank two. A coupled system of equations is formulated for the mesh solution at both Delaunay and Voronoi nodes. Numerical experiments confirm the accuracy, robustness, and efficiency of the proposed approach.