Why do we need Voronoi cells and Delaunay meshes? Essential properties of the Voronoi finite volume method

TL;DR

This paper highlights the importance of Voronoi finite volume methods and Delaunay meshes in accurately preserving the stability and qualitative properties of solutions to parabolic and elliptic problems across various resolutions and time scales.

Contribution

It provides a concise overview of the key properties of Voronoi FVM, demonstrates their advantages, and advocates for their increased adoption in practical applications.

Findings

Voronoi FVM preserves essential stability properties.

Delaunay meshes accurately approximate problem geometry.

Voronoi FVM maintains qualitative solution properties across resolutions.

Abstract

Unlike other schemes that locally violate the essential stability properties of the analytic parabolic and elliptic problems, Voronoi finite volume methods (FVM) and boundary conforming Delaunay meshes provide good approximation of the geometry of a problem and are able to preserve the essential qualitative properties of the solution for any given resolution in space and time as well as changes in time scales of multiple orders of magnitude. This work provides a brief description of the essential and useful properties of the Voronoi FVM, application examples, and a motivation why Voronoi FVM deserve to be used more often in practice than they are currently.

Equations

$$\forall {\mathbf{v}}_i,{\mathbf{E}}_j^d\in {\mathcal{T}}_h(\Omega ),\text{leavevmode}{\mathbf{v}}_i\notin {\mathrm{B}}^{\mathrm{o}}({\mathbf{E}}_j^d).$$

$$V_i:=\{\mathbf{x}\in {\mathbb{R}}^d:\Vert \mathbf{x}-{\mathbf{v}}_i\Vert <\Vert \mathbf{x}-{\mathbf{v}}_j\Vert \text{leavevmode}\forall {\mathbf{v}}_j\in {\mathcal{T}}_h(\Omega )\}.$$

$$\forall {\mathbf{v}}_i\in {\mathcal{T}}_h(\Omega ),\forall {\mathbf{E}}_j^m\in \partial {\Omega }_i,\text{leavevmode}{\mathbf{v}}_i\notin {\mathrm{B}}^{\mathrm{o}}({\mathbf{E}}_j^m)(1\le m<d).$$

$${\tilde{G}}_{i+1}=\left(\begin{array}{cc}[l]{\tilde{G}}_i& 0_{\frac{i(i-1)}{2}}\\ -{\mathbb{I}}_{i+1}& 1_{i+1}\end{array}\right),i=1,2,\dots ,$$

$${\tilde{G}}_1\left(\begin{array}{c}{\mathbf{x}}_1\\ {\mathbf{x}}_2\end{array}\right)=\left(\begin{array}{cc}-1& 1\end{array}\right)\left(\begin{array}{c}{\mathbf{v}}_1\\ {\mathbf{v}}_2\end{array}\right)={\mathbf{v}}_2-{\mathbf{v}}_1.$$

$${\tilde{G}}_2$$

$${\tilde{G}}_3$$

$${\tilde{G}}_4$$

$${\left[|{\tilde{G}}_d{\mathbf{E}}_i^d|\right]}^{-1}{\tilde{G}}_d\mathbf{u}$$

$${\tilde{G}}_d^T[{\sigma }_{{\mathbf{E}}_i^d}]={\tilde{G}}_d^{\mathrm{T}}{\left[|{\tilde{G}}_d{\mathbf{E}}_i^d|\right]}^{-1}{\left\{\left[|{\tilde{G}}_d{\mathbf{E}}_i^d|\right]\left[{\sigma }_{{\mathbf{E}}_i^d}\right]\right\}}^{\frac{1}{2}}$$

$$\epsilon {\int }_{{\mathbf{E}}_i^d}\nabla \cdot \nabla u\mathrm{d}V={\oint }_{\Sigma }\nabla u\mathrm{d}\mathbf{S}\approx \epsilon {\tilde{G}}_d^{\mathrm{T}}[{\sigma }_i]{\left[|{\tilde{G}}_d{\mathbf{E}}_i^d|\right]}^{-1}{\tilde{G}}_d\mathbf{u}:=\epsilon G_d^{\mathrm{T}}{G}_d\mathbf{u}$$

$$\alpha (u)u+\beta (u)\frac{\partial u}{\partial n}+\gamma (u)=0$$

$$\epsilon \frac{\partial u}{\partial n}\partial V_{s,i}=-\partial V_{s,i}\frac{\epsilon }{\beta (u)}\left(u_i\alpha (u_i)+\gamma (u_i)\right).$$

$$P^{-\mathrm{T}}{P}^{-1},P={\mathbf{v}}_j-{\mathbf{v}}_0,j=1,\dots ,d.$$

$$A=G_d^{\mathrm{T}}[\epsilon ]G_d$$

$$u_1=1-\epsilon ,u_j={\epsilon }_j,\text{with}\text{leavevmode}0<{\epsilon }_j\le \epsilon \ll 1,j=i,2,\dots ,k,$$

$$\{\begin{array}{ll}{W}_{ij}=\frac{{\sigma }_{ij}}{l_{ij}}=1,& j=2,\dots ,k,\\ W_{i1}=\frac{{\sigma }_{i1}}{l_{i1}}=-1,& \end{array}$$

$$A_{\text{non-Del}}\mathbf{u}{|}_i=\sum _{j=2}^k(u_i-u_j)-(u_i-u_1)=\sum _{j=2}^k({\epsilon }_i-{\epsilon }_j)-({\epsilon }_i-(1-\epsilon ))\approx 1>0,$$

$$A_{\text{Del}}\mathbf{u}{|}_i=\sum _{j=2}^k(u_i-u_j)+(u_i-u_1)=\sum _{j=2}^k({\epsilon }_i-{\epsilon }_j)+({\epsilon }_i-(1-\epsilon ))\approx -1<0.$$

$$\delta \tilde{u}{|}_i=(-A_{\text{Del}}^{-1}f){|}_i>0$$

$$\{\begin{array}{ll}u(y{)}^{\prime \prime }+cu(y)=0,& y\in (0,1),c=const,\\ u(0)=1,& \\ u(1)=0.& \end{array}$$

$$M_{ii}=\frac{2|{\omega }_i|}{(d+1)(d+2)},$$

$$M_{11}=M_{22}=\frac{2\times \left(6\times \frac{h^2}{2}\right)}{12}=\frac{1}{2}{h}^2.$$

$$M_{11}=\frac{2\times \left(8\times \frac{h^2}{2}\right)}{12}=\frac{2}{3}{h}^2,\text{whereas}{M}_{22}=\frac{2\times \left(4\times \frac{h^2}{2}\right)}{12}=\frac{1}{3}{h}^2,$$

$$1>1-\epsilon >u_b>\epsilon >0,u_w=1-u_b.$$

$$\{\begin{array}{l}\frac{\partial u}{\partial t}-\nabla \cdot \left(f(|\nabla v|)\left(\nabla u+\frac{\nabla w}{{\varphi }^{\prime \prime }(u)}\right)\right)=0\text{in Ω},\\ u(0,\cdot )=u_0(\cdot ),\\ v={\varphi }^{\prime }(u)+w,\\ w(t,\mathbf{x})={\int }_{\Omega }\mathcal{K}\left(|\mathbf{x}-\mathbf{y}|\right)\left(1-2u(t,\mathbf{y})\right)\mathrm{d}\mathbf{y}.\end{array}$$

$$F_{gl}(u)={\int }_{\Omega }\left(\varphi (u)+\kappa u\left(1-u\right)+\frac{\lambda }{2}|\nabla u{|}^2\right)\mathrm{d}\mathbf{x}.$$

$$F(u)={\int }_{\Omega }\left(\varphi (u)+u{\int }_{\Omega }\mathcal{K}\left(|x-y|\right)\left(1-u(y)\right)\mathrm{d}y\right)\mathrm{d}\mathbf{x}.$$

$$\{\begin{array}{ll}-{\sigma }^2\nabla \cdot \nabla w+w=m\left(1-2u\right)& \text{in Ω},\\ \nu \cdot \nabla w=0& \text{on ∂Ω}.\end{array}$$

$$\left(G^T{\sigma }^{2}G+[V]\right){\mathbf{w}}^{+}$$

$$[V]\frac{\mathbf{u}-{\mathbf{u}}^{+}}{\tau }$$

$$u_0(\mathbf{x})=\epsilon +\frac{x_1}{128\left(1+2\epsilon \right)},\Omega =(0,128)\times (0,128),m=8,\epsilon =10^{-7},\sigma =\sqrt{f}=2.$$

Full text

© 2019 The Authors. Preprint of the Work accepted for publication in CMMP by Pleiades Publishing (http://pleiades.online).

\usetkzobj

all

1 Introduction

Many application problems result in highly nonlinear elliptic or parabolic partial differential equations. Many of them have been studied for a long time and analytic results are available, which are often based on a variational formulation. Using entropy or free energy functionals allows the derivation of solution properties and decay properties of the free energy without restrictive smallness assumption on the free energy or the boundary conditions. Typically, a priori bounds for the solutions, such as positivity and boundedness, have to be derived and decay of the free energy along trajectories can be shown.

In the discrete case, it is less obvious. Often the discretization parameters in space and time, $h$ and $\tau$, which are in fact limits, are regarded as finite and universally achievable. However, this is not always the case in situations with limited smoothness, such as in states further away from the thermodynamic equilibrium or in case of higher free energy densities, and particularly not, if the approximation of the non-smooth limits is required. Solutions may have plateaus and internal or boundary layers connecting the plateaus and the boundary values. Energy densities are often very large if compared to a typical human environment. A well-known example is the charge transport in semiconductors, where the large difference in force results from the fact that the difference of the electrostatic potential of merely $1\text{\,}\mathrm{V}$ corresponds to a height difference of $4\;000\text{\,}\mathrm{km}$ in the gravitational field of the earth (mechanical heat equivalent relation at $300\text{\,}\mathrm{K}$); accordingly, $100\text{\,}\mathrm{V}$, which is not much for such problems, correspond to a height difference from the earth to somewhere behind the moon. Many questions in electro-chemistry, chemical reactions, and other areas result in the same type of problems and require the preservation of qualitative properties in their discrete form.

In such cases, the classic discretization error has only a limited significance with respect to the guarantee of a stable discreet solution. It is often even physically irrelevant whether a stationary limit value at the smallest time scale of $1\text{\,}\mathrm{s}$ of the problem is reached at $10^{19}$ or $10^{20}$, since the relevant human response time and life span are $1\text{\,}\mathrm{s}$ and $100\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$, respectively (100\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$\approx 100\times 400\times 25\times$3\;600\text{\,}\mathrm{s}$=$3.6\text{\times}{10}^{9}\text{\,}\mathrm{s}). In contrast, of great interest are stationary states and qualitative statements, such as whether design A provides a better solution in the $\mathrm{GHz}$ range then design B. Therefore, preservation of the structural properties of the analytical problems in space and time for finite and not too small discretization parameters is an essential.

Moreover, often, the provided data is not precise and internal layers are separating volumes with different properties (generation on one side of the layer, recombination on the other). To estimate these volumes within the order of percent of the total volume may be sufficient. What is often needed is the position of the layer but not its resolution. As long as the discretization scheme inherits the stability properties of the analytic problem independently of the spatial and time step sizes, the requirement to resolve everything may be not adequate.

It is not a priori clear which combination of space-time discretization will provide unconditionally stable discrete equations sufficiently well approximating the original problem or whether such a scheme is efficient. Voronoi-Delaunay meshes and implicit Euler time discretization play a special role here.

Stable discretizations for reduced problems are dating back to 1950s [3, 19, 29] and are still relevant today in a generalized setting. Voronoi cells [36] and their dual, Delaunay triangulation [6], have been used as early as [27] (e.g., the classical book by Varga [35]). Delaunay triangulations are optimal for a variety of criteria, e.g., maximizing the smallest angle [25], minimizing the roughness of the mesh [28] (the squared $L^{2}$ norm of the gradient of the piecewise linear interpolating surface) in two dimensions, and other (see, e.g., [31, Sect. 4] and the references therein). It is hard to follow the history of the development of Delaunay meshes in detail but the central points have been related to $C^{1}$ surface interpolation and to applications such as neutron transport and charge transport in semiconductors. Admissible meshes [7] summarize the minimal mesh attributes required to prove convergence.

Vienna may have been the birthplace of the term boundary conforming Delaunay meshes (e.g., comp. [8]), which are required in case of interfaces with jumps in material properties (e.g., diffusion coefficients), or nonlinear boundary conditions of the third kind, or some other special cases to ensure the boundary and internal interfaces of the problem domain to be a part of the Voronoi / Delaunay mesh, whereas the constraint Delaunay property [5], which gained popularity recently [30, 32], is not really helpful. Local orthogonal coordinates would be the best choice and, although they will often not exist globally, locally they might be the key to combine stability, good approximation properties, and low complexity of the discrete problem.

The following list summarizes some desired properties of a discrete problem in relation to its analytic counterpart:

• •

the same qualitative stability properties as the original, analytic problem,

• •

no artificial smallness or smoothness assumptions,

• •

a weak form of the discrete problem, known analytic test functions and techniques,

• •

hence, a weak discrete energy functional,

• •

some weak convergence for distributions on the right-hand side,

• •

internal layer positions as bounded jumps (unresolved layers) in the solution.

Some examples, where the list is partially fulfilled, are: transport of electrons and holes in semiconductors and phase separation (non-Cahn-Hilliard) [10, 11, 14, 15, 17, 12].

2 Delaunay meshes, Voronoi volumes, and boundary conforming Delaunay

In the following, we consider a discretization (mesh) $\mathcal{T}_{h}(\Omega)$ of a domain $\Omega\subset\mathbb{R}^{d}$ and assume that $\Omega$ is a union of finite $d$-dimensional subdomains, $\Omega=\cup_{i}\Omega_{i}$. Each subdomain contains only one material and consists of $d$-dimensional simplices $\bm{E}^{d}_{j}$, $\Omega_{i}=\cup_{j}\bm{E}^{d}_{j}$. The vertices of the mesh simplices are denoted by $\bm{v}_{k}$. Consequently, all interfaces between materials and boundaries of the problem domain are represented by low-dimensional simplices $\bm{E}^{m}_{k}$, $1\leq m<d$. It is assumed that the domain geometry representation contains all critical lines and points.

2.1 Delaunay triangulations

A discretization (mesh) $\mathcal{T}_{h}(\Omega)$ of a given $d$-dimensional domain $\Omega$ by simplices is called a Delaunay mesh or Delaunay triangulation if it fulfills the empty circumball property, that is, if no mesh vertex is inside the circumball of any simplex, i.e.,

[TABLE]

For a given mesh, a mesh simplex is called Delaunay if its circumball is empty (Fig. 1). A mesh edge is called Delaunay if there exist an empty sphere containing this edge but no other mesh vertices inside. Consequently, a simplex is a Delaunay simplex iff all of its edges are Delaunay. Note that the Delaunay edges are not necessarily the shortest edges possible (see Fig. 1, right) but a Delaunay triangulation maximises the minimum angle of a triangulation [25] and it is a minimal roughness triangulation in two dimensions [28].

A Delaunay triangulation is unique if no $d+2$ vertices lie on a common empty $d$-sphere.

2.2 Voronoi volumes, surfaces, and cells

The Voronoi volume $V_{i}$ of a vertex $\bm{v}_{i}$ is the region of $\mathbb{R}^{d}$ consisting of all points closer to $\bm{v}_{i}$ than to any other mesh vertex,

[TABLE]

Its boundary $\partial V_{i}$ is the corresponding Voronoi surface (Fig. 2, left).

The Voronoi decomposition for a set of vertices is dual to its Delaunay triangulation. Each Voronoi facet corresponds to an edge of the Delaunay triangulation.

The intersection of a Voronoi volume $V_{i}$ with a simplex $\bm{E}^{d}_{j}$, $V_{ij}=V_{i}\cap\bm{E}^{d}_{j}$, is the Voronoi volume element (cell) of $\bm{v}_{i}$ with respect to $\bm{E}^{d}_{j}$ (Fig. 2, right). Being an intersection of half-spaces, Voronoi cells and facets are always convex and their measures are continuous functions of vertex positions.

2.3 Boundary conforming Delaunay triangulation

Introducing boundary or interface conditions into the mesh might violate the Delaunay property. For instance, consider a homogeneous Neumann boundary or interface condition in Fig. 2 along the edges $\bm{v}_{1}\bm{v}_{3}$ and $\bm{v}_{3}\bm{v}_{2}$ (Fig. 3 (left)). In cell centered schemes, a common technique to treat the Neumann boundary is to introduce a ghost vertex $\bm{v}_{4}^{*}$ as a reflection of $\bm{v}_{4}$ over $\bm{v}_{1}\bm{v}_{3}$. If the ghost vertex $\bm{v}_{4}^{*}$ is inside the circumball of $\triangle\bm{v}_{1}\bm{v}_{3}\bm{v}_{4}$ (which will be the case iff $\bm{v}_{4}$ is inside the circumball of $\bm{v}_{1}\bm{v}_{3}$) then $\triangle\bm{v}_{1}\bm{v}_{4}\bm{v}_{3}$ becomes non-Delaunay. The problem extends to higher dimensions as well and motivates the following definition.

A Delaunay mesh is called boundary conforming Delaunay if no mesh vertex is inside the smallest (diametrical) circumball of any low-dimensional boundary or interface simplex, i.e.,

[TABLE]

The property of an empty diametrical ball is sometimes referred to as the Gabriel property [9]. In brief, a boundary conforming Delaunay mesh is a mesh, where each simplex is Delaunay and each boundary simplex is Gabriel.

In the situation of Fig. 3 (left), boundary conformity can be enforced by splitting $\triangle\bm{v}_{1}\bm{v}_{3}\bm{v}_{4}$ with a new vertex on the boundary edge $\bm{v}_{1}\bm{v}_{3}$. Taking the intersection of the diametrical disks of $\bm{v}_{1}\bm{v}_{4}$ and $\bm{v}_{4}\bm{v}_{3}$ with the boundary edge $\bm{v}_{1}\bm{v}_{3}$, which is exactly the projection of $\bm{v}_{4}$ onto $\bm{v}_{1}\bm{v}_{3}$ (Fig. 3, center), guarantees that all other mesh simplices remain Delaunay while restoring the boundary conformity of the mesh (Fig. 3, right).

Details on the boundary conforming Delaunay mesh generation can be found in [33] and the references therein.

3 Voronoi finite volumes and their essential properties

If applying the Gauss theorem to a Voronoi cell or its intersection with $\partial\Omega$, each facet is related to an edge and therefore is a set of simplices. It is more convenient to assemble the finite volume equations using simplices of the dual Delaunay mesh, in particular it is more convenient to define material properties per simplex than per Voronoi volume.

3.1 Voronoi Finite Volume Equations

The adjacency matrices $\tilde{G}_{i}$ of $i$-dimensional simplices mapping vertices to edges are defined recursively through the relation

[TABLE]

with $\mathbb{I}_{j}$ denoting the $j\times j$ identity matrix, $\bm{0}_{j}=(\overbrace{0,\dotsc,0}^{j}{)}^{\mathrm{T}}$, $\bm{1}_{j}=(\overbrace{1,\dotsc,1}^{j}{)}^{\mathrm{T}}$, and $\tilde{G}_{1}=(-1,1)$ describing the difference along an interval,

[TABLE]

It holds,

[TABLE]

$\tilde{G}_{d}$ maps vertices to edges. Correspondingly, $\tilde{G}_{d}^{\mathrm{T}}$ maps edges to vertices and $\tilde{G}_{d}^{\mathrm{T}}\tilde{G}_{d}$ is the graph-Laplacian of the $d$-dimensional simplex.

Denote the length of an edge $\bm{E}_{i}^{1}$ by $\lvert\tilde{G}_{1}\bm{E}{{}_{i}^{1}}\rvert$ and the diagonal matrix of all edge lengths for a simplex $\bm{E}^{d}_{i}$ by $[\,\lvert\tilde{G}_{d}\bm{E}^{d}_{i}\rvert\,]$. Then,

[TABLE]

is the approximation of $\nabla u$ in directions $\partial u/\partial n_{k}$ orthogonal to the surfaces of the Voronoi cell related to the simplex $\bm{E}^{d}_{i}$, transforming linear functions defined on each edge into the correct constant. Further,

[TABLE]

is the symmetric form of the edge to vertex map, including a geometric weight of the Voronoi facet intersected by the edge times the edge length. Hence, the metric factor $\{\cdot\}$ has the meaning of the determinant times two (the edge length is twice the hight of the $d$-dimensional simplex formed by the intersection of the simplex and $\sigma_{i}$) and not that of the volume intersected with the corresponding part of each edge of the simplex $\bm{E}^{d}_{i}$. It has to be handled in all terms of the operator, boundary conditions consistently. Hence,

[TABLE]

is the approximation of the volume integrated Laplace operator on $\bm{E}^{d}_{i}$, where $\varepsilon=const$ on $\bm{E}^{d}_{i}$ and $\sigma_{j}\geq 0$ for each edge $j$.

The boundary condition

[TABLE]

is integrated over Voronoi surfaces belonging to surface vertices and the surface is defining the normal direction and coincides with the edge direction in case of a boundary conforming mesh,

[TABLE]

The normal derivative causes a change of the corresponding diagonal element and the right-hand side. Dirichlet values follow for $\beta\approx\epsilon_{\mathrm{machine}}^{2}$ and are easily homogenized.

A boundary conforming Delaunay mesh allows to split the Voronoi cells along the boundary, which coincides with the simplices, and to introduce the same surface integration on the $d-1$-dimensional Voronoi cells. In case of boundary layers, the vertices on the boundary can assume the extrema enforced, e.g., by the boundary condition. This is not the case if the Voronoi cells are used to approximate $\Omega_{i}^{d}$.

The essential assumption for the FVM is the validity of the Gauss theorem. The largest function space where the Gauss theorem hods is the space of functions with bounded variation. It is a very weak assumption.

Allowing an arbitrary jump in the diffusion coefficient of adjacent regions results in additional angle restrictions.

3.2 Compensation and the S-matrix property

A major advantage of the Voronoi FVM is that the system matrix is a non-singular M-matrix. An M-matrix (Minkowski matrix) is a Z-matrix (off-diagonal entries are non-positive) with a positive diagonal whose column sums are all non-negative. If at least one column sum of an M-matrix is positive, then it is inverse-positive, i.e., it is invertible and its inverse is a positive matrix. A symmetric M-matrix is called an S-matrix (Stieltjes Matrix).

Consider the $\mathbb{R}^{3}$ case and two possible Delaunay tetrahedra in Figs. 4 and 4. The free fourth vertex opposing each triangular facet determines the four Delaunay conditions. In Fig. 4 (left) all circumcenters (hereafter $C$ with an index describing the corresponding, possibly low-dimensional simplex) are inside the corresponding defining simplex of dimension one, two, or three, respectively.

In Fig. 4 (right) the tetrahedron’s circumcenter $C_{1234}$ and the circumcenter $C_{234}$ of $\triangle 234$ are outside of their defining simplex. Each of the Voronoi faces is defined by two triangles, e.g., the one defined by the vertices 1 and 2 consists of $\triangle C_{12}C_{123}C_{1234}$ and $\triangle C_{21}C_{124}C_{1234}$ which are orthogonal to the edge $12$ and change the sign accordingly to the positions of the related triangle’s and tetrahedron’s circumcenters and have a positive projection on edges 12 and 21 for interior circumcenters. The neighboring tetrahedron sharing $\triangle 124$ as a facet has the same circumcenters of the shared edges and the same direction of the normal to $\triangle 124$. The Voronoi facet is bounded by a polygon connecting all circumcenters of tetrahedra sharing one edge and is located in the plane orthogonal to that edge. Summing up all signed contributions of tetrahedra sharing an edge yields the signed Voronoi facet area related to this edge. The Delaunay condition guaranties that the distance between the circumcenters of two neighboring tetrahedra is non-negative and, thus, the signed Voronoi facet area is non-negative.

This effect is called compensation: the negative contribution is compensated by the positive contribution of the neighboring tetrahedron. This relation holds in any dimension and guarantees that the finite volume matrix is an S-matrix or an M-matrix.

If more than $m\geq d+2$ vertices lie on a common $d$-sphere, $m-d$ circumcenters will coincide. A boundary conforming Delaunay mesh guarantees that all circumcenters are inside the corresponding defining simplices, the per simplex defined oriented parts of Voronoi cells and surfaces are non-negative, and the volume of $\Omega_{i}$ as well as the surface measures are preserved. This allows an easy handling of all third kind boundary conditions.

For the P1 (linear) finite elements method (FEM), the Laplace matrix in $d$ dimensions is the projection of each of the $d-1$-dimensional facets on all others,

[TABLE]

In two dimensions, it is identical to the Voronoi finite volume matrix $G_{d}^{\mathrm{T}}G_{d}$. In more than two dimensions, it is not the case anymore (except in very special cases), the compensation of negative contributions works differently, and boundary conforming Delaunay meshes are neither sufficient nor necessary for the S-matrix property for the FEM (for details see [4, 21, 37, 24] and the references therein). The so-called non-obtuse angle condition (all angles between the surface normals are non-obtuse), which guarantees the S-matrix property for the FE Laplace matrix, is too restrictive and difficult to fulfill in practice.

3.3 Maximum principle and dissipativity

The finite volume solution satisfies the weak discrete maximum principle. Let $\bm{u}$ fulfill $G_{d}^{d}G_{d}\bm{u}=\bm{0}$ with some Dirichlet boundary values $\bm{u}_{D}$, $\hat{u}=\max u_{D}$, $\check{u}=\min u_{D}$, and assume that $\max u=U>\hat{u}$. Testing the equation with $(\bm{u}-\hat{u})^{\mathrm{+T}}$ results in $(\bm{u}-\hat{u})^{\mathrm{+T}}G^{\mathrm{T}}_{d}G_{d}\bm{u}>0$ because $G_{d}\bm{u}$ and $G_{d}{(\bm{u}-\hat{u})}^{+}$ have either the same signs (and at least one edge value is non zero) or $G_{n}{(\bm{u}-\hat{u})}^{+}=\bm{0}$, hence a contradiction.

Dissipativity follows from the same argument. In the discrete case both are equivalent in this sense.

Dissipativity of zero-order terms: the Voronoi volumes form a diagonal matrix $[V]$, larger support can not be handled without smallness assumptions of the variation of the solution on the support.

3.4 The effect of non-Delaunay meshes for the Laplace operator

Since Voronoi cells have non-negative surface measures, assembling

[TABLE]

results in an S-matrix: $A_{ii}>0$, $A_{ij}\leq 0$, $i\neq j$, and ${\bm{1}}^{\mathrm{T}}_{d}G_{d}^{\mathrm{T}}=\bm{0}_{d}$ (column and row sum zero at vertices without boundary conditions). Clearly, $A^{-1}$ is positive for irreducible $A$ or connected meshes. Bounded $\varepsilon$ and boundary conditions add only positive values to the diagonal at least at one vertex. Hence, the solution is bounded by the Dirichlet values and positive for zero Dirichlet values and positive right-hand sides.

If the underlying mesh is non-Delaunay, summing up the surface contributions per edge may result in negative surfaces around this edge. It could be an arbitrary edge between two vertices not corresponding to a Voronoi face (the circumsphere of one simplex is containing the circumcenter of another). It follows $A_{ij}>0$ and positivity is not guarantied anymore, at least one may observe some local extrema at vertices $i$ and $j$ for significant deviations from the Delaunay condition. Due to the construction, $G_{d}\bm{1}_{d}=\bm{0}_{d}$. Setting $A_{ij}=0$ consistently with $A_{\mathrm{ii}}$ corresponds to introducing a local Neumann boundary condition and an inconsistent discretization.

Assume that a solution $u$ has a moving internal layer with

[TABLE]

representing the order of the variation of $u$ outside of the jump (Fig. 5). Let $A_{\text{non-Del}}$ be of type $\tilde{G}^{\mathrm{T}}[W]\tilde{G}$ and the sum over all simplices at vertex $i$ (the sum is not indicated, the summation order over edges or simplices of given type can be chosen for each purpose and the context should distinguish between a local and a global meaning).

Let $l_{ij}$ be the edge length for all neighbors $\bm{v}_{j}$ of $\bm{v}_{i}$ and

[TABLE]

i.e., the edge $\bm{v}_{i}\bm{v}_{1}$ is non-Delaunay. Hence, $A_{\text{non-Del}}{\bm{u}}|_{i}$ at vertex $i$ results in

[TABLE]

whereas in the Delaunay case

[TABLE]

Hence, the right-hand side at $i$ is very large compared to that at vertices inside the plateau of level $\varepsilon$ or $1-\varepsilon$. Because $A\bm{u}$ in Newton’s method (or the corresponding linear residual correction method) is a part of the function $f$ and $A_{\text{Del}}$ contributes to the Jacobian ($A$) we have

[TABLE]

and therefore $u^{+}=u+\delta u$, the next approximation of $u_{i}$ in the Newton process, is increased towards the mean value of its next neighbors by the diffusion part of the Jacobian.

The factor $-1$ and the huge (by modulus) right-hand side is equivalent to a large source of the wrong sign in the non-Delaunay case. Additionally, $A_{\text{non-Del}}$ can be indefinite. Clearly, a constant $A$ cannot be the reason for the jump but other terms may create jumps in case of a small diffusion. A small $A_{\text{non-Del}}$ can be sufficient to violate the bounds for $u$ and, thus, causes serious problems which can not be cured by simple limiters for $u$ without creating new issues depending on the details of the problem.

In a more general setting, this example corresponds to a locally negative dissipation rate and therefore to the local free energy production by the discretization scheme. The concept of a decaying free energy along trajectories will be violated in general and the possible steady state will deviate in dependence of the mesh from the minimal energy solution because the free energy produced by the numerical scheme has to be dissipated somewhere else.

3.5 Continuous dependence on vertex positions

Consider the following one-dimensional boundary value problem:

[TABLE]

This problem and its discretization on a given one-dimensional mesh in $y$ can be extended orthogonally in $x$ direction with arbitrary step sizes and homogeneous Neumann boundary conditions while keeping the profile of $u$ in $y$ direction (the two-dimensional solution will be invariant in $x$ direction).

The Delaunay triangulation for the resulting regular rectangular vertex set is not unique because each four neighboring vertices sharing one interval in $x$ and $y$ directions lie on one circle (Fig. 6, left: both diagonal edge directions are possible). Therefore, in this case, the unique assignment of the diagonal edges can be done only by a special rule, e.g., ‘lower left to upper right’. Using exact arithmetics to improve the Delaunay test will not help, since it will still result in non-uniqueness in case of exactly represented coordinates or randomly assign an edge depending on the rounding error in coordinates if the four vertices are co-circular by construction but have rounded coordinates in the floating point representation. For instance, rotating meshes in Fig. 6 by $\pi/4$ will result in meshes where four co-circular vertices will not be co-circular anymore in exact arithmetics, because $\sin\pi/4=\cos\pi/4=\sqrt{2}/2$ cannot be represented exactly by floating point numbers.

For the FVM, the choice of the diagonal edges is not important because the Voronoi facets corresponding to the diagonal edges have a measure of zero. The Voronoi surfaces and the Voronoi volumes are scaled correspondingly to the step size in $x$ direction, thus, preserving the invariance of the solution in $x$ direction. Moreover, small perturbations of the vertex coordinates will result in a comparably small change of facet measures (recall that Voronoi cell and facet measures are continuous functions of vertex positions). The same is true for the matrices of type $G^{\mathrm{T}}G$ in any dimension $d$ (e.g., the Laplace and the mass parts). Hence, a very small change of vertex positions will result in a very small change of the numerical solution.

For the FEM, the situation is different since topology comes into play as well. On one hand, for meshes in Fig. 6, the Laplace part of the P1 finite element discretization is invariant to the choice of the diagonal edges, since the finite element Laplace part is equal to that of the finite volumes for $d=2$ and the latter depends only on the vertex positions. On the other hand, the mass matrix $M$ depends on topology, e.g.,

[TABLE]

where $\omega_{i}$ is the patch associated with the vertex $\bm{v}_{i}$ and $\lvert\omega_{i}\rvert$ is its volume (area). In Fig. 6 (left), the mesh is translation-invariant (up to the first and the last vertex) and each patch has the same size. Hence,

[TABLE]

In Fig. 6 (right), the mesh is still Delaunay but some of the diagonal edges are chosen differently. The patches of $\bm{v}_{1}$ and $\bm{v}_{2}$ are of different sizes and so are the corresponding entries of the mass matrix,

[TABLE]

so that the mass part is not translation-invariant, while the Laplace part is. As a consequence, the numerical solution will loose its invariance in the $x$ direction. Moreover, the inconsistency of the mass matrix may also lead to lower than expected order of convergence. Lumped-mass quadrature cannot correct this issue in general because the mass assignment to vertices depends on the mesh topology and, thus, is already inconsistent before mass lumping is carried out. An interested reader can find a discussion on the lost of accuracy of linear and higher-order finite elements in dependence of the mesh topology in recent work of Kopteva [22, 23].

4 Application examples

4.1 Phase separation

Consider diffusing black and white particles which sum up to 1 everywhere (Fermi statistics), a short range interaction potential (black attracts black), and the diffusion acting as repulsive force (that is the qualitative meaning of the well-known Einstein relation). Gradients can be increased by increasing the attractive interaction potential relative to the diffusion. The concentrations $u_{b}$ of the black and $u_{w}$ of the white particles fulfill

[TABLE]

If the diffusive repulsion is dominating, then the solution is a perfect gray distribution. Otherwise, the black and the white particles are separated by a minimal surface.

In the following we give a brief mathematical description for this non-local phase segregation problem; a detailed discussion can be found in [10] and the references therein. The problem is given by

[TABLE]

where $\Omega\subset\mathbb{R}^{N}$, $1\leq N\leq 3$ is a bounded Lipschitz domain, $\phi$ a convex function, the kernel $\mathcal{K}$ represents non-local attracting forces, and $w$ and $v$ are the interaction and chemical potentials. The initial value $u_{0}$ satisfies $0<u_{0}<1$ and the system is completed with homogeneous Neumann boundary conditions.

The system was rigorously derived in [16] for the case of a constant $f$ and can be seen as a non-local variant of the Cahn-Hilliard equation associated with the local Ginzburg-Landau free energy

[TABLE]

The Lyapunov functional for 1 turns out to be

[TABLE]

We consider the specialized model with a constant $f$, $\phi(u)=u\ln u+(1-u)\ln(1-u)$, and $\mathcal{K}$ being the Green’s function of the elliptic boundary value problem

[TABLE]

with constant $\sigma,m>0$.

In the short time regime the model has been used for image reconstruction and denoising as well and it is not restricted to just two species [11, 12, 18].

An appropriate numerical scheme for this problem is

[TABLE]

with the notation $z:=z(t)$, $z^{+}:=z(t+\tau)$, $[V]$ being the diagonal matrix of Voronoi volumes, and $\left[\frac{1}{\Phi^{\prime\prime}_{a}}\right]$ depending on ${\bm{u}^{+}}$ [10]. This scheme is shown to be dissipative [10, Theorem 2.7]. Moreover, if $\bm{v}^{+}\neq c$, $0<u<1$, then the implicit Euler scheme preserves the property $0<u^{+}<1$ [10, Proposition 2.9]

The time step size control in the computation is based on the predictor corrector differences of the free energy and discrete dissipation rate. Jacobian matrices are derived as analytic derivatives on simplices, and Newton’s method is applied without any damping. The initial value is

[TABLE]

The dissipative scheme approximates the expected minimal surface. After $10^{20}$ in time, the initial gray distribution is transformed through a striped pattern with bubbles into a black and white square halves divided by a straight vertical line in the middle of the square. Figures 7 and 8 show a typical evolution process: thin stripes (Fig. 7, left) decay into circles with a typical diameter related to the stripe width (Fig. 7, center). The symmetry is broken because of the rounding, small perturbations in the initial data, and because $\Omega$ is large compared to $\sigma$ and not compatible with arbitrary stripe width. The longest time scale is defined by the time which is required to straighten the borderline between these halves to a straight line.

If the discretization is non-dissipative, then the asymptotic minimal surface will be replaced by a much larger surface. The free energy is maximal for the initial value. Close to the end it decays and, after some orders of magnitude in time, looks like a constant. The remaining gradient in the solution, very small almost everywhere, is sufficient to produce some free energy at some position in the mesh in the non-dissipative discretization. This causes larger gradients and larger free energy production. In the mostly dissipative volume it is resulting in a much larger separating surface because the strong dissipation along the surface is needed to compensate the production of free energy. At the end, the computed steady state can be far from that of the true minimal energy.

4.2 Active pixel DePFET detectors

Active pixel DePFET detectors are silicon sensors based on depleted $p$-channel field-effect transistors (DePFETs) [20, 26], which are p-channel field-effect transistors (FETs) integrated onto a fully depleted silicon bulk. Entering particles generate electrons in the bulk, which are then collected at the internal gate below the transistor channel. This increases the conductivity of the FET proportionally to the amount of the stored charge, providing a measure for the photon energy. The unique advantages of detectors based on this technology are the amplification and storage of the signal charge directly inside the pixels, non-destructive readouts, very low noise, and high power-efficiency, which makes them highly suitable for X-ray astronomy and particle physics [26, 34].

An example of such a detector is the Wide Field Imager (WFI) for the Advanced Telescope for High-Energy Astrophysics (Athena) [1, 2] of the European Space Agency. The WFI is a solid-state imaging spectrometer for X-ray images in the 0.1–$15\text{\,}\mathrm{keV}$ energy band with the simultaneous spectrally and time-resolved photon detection. It is developed by a collaboration of the Halbleiterlabor of the Max-Planck-Society (MPG HLL), Erlangen Centre for Astroparticle Physics, the Institute for Astronomy and Astrophysics Tübingen, the Dresden University of Technology, PNSensor GmbH, the University of Leicester, and the L’Institut de Recherche en Astrophysique et Planétologie Toulouse [2]. The numerical simulation was carried out in collaboration of the MPG HLL with the Weierstrass Institute for Applied Analysis and Stochastics.

In semiconductors, in comparison to the previous example, the particle interaction is reversed: particles of different types (holes and electrons) attract each other whereas diffusion is the repelling force for the particles of the same type. The mathematical model is given by the van Roosbroeck equations

[TABLE]

where $w$ is the electrostatic potential, $n=n_{i}e^{w-\phi_{n}}$ and $p=n_{i}e^{\phi_{p}-w}$ are the electron and hole densities, $C$ is doping, $\phi_{n}$ and $\phi_{p}$ are quasi-Fermi potentials and $R=r(\bm{x},n,p)\left(n_{i}^{2}-np\right)$ is the recombination/generation rate with $r(\bm{x},n,p)>0$.

For the numerical simulation, it is essential for the discrete system to preserve essential properties of the analytical system, such as the maximum principle, positivity of the carrier concentrations, and current conservation. The Voronoi FVM allows to carry over this properties to the discretized equations. Numerical computations are done with the well-known Scharfetter–Gummel scheme (for a rigorous treatment see, e.g., [14]). Fig. 9 shows the simulation of the depletion front propagation in a pixel of the sensor. Because of the symmetry of the domain the computation is carried out on one half a pixel.

5 Concluding remarks

Delaunay-Voronoi meshing offers a chance to approximate the geometry of a problem while preserving the essential qualitative properties of the analytic solution. It is difficult to imagine another way to reach more than just a short time approximations with discrete schemes that violate the essential stability properties of the analytic problem such as solution bounds, dissipativity for all concentrations, and any space and time step sizes.

The knowledge is not always a strictly increasing function in time and the relations between analytic properties of (nonlinear) PDE systems, finite volume methods, Delaunay meshes, and mesh generation may have been stronger in the past. During the last decades, boundary conforming Delaunay meshes seem to be an uncommon species but there is good hope that they will experience a new revival [13].