6.2. Finite Volume Method on Triangular Grids 263

.a2

. . (a1 - a4 are cocyclic)

a3 a5

.

. a1

a4

Figure 6.6. Degenerate and regular Voronoi vertex.

Now the elements „¦i (control volumes) of the partition of „¦ required for

the de¬nition of the ¬nite volume method can be introduced as follows:

„¦i := „¦i © „¦ , i ∈ Λ.

˜

As a consequence, the domains „¦i need not necessarily be convex if „¦ is

nonconvex (cf. Figure 6.5).

Furthermore, the following notation will be used:

j ∈ Λ \ {i} : ‚„¦i © ‚„¦j = … , i ∈ Λ,

Λi :=

for the set of indices of neighbouring nodes,

‚„¦i © ‚„¦j , j ∈ Λi , for a joint piece of the

“ij :=

boundaries of neighbouring control volumes,

mij for the length of “ij .

The dual graph of the Voronoi diagram is de¬ned as follows:

Any pair of points ai , aj such that mij > 0 is connected by a straight-line

segment. In this way, a further partition of „¦ with an interesting property

results.

Theorem 6.4 If all Voronoi vertices are regular, then the dual graph coin-

cides with the set of edges of a triangulation of the convex hull of the given

point set.

This triangulation is called a Delaunay triangulation.

If among the Voronoi vertices there are degenerate ones, then a tri-

angulation can be obtained from the dual graph by a subsequent local

triangulation of the remaining m-polygons (m ≥ 4). A Delaunay triangula-

tion has the interesting property that two interior angles subtended by any

given edge sum to no more than π. In this respect Delaunay triangulations

satisfy the ¬rst part of the angle condition formulated in Section 3.9 for

the maximum principle in ¬nite element methods.

Therefore, if „¦ is convex, then we automatically get a triangulation to-

gether with the Voronoi diagram. In the case of a nonconvex domain „¦,

certain modi¬cations could be required to achieve a correct triangulation.

.

264 6. Finite Volume Method

. .

.. This edge has to be

. removed from the

.. Delaunay triangulation.

Figure 6.7. Delaunay triangulation to the Voronoi diagram from Figure 6.5.

The implication

’

Voronoi diagram Delaunay triangulation ,

which we have just discussed, suggests that we ask about the converse

statement. We do not want to answer it completely at this point, but we

give the following su¬cient condition.

Theorem 6.5 If a conforming triangulation of „¦ (in the sense of ¬nite

element methods) consists of nonobtuse triangles exclusively, then it is

a Delaunay triangulation, and the corresponding Voronoi diagram can be

constructed by means of the perpendicular bisectors of the triangles™ edges.

We mention that the centre of the circumcircle of a nonobtuse triangle is

located within the closure of that triangle.

In the analysis of the ¬nite volume method, the following relation is

important.

Lemma 6.6 Given a nonobtuse triangle K with vertices aik , k ∈ {1, 2, 3},

then for the corresponding parts „¦ik ,K := „¦ik © K of the control volumes

„¦ik , we have

1 1

|K| ¤ |„¦ik ,K | ¤ |K| , k ∈ {1, 2, 3} .

4 2

The Donald diagram

In contrast to the Voronoi diagram, where the construction starts from a

given point set, the starting point here is a triangulation Th of „¦, which is

allowed to contain obtuse triangles.

Again, let K be a triangle with vertices aik , k ∈ {1, 2, 3}. We de¬ne

„¦ik ,K := x ∈ K »j (x) < »k (x), j = k ,

where »k denote the barycentric coordinates with respect to aik (cf. (3.51)).

Obviously, the barycentre satis¬es aS = 1 (ai1 + ai2 + ai3 ), and (see, for

3

comparison, Lemma 6.6)

3 |„¦ik ,K | = |K| , k ∈ {1, 2, 3} . (6.4)

This relation is a simple consequence of the geometric interpretation of

the barycentric coordinates as area coordinates given in Section 3.3. The

6.2. Finite Volume Method on Triangular Grids 265

.

ai2

„¦i2 ,K

.

. aS

„¦i1 ,K ai1

„¦ i3 ,K

.

ai3

Figure 6.8. The subdomains „¦ik ,K .

required control volumes are de¬ned as follows (see Figure 6.8):

i ∈ Λ.

„¦i := int „¦i,K ,

K:‚K ai

The family {„¦i }i∈Λ is called a Donald diagram.

The quantities “ij , mij , and Λi are de¬ned similarly as in the case of

the Voronoi diagram. We mention that the boundary pieces “ij are not

necessarily straight, but polygonal in general.

6.2.2 Finite Volume Discretization

The model under consideration is a special case of equation (6.1). Instead

of the matrix-valued di¬usion coe¬cient K we will take a scalar coe¬cient

k : „¦ ’ R, that is, K = kI. Moreover, homogeneous Dirichlet boundary

conditions are to be satis¬ed. So the boundary value problem reads as

follows:

’∇ · (k ∇u ’ c u) + r u =f in „¦ ,

(6.5)

u =0 on ‚„¦ ,

with k, r, f : „¦ ’ R, c : „¦ ’ R2 .

The Case of the Voronoi Diagram

Let the domain „¦ be partioned by a Voronoi diagram and the correspond-

ing Delaunay triangulation. Due to the homogeneous Dirichlet boundary

conditions, it is su¬cient to consider only those control volumes „¦i that

are associated with inner nodes ai ∈ „¦. Therefore, we denote the set of

indices of all inner nodes by

Λ := i ∈ Λ ai ∈ „¦ .

In the ¬rst step, the di¬erential equation (6.5) is integrated over the single

control volumes „¦i :

’ ∇ · (k ∇u ’ c u) dx + i ∈ Λ.

r u dx = f dx , (6.6)

„¦i „¦i „¦i

266 6. Finite Volume Method

The application of Gauss™s divergence theorem to the ¬rst integral of the

left-hand side of (6.6) yields

∇ · (k ∇u ’ c u) dx = ν · (k ∇u ’ c u) dσ .

„¦i ‚„¦i

Due to ‚„¦i = ∪j∈Λi “ij (cf. Figure 6.9), it follows that

∇ · (k ∇u ’ c u) dx = νij · (k ∇u ’ c u) dσ ,

„¦i “ij

j∈Λi

where νij is the (constant) outer unit normal to “ij (with respect to „¦i ).

In the next step we approximate the line integrals over “ij .

.

.

. ai

“ij

aj

.

Figure 6.9. The edge “ij .

First, the coe¬cients k and νij · c are approximated on “ij by constants

µij > 0, respectively γij :

k|“ij ≈ µij = const > 0 , νij · c|“ij ≈ γij = const .

In the simplest case, the approximation can be realized by the correspond-

ing value at the midpoint a“ij of the straight-line segment “ij . A better

choice is

±

1 νij · c dσ , mij > 0 ,

mij “ij

γij := (6.7)

νij · c(a“ij ) , mij = 0 .

We thus obtain

∇ · (k ∇u ’ c u) dx ≈ [µij (νij · ∇u) ’ γij u] dσ .

„¦i “ij

j∈Λi

The normal derivatives are approximated by di¬erence quotients; that is,

u(aj ) ’ u(ai )

νij · ∇u ≈ with dij := |ai ’ aj | .

dij

This formula is exact for such functions that are linear along the straight-

line segment between the points ai , aj . So it remains to approximate the

integral of u over “ij . For this, a convex combination of the values of u at

6.2. Finite Volume Method on Triangular Grids 267

the nodes ai and aj is taken:

u|“ij ≈ rij u(ai ) + (1 ’ rij ) u(aj ) ,

where rij ∈ [0, 1] is a parameter to be de¬ned subsequently. In general, rij

depends on µij , γij , and dij .

Collecting all the above approximations, we arrive at the following

relation:

∇ · (k ∇u ’ c u) dx

„¦i

u(aj ) ’ u(ai )

≈ ’ γij [rij u(ai ) + (1 ’ rij ) u(aj )]

µij mij .

dij

j∈Λi

To approximate the remaining integrals from (6.6), the following formulas

are used:

r u dx ≈ with mi := |„¦i | ,

r(ai ) u(ai ) mi =: ri u(ai ) mi ,

„¦i

≈

f dx f (ai ) mi =: fi mi .

„¦i

Instead of ri := r(ai ) or fi := f (ai ), the approximations

1 1

ri := r dx respectively fi := f dx (6.8)

mi mi

„¦i „¦i

can also be used. Denoting the unknown approximate values for u(ai ) by

ui , we obtain the following linear system of equations:

ui ’ uj

+ γij [rij ui + (1 ’ rij ) uj ] mij + ri ui mi

µij

dij (6.9)

j∈Λi

= f i mi , i ∈ Λ .

This representation clearly indicates the a¬nity of the ¬nite volume method

to the ¬nite di¬erence method. However, for the subsequent analysis it is

more convenient to rewrite this system of equations in terms of a discrete

variational equality.

Multiplying the ith equation in (6.9) by arbitrary numbers vi ∈ R and

summing the results up over i ∈ Λ, we get

±

ui ’ uj

+ γij [rij ui + (1 ’ rij ) uj ] mij + ri ui mi

vi µij

dij

i∈Λ j∈Λi

= fi vi mi .

i∈Λ

Further, let Vh denote the space of continuous functions that are piecewise

linear over the (Delaunay) triangulation of „¦ and that vanish on ‚„¦. Then

the values ui and vi can be interpolated in Vh ; that is, there are unique

268 6. Finite Volume Method

uh , vh ∈ Vh such that uh (ai ) = ui , vh (ai ) = vi for all i ∈ Λ. The following

discrete bilinear forms on Vh — Vh can then be de¬ned:

mij

µij (ui ’ uj )

a0 (uh , vh ) := vi ,

h

dij

i∈Λ j∈Λi

[rij ui + (1 ’ rij ) uj ] γij mij ,

bh (uh , vh ) := vi

i∈Λ j∈Λi

dh (uh , vh ) := ri ui vi mi ,

i∈Λ

a0 (uh , vh )

ah (uh , vh ) := + bh (uh , vh ) + dh (uh , vh ) .

h

Finally, for two continuous functions v, w ∈ C(„¦), we set

w, v := wi vi mi ,

0,h

i∈Λ

where vi := v(ai ), wi := w(ai ).

Remark 6.7 ·, · 0,h is a scalar product on Vh . In particular, the following

norm can be introduced:

vh ∈ Vh .

vh := vh , vh 0,h , (6.10)

0,h

In (3.136) a discrete (L2 -) norm for a general ¬nite element space vh has

been de¬ned using the same notation. This multiple use seems to be accept-

able, since for regular triangulations both norms are equivalent uniformly

in h (see Remark 6.16 below).

Now the discrete variational formulation of the ¬nite volume method is

this:

Find uh ∈ Vh such that

for all vh ∈ Vh .

ah (uh , vh ) = f, vh (6.11)

0,h

Up to now, the choice of the weighting parameters rij has remained open.

For this, two cases can be roughly distinguished:

(1) There exists a pair of indices (i, j) ∈ Λ — Λ such that µij |γij |dij .

|γij |dij .

(2) There is no such pair (i, j) with µij

In the second case, an appropriate choice is rij ≡ 1 . To some extent, this

2

can be seen as a generalization of the central di¬erence method to nonuni-

form grids. The ¬rst case corresponds to a locally convection-dominated

situation and requires a careful selection of the weighting parameters rij .

This will be explained in more detail in Section 9.3.

In general, the weighting parameters are of the following structure:

γij dij

rij = R , (6.12)

µij

6.2. Finite Volume Method on Triangular Grids 269

γ d

where R : R ’ [0, 1] is some function to be speci¬ed. The argument ij ijij

µ

is called the local P´clet number. Typical examples for this function R are

e

1

R(z) = [sign (z) + 1] , full upwinding,

2

(1 ’ „ )/2 , z < 0 , 2

„ (z) := max 0, 1 ’

R(z) = ,

(1 + „ )/2 , z ≥ 0 , |z|

1 z

R(z) = 1 ’ 1’ z , exponential upwinding .

e ’1

z

All these functions possess many common properties. For example, for all

z ∈ R,

[1 ’ R(z) ’ R(’z)] z

(P1) = 0,

R(z) ’ 1 z ≥

(P2) 0, (6.13)

2

1 ’ [1 ’ R(z)] z ≥

(P3) 0.

1

Note that the constant function R = satis¬es the conditions (P1) and

2

(P2) but not (P3).

The Case of the Donald Diagram

Let the domain „¦ be triangulated as in the ¬nite element method. Then,

following the explanations given in the second part of Section 6.2.1, the

corresponding Donald diagram can be created.

The discrete bilinear form in this case is de¬ned by

ah (uh , vh ) := k ∇uh , ∇vh + bh (uh , vh ) + dh (uh , vh ) ;

0

that is, the principal part of the di¬erential expression is discretized as in

the ¬nite element method, where bh , dh , and Vh are de¬ned as in the ¬rst

part of this section.

6.2.3 Comparison with the Finite Element Method

As we have already seen in Example 6.1, it may happen that a ¬nite volume

discretization coincides with a ¬nite di¬erence or ¬nite element discretiza-

tion. We also mention that the control volumes from that example are

exactly the Voronoi polygons to the grid points (i.e., to the nodes of the

triangulation).

Here we will consider this observation in more detail. By {•i }i∈Λ we de-

note the nodal basis of the space Vh of continuous, piecewise linear functions

on a conforming triangulation of the domain „¦.

Lemma 6.8 Let Th be a conforming triangulation of „¦ (in the sense of

¬nite element methods), all triangles of which are nonobtuse, and consider

the corresponding Voronoi diagram in accordance with Theorem 6.5. Then,

for an arbitrary triangle K ∈ Th with vertices ai , aj (i = j), the following

270 6. Finite Volume Method

relation holds:

mKij

∇•j · ∇•i dx = ’ ,

dij

K

where mK is the length of the segment of “ij that intersects K.

ij

Proof: Here we use some of the notation and the facts prepared at the

beginning of Section 3.9. In particular, ±K denotes the interior angle of K

ij

that is located in opposite the edge with vertices ai , aj . Next, the following

equality is an obvious fact from elementary geometry: 2 sin ±K mK = ij ij

K

cos ±ij dij . It remains to recall the relation

1

∇•j · ∇•i dx = ’ cot ±K

ij

2

K

2

from Lemma 3.47, and the statement immediately follows.