Applied Computational Fluid Dynamics Techniques - Wiley Episode 1 Part 4 pot

GRID GENERATION 61

D1. Assume given a boundary point distribution.

D2. Generate a Delaunay triangulation of the boundary points.

D3. Using the information stored on the background grid and the

sources, compute the desired element size and shape for the points

of the current mesh.

D5. nnewp=0 ! Initialize new point counter

D6. do ielem=1,nelem ! Loop over the elements

Define a new point inewp at the centroid of ielem;

Compute the distances dispc(1:4) from inewp to the four nodes

of ielem;

Compare dispc(1:4) to the desired element size and shape;

If any of the dispc(1:4) is smaller than a fraction

α of the desired element length: skip the element (goto D6);

Compute the distances dispn(1:nneip) from inewp to the new

points in the neighbourhood;

If any of the dispn(1:nneip) is smaller than a fraction

β of the desired element length: skip the element (goto D6);

nnewp=nnewp+1 ! Update new point list

Store the desired element size and shape for the new point;

enddo

D7. if(nnewp.gt.0) then

Perform a Delaunay triangulation for the new points;

goto D.5

endif

The procedure outlined above introduces new points in the elements. One can also

introduce them at edges (George and Borouchaki (1998)). In the following, individual aspects

of the general algorithm outlined above are treated in more detail.

3.7.1. CIRCUMSPHERE CALCULATIONS

The most important ingredient of the Delaunay generator is a reliable and fast algorithm for

checking whether the circumsphere (-circle) of a tetrahedron (triangle) contains the points to

be inserted. A point xp lies within the radius Ri of the sphere centred at xc if

d2

p = (xp − xc) · (xp − xc)<R2

i . (3.31)

This check can be performed without any problems unless |dp − Ri|

2 is of the order of

the round-off of the computer. In such a case, an error may occur, leading to an incorrect

rejection or acceptance of a point. Once an error of this kind has occurred, it is very difficult

to correct, and the triangulation process breaks down. Baker (1987) has determined the

following condition:

Given the set of points P := x1, x2,..., xn with characteristic lengths dmax = max|xi −

xj | ∀i

= j and dmin = min|xi − xj | ∀i

= j , the floating point arithmetic precision required

for the Delaunay test should be better than

 =

dmin

dmax 2

. (3.32)

62 APPLIED COMPUTATIONAL FLUID DYNAMICS TECHNIQUES

Consider the generation of a mesh suitable for inviscid flow simulations for a typical transonic

airliner (e.g. Boeing-747). Taking the wing chord length as a reference length, the smallest

elements will have a side length of the order of 10−3L, while far-field elements may be

located as far as 102L from each other. This implies that  = 10−10, which is beyond the

10−8 accuracy of 32-bit arithmetic. For these reasons, unstructured grid generators generally

operate with 64-bit arithmetic precision. When introducing points, a check is conducted for

the condition

|dp − Ri|

2 < , (3.33)

where  is a preset tolerance that depends on the floating point accuracy of the machine. If

condition (3.33) is met, the point is rejected and stored for later use. This ‘skip and retry’

technique is similar to the ‘sweep and retry’ procedure already described for the AFT. In

practice, most grid generators work with double precision and the condition (3.33) is seldom

met.

A related problem of degeneracy that may arise is linked to the creation of very flat

elements or ‘slivers’ (Cavendish et al. (1985)). The calculation of the circumsphere for a

tetrahedron is given by the conditions

(xi − xc) · (xi − xc) = R2, i = 1, 4, (3.34)

yielding four equations for the four unknowns xc, R. If the four points of the tetrahedron lie

on a plane, the solution is impossible (R → ∞). In such a case, the point to be inserted is

rejected and stored for later use (skip and retry).

3.7.2. DATA STRUCTURES TO MINIMIZE SEARCH OVERHEADS

The operations that could potentially reduce the efficiency of the algorithm to O(N1.5) or

even O(N2) are:

(a) finding all tetrahedra whose circumspheres contain a point (step B3);

(b) finding all the external faces of the void that results due to the deletion of a set of

tetrahedra (step B5);

(c) finding the closest new points to a point (step D6);

(d) finding for any given location the values of generation parameters from the background

grid and the sources (step D3).

The verb ‘find’ appears in all of these operations. The main task is to design the best

data structures for performing the search operations (a)–(d) as efficiently as possible. As

before, many variations are possible here, and some of these data structures have already

been discussed for the AFT. The principal data structure required to minimize search

overheads is the ‘element adjacent to element’ or ‘element surrounding element’ structure

esuel(1:nfael,1:nelem) that stores the neighbour elements of each element. This

structure, which was already discussed in Chapter 2, is used to march quickly through the

grid when trying to find the tetrahedra whose circumspheres contain a point. Once a set

of elements has been marked for removal, the outer faces of this void can be obtained by

interrogating esuel. As the new points to be introduced are linked to the elements of the

current mesh (step D6), esuel can also be used to find the closest new points to a point.

Furthermore, the equivalent esuel structure for the background grid can be used for fast

interpolation of the desired element size and shape.