The efficient generation of unstructured control
volumes in 2D and 3D

arXiv:cs/0202038v1 [cs.CG] 26 Feb 2002

Jacek Leszczynski, Sebastian Pluta
Technical University of Czestochowa,
Institute of Mathematics and Computer Science,
ul. Dabrowskiego 73, 42-200 Czestochowa, Poland
jale@k2.pcz.czest.pl, pluta@matinf.pcz.czest.pl

Abstract. Many problems in engineering, chemistry and physics require
the representation of solutions in complex geometries. In the paper we
deal with a problem of unstructured mesh generation for the control
volume method. We propose an algorithm which bases on the spheres
generation in central points of the control volumes.

1

Introduction

An unstructured mesh generator can solve many of the problems associated
with the structured meshes. The unstructured meshes are suitable for complex
geometries, especially when we use the control volume method [10]. The control
volume method on the one hand, is based upon the simple physical principle of a
flux balance. Nevertheless, the generated meshes have to fulfill some restrictions
for such method. One of the conditions says than lines connecting the central
point in the control volume with neighbour central points have to be perpendicular to the sides of the neighbour volumes. However, the neighbourhood of each
neighbour or cell in an unstructured mesh must be defined explicitly. This is
one disadvantage of unstructured meshes because it has to reserve large storage
of the computer memory. Nevertheless, the advantages like easy for handling
adaptability in time, ability to generate meshes about arbitrary geometries. We
can divide such meshes into the cell-vertex control volumes and the cell-centered
control volumes. In the cell-vertex volumes, all elements containing the relevant
point are applied as a control volume, so the control volumes overlap. In the cellcentered control volumes the elements are subdivided and the control volumes
do not overlap.
In this paper we propose a generator for unstructured cell-centered volumes,
which is useful for the control volume method [10]. The mesh is constructed for
a number of points randomly located in a domain and the domain border.

2

Mathematical background

We turn our attention in two and three dimensional space. Let us consider a domain Ω ⊂ R3 with smooth its boundary ∂Ω. In such domain and boundary we

generate randomly a set of points Pi (a, b, c), where i = 1..N . The points establish cell-centers, around of which we determine convex polygons. The polygons
are the non-overlapping control volumes which are applied in the control volume
method [10]. We can formulate two groups of basic assumptions necessary for
the meshes construction as follows:
1. global conditions:
– such domain has to be consistent,
– edges of some polygon have the same lengths to the corresponding edges
of neighbour polygons,
2. local conditions:
– inside a polygon one may found only one random point called the cellcenter or the point,
– an edge of the polygon has to be perpendicular to the line connecting two
points.
If we fulfill the global assumptions, we may discretize the domain Ω through its
division into the triangles in R2 space or into the tetrahedrons in R3 space. Such
discretization is a standard solution, which one may find in convex geometry
[4] through Delaunay triangulation [1,8,11]. Regarding to the global and local
conditions we take into consideration a fact, that we cannot generate polygons
in R2 which overlap neighbour polygons. The fact also exist for the generated
tetrahedrons in R3 space.
In this paper we propose a novel algorithm for generation of control volumes.
We start the algorithm in two dimensions because one can see how it works. We
consider circles with unknown radiuses which centers are located in the points,
which are randomly located in the domain Ω. Let us assume, that within the
triangle created by the points the circles only in the one point crosses. The point
is a vertice of control volumes generated around the point. Following that we
need to solve the system of equations
(b
x − al )2 + (b
y − bl )2 = rl2 , (l = 1, 2, 3).

(1)

Assuming, that the radiuses
 rl of the circles are unknown, we obtain a point
Q x
b(r1 , r2 , r3 ), yb(r1 , r2 , r3 ) which coordinates are solution of the system (1)
−(b2 −b3 )·r12 +(b1 −b3 )·r22 −(b1 −b2 )·r32 +(a21 +b21 )·(b2 −b3 )
+
2·[(a1 −a2 )·(b1 −b3 )−(a1 −a3 )·(b1 −b2 )]
(a22 +b22 )·(b1 −b3 )+(a23 +b23 )·(b1 −b2 )
− 2·[(a1 −a2 )·(b1 −b3 )−(a1 −a3 )·(b1 −b2 )]
(a −a )·r12 −(a1 −a3 )·r22 +(a1 −a2 )·r32 −(a21 +b21 )·(a2 −a3 )
+
= 2 3 2·[(a
1 −a2 )·(b1 −b3 )−(a1 −a3 )·(b1 −b2 )]
(a22 +b22 )·(a1 −a3 )−(a23 +b23 )·(a1 −a2 )
+ 2·[(a1 −a2 )·(b1 −b3 )−(a1 −a3 )·(b1 −b2 )]

x
b(r1 , r2 , r3 ) =
yb(r1 , r2 , r3 )

(2)

The indexes in formulas (1) and (2) correspond to the local case.
In the global case we have random points Pi (a, b, c), i = 1..N which are
generated inside the domain Ω and its boundary ∂Ω. We use Delaunay triangulation [1,8,11] that to find neighbourhood of point Pi defined by several points
Pj(i,1) . . . Pj(i,k) . . . Pj(i,M(i)) . The index j(i, k) is a function which establishes

a point number having neighbourhood with the i-th point. The temporal index
k varies from 1 to M (i). The M (i) is a number of points corresponding to the
i-th point. Three points Pi , Pj(i,k) , Pj(i,k+1) define a triangle. When the temporal index k + 1 exceeds M (i) then we have Pj(i,k+1) = Pj(i,1) . Moreover, we
take into consideration a condition j(i, k) 6= i. We extend a definition of formula (2) for the global case putting (ai , bi ), (aj(i,k) , bj(i,k) ), (aj(i,k+1) , bj(i,k+1) )
and (ri , rj(i,k) , rj(i,k+1) ) respectively into variables with established indexes. In
our approach, we are looking for the radiuses ri . Taking into consideration of
such fact we formulate a range of radiuses variation as


max (Li,j(i,k) − rmax j(i,k) ) < ri < rmax i ,
(3)
k

where Li,j(i,k) is a distance between two neighbour points defined as
Li,j(i) =

q
(ai − aj(i) )2 + (bi − bj(i) )2 .

(4)

Maximal radius of a circle is formed as


rmax i = min hi,j(i,k) ,

(5)

k

where hi,j(i,k) is a distance defined in acute triangles as
|(b
−b
)a −(a
−a
)b |
hi,j(i,k) = √j(i,k) j(i,k+1) i 2 j(i,k) j(i,k+1) i2 +

+

(bj(i,k) −bj(i,k+1) ) +(aj(i,k) −aj(i,k+1) )
|(aj(i,k) −aj(i,k+1) )bj(i,k+1) −(bj(i,k) −bj(i,k+1) )aj(i,k+1) |

√

.

(6)

(bj(i,k) −bj(i,k+1) )2 +(aj(i,k) −aj(i,k+1) )2

In acute angles of right and obtuse triangles we have


hi,j(i,k) = min Li,j(i,k) , Li,j(i,k+1) .
k

(7)

The symbol rmax j(i,k) represents a maximal radius of the j-th point being
a neighbour of the i-th point. Using extended formula (2) and additional expressions (3)÷(7) one can construct an unstructured mesh. Fig. 1 shows an idea
of the control volume generation using the expression (2) and the limit (3). We
can observe that the local condition - an edge of the polygon has to be perpendicular to the line connecting two points - is fulfilled. The general algorithm uses
some optimisation that to find the points (2) being the vertices of the polygons
(control volumes). In this paper we neglect the construction of an objective function because it is not the subject of our consideration. We applied Rosenbrock’s
method [6] for the local solution. For the global solution we use the soft selection
method [3] which is a global optimisation method. Prediction of the unstructured
mesh under such way is not only the one way of the mesh construction. It is possible to construct another type of mesh without restriction that expression (2)
is a solution of eqn. (1). But we have to restrict the expression (3) because it
guaranties right mesh construction for the control volume method [10]. Let us
assume that a radius of the circle overlaps other radiuses of the neighbour circles.

Fig. 1. Idea of generation of control volumes.
Of course we have to restrict the condition (3). Corresponding to previous and
present mesh construction we have

Li,j(i,k) − ri + rj(i,k) ≤ 0.
(8)
Following eqn. (8) the point coordinates given by expression (2) are not the
solution of the system (1) but they become the coordinates of the point predicted
by lines connecting two points of crossing circles. From the other hand, we can
assume that a radius of the circle does not overlap other radiuses of the neighbour
circles. In such situation we have

Li,j(i,k) − ri + rj(i,k) > 0.
(9)
One can say, that conditions (8) and (9) allow us to obtain different forms of
unstructured meshes. However, in local solution is possible some combination of
expressions (8) and (9) respectively.
Similar to previous results we consider a set of random points Pi (a, b, c) in
R3 space. Analogously to the previous considerations, we use spatial Delaunay
triangulation that to find neighbourhood of point Pi which is defined by several
points Pj(i,k,l) . The index j(i, k, l) is a function which establishes a point number having neighbourhood with the i-th point. The temporal index k describes

a number of tetrahedrons having the common i-th point. It varies from 1 to
N (i). The next index l varies from 1 to 3 and it is a point number corresponding to the k-th tetrahedron. The four points Pi , Pj(i,k,1) , Pj(i,k,2) , Pj(i,k,3) define
a k-th tetrahedron. We search radiuses of spheres which centers are located in
the points. In each tetrahedron crosses the four radiuses and they create a point
which is the vertice of the control volume. Solving the following system
(b
x − al )2 + (b
y − bl )2 + (b
z − cl )2 = rl2 , (l = 1..4),

(10)

we can find the local vertices coordinates
x
b(r1 , r2 , r3 , r4 ) =

where

Wx
Wy
Wz
, yb(r1 , r2 , r3 , r4 ) =
, zb(r1 , r2 , r3 , r4 ) =
, (11)
W
W
W
△1 2(b1 − b2 ) 2(c1 − c2 )
Wx = △2 2(b1 − b3 ) 2(c1 − c3 ) ,
△3 2(b1 − b4 ) 2(c1 − c4 )

2(a1 − a2 ) 2(b1 − b2 ) 2(c1 − c2 )
W = 2(a1 − a3 ) 2(b1 − b3 ) 2(c1 − c3 ) ,
2(a1 − a4 ) 2(b1 − b4 ) 2(c1 − c4 )
2(a1 − a2 ) △1 2(c1 − c2 )
Wy = 2(a1 − a3 ) △2 2(c1 − c3 ) ,
2(a1 − a4 ) △3 2(c1 − c4 )

2(a1 − a2 ) 2(b1 − b2 ) △1
Wz = 2(a1 − a3 ) 2(b1 − b3 ) △2
2(a1 − a4 ) 2(b1 − b4 ) △3

and
△1 = r22 − r12 + a21 − a22 + b21 − b22 + c21 − c22 ,

△2 = r32 − r12 + a21 − a23 + b21 − b23 + c21 − c23 ,
△3 = r42 − r12 + a21 − a24 + b21 − b24 + c21 − c24 .

Similar to previous mesh construction, we define lower and upper limits of the
radiuses variation


max (hi,j(i,k,l) − rmax j(i,k) ) < ri < rmax i .
(12)
k,l

The symbol hi,j(i,k,l) represents height of a triangle being the wall of the tetrahedron and for the acute triangles is defined by the following formula
q
·CC12 )2 (AA32 ·CC12 −AA12 ·CC32 )2 +
hi,j(i,k,l) = (BB12 ·CC32 −BB32AA
+
2 +BB 2 +CC 2
32
32
32
q
,
(13)
(AA ·BB −AA ·BB )2
12

32

32

12

2 +CC 2
AA232 +BB32
32

where
AA12 = ai − aj(i,k,l) ,

BB12 = bi − bj(i,k,l) ,
CC12 = ci − cj(i,k,l) ,

AA32 = aj(i,k,l+1) − aj(i,k,l) ,

BB32 = bj(i,k,l+1) − bj(i,k,l) ,
CC32 = cj(i,k,l+1) − cj(i,k,l) .

In acute angles of right and obtuse spatial triangles we have


hi,j(i,k,l) = min Li,j(i,k,l) , Li,j(i,k,l+1) ,
l

where Li,j(i,k,l) is a spatial distance defined between two points as
q
2
2
2
Li,j(i,k,l) =
ai − aj(i,k,l) + bi − bj(i,k,l) + ci − cj(i,k,l) .

(14)

(15)

When the temporal index l + 1 exceeds 3 then we put Pj(i,k,l+1) = Pj(i,k,1) .
Maximal radius of a sphere is defined as


rmax i = min Hi,j(i,k) .
(16)
k

However, the function Hi,j(i,k) is height of a tetrahedron created by the vertices
Pi , Pj(i,k,1) , Pj(i,k,2) , Pj(i,k,3) and for the acute tetrahedron we have
Hi,j(i,k) =

|A · ai + B · bi + C · ci + D|
√
A2 + B 2 + C 2

(17)

where
A = m21 · n31 − m31 · n21 ,

C = d21 · m31 − d31 · m21 ,

B = d31 · n21 − d21 · n31 ,

D = −A · aj(i,k,1) − B · bj(i,k,1) − C · cj(i,k,1)

and
d21 = aj(i,k,2) − aj(i,k,1) , m21 = bj(i,k,2) − bj(i,k,1) , n21 = cj(i,k,2) − cj(i,k,1) ,
d31 = aj(i,k,3) − aj(i,k,1) , m31 = bj(i,k,3) − bj(i,k,1) , n31 = cj(i,k,3) − cj(i,k,1) .
But for the right and obtuse tetrahedrons we have


Hi,j(i,k) = min hi,j(i,k,l) .
l

(18)

The symbol rmax j(i,k) represents a maximal radius of j-th point being a neighbour of the i-th point.

3

Example of calculations

In practical mesh construction when expression (3) is fulfilled we obtain a unique
solution. Fig. 2 shows an example of mesh created under restriction that formula (2) is a solution of the system (1). Fig. 2a presents triangulation of some
points in which generated circles. Within the triangle crosses the circles only
in one point. According to the previous description we draw control volumes
presented by Fig. 2b. Nevertheless, we can construct a mesh without restriction
that formula (2) is a solution of the system (1). We can apply conditions (8)
or (9) respectively. In this case we generate a mesh much more easier than in
previous case. Expressions (8) and (9) extend our considerations for different
generation of meshes.

Fig. 2. Example of control volumes taking into consideration that formula (2)
is a solution of the system (1): a) random points after triangulation; b) forms of
control volumes.

4

Concluding remarks

In this paper we elaborate a new method of generation of unstructured meshes.
The algorithm is prepared for the control volume method in two- and three
dimensional space. For random location of points we use Delaunay triangulation
before the mesh generation. Our method bases in 3D on circles generation in
the points. Within one triangle created by three points crosses the circles in one
point. The point is a vertice of a control volume. We extend our considerations for
the triangles in which are overlapping or non-overlapping circles. The extension
allow us to generate meshes much more easier than in previous case. In opposite
to the Thyssen polygons we can generate control volumes for the right and
obtuse triangles. However, our method is suitable and easy to use in generation
of control volumes in three dimensional space.

References
1. Chiang Y. J., Tamassia R.: Dynamic algorithms in computational geometry, Proc.
IEEE 80(9), (1992) 1412-1432
2. Elmahi I, Benkhaldoun F.: Finite volume simulation of interval of a droplet flame
ignition on unstructured meshes, Journal of Computational an Applied Mathematics, 103, (1999) 187-203
3. Galar R.: Soft selection in Random Global Adaptation in Rn . A biocybernetic
model of development, Mongraph 17, Wroclaw University of Technology, Wroclaw
(Poland), (1990) (in polish)
4. Gruber P.M., Wills J.M.: Handbook on Convex Geometry, vol. A, B, NorthHolland, (1993)
5. Jorgenson P. C. E., Pletcher R. H.: An implicit numerical scheme for the simulation
of interval viscous flow on unstructured grids, Computers & Fluids 5, (1996) 447466
6. Kunzi H. P., Tzschach H. G., Zehnder C. A.: Numerical methods of mathematical
optimization with Algol and Fortran programmes, Academic Press (1971)
7. Neises J., Steinbach I.: Finite element integration for the control volume method,
Communications in numerical methods in engineering 12, (1996) 543-555
8. Orkisz J.: Finite difference method, in: Kleiber M. (ed.), Computational methods
in solid mechanics, PWN, Warsaw (1995) (in polish)
9. O’Rourke J.: Computational geometry in C, Cambridge University Press (1995)
10. Patankar V. S.: Numerical heat transfer and fluid flow, Hemisphere, New York,
(1980)
11. Renka R. J., Algorithm 772: STRIPACK: Delaunay triangulation and Voronoi
diagram on the surface of a sphere, ACM Trans. Math. Soft. 23(3), (1997) 416434

