Anisotropic 3D Grid Generation with Intersection-Based Geometry Interface
Ilja Schmelzer IAAS, Mohrenstr. 39 D-10117 Berlin November 29, 1993
In this paper we present a new interface for geometry description. This interface is based on four intersection routines. It allows to use geometry description without explicit boundary description, for example using characteristic functions and boolean operations. This seams to be useful especially for time dependent geometries. We show how to create and change this geometry description, and how to use this interface for grid generation. We discuss grid quality criteria and possibilities to create and manage anisotropic grids. At last we describe an algorithm for 3D grid generation based on this interface which was implemented by the author. This algorithm combines an octree-like renement process, shift of nodes to the boundary and Delaunay grid generation.
Abstract
1 Introduction
Usually a grid generator requires the full information about the boundary grid, especially a list of all boundary vertices, boundary lines and boundary
This research was supported in part by the Institute for Mathemetics and it's Applications with funds provided by the National Science Foundation.
1
faces and a description of each of them. Such a boundary description can be created, for example, by standard CAD systems. But there are situations where it is easier or more natural to dene a geometry in another way, for example using a characteristic function which is greater zero inside and later zero outside the region. Then we have no list of boundary parts. We dene here an alternative interface for the geometry description which makes it possible to use such incomplete geometry descriptions. We describe how the interface can be created and manipulated for dierent types of input data. Our geometry interface is reduced compared with the standard CAD geometry descriptions. For example, we have no lists of regions, boundary faces, edges and vertices. Therefore topological errors cannot be excluded in general. We discuss possibilities to avoid such errors. Another problem of grid generation is anisotropic renement. Most of grid generators for unstructured grids allows local renement, but create only isotropically rened grids. But especially for 3D problems anisotropic renement can essentially reduce the number of nodes in the grid. We discuss some problems connected with anisotropic grids, especially grid quality measures. Unfortunately, many usual techniques for isotropical grids cannot be used for anisotropic grid generation. At last we describe an algorithm for 3D anisotropic grid generation based on our alternative geometry interface. The algorithm starts with an octreelike renement process. Then nodes near the boundary will be shifted to the boundary. Then we create the Delaunay grid for the resulting point set. At last, some corrections near the boundary are necessary. The algorithm is implemented in ANSI C. The speed is approximately linear, with 3-5 ms per node CPU time on a VAX 4000/90 workstation.
2 problems of usual geometry description concepts
For a specic, xed geometrical conguration a geometry description can be generated by usual CAD systems. But such CAD systems cannot handle a geometry changing in time. This may be caused by mechanical deformation or a chemical reaction on a surface. Then the geometry description has to be changed automatically at every time step. The available input data may 2
be very dierent. There may be a (scalar or direction dependent) reaction rate on a boundary grid, a concentration in a region, a velocity eld and so on. This usually has to be done by the application programmer, because only he knows the velocities, reaction rates and so on which are necessary to compute the new boundary position. At rst let's consider the process of changing the "standard" geometry description. This geometry description contains: a list of all geometrical objects, that means regions, boundary faces, boundary lines and boundary vertices. a set of coincidence and neighbourhood relations between the geometrical objects. a description of every boundary object. This will be usually a discretization of this object (boundary grid for a face, polygon for an edge, coordinates of a vertex). In CAD systems also other possibilities (spline approximations or analytical denitions) are usual, but we don't consider them here. Assume we have computed some velocity eld on the boundary grid. Then we shift the boundary grid. It is easy to see that we obtain a lot of problems: caustics: Also if the initial grid and the velocity eld are smooth, the shifted grid often will not be smooth. There may be caustics. A typical caustic in 2D consists of two "turning points" and one self-intersection. But there may be also other types of caustics. In 3D caustics are very complex global objects. Instead of the "turning points" in 2D we have now lines with dicult global behaviour. It is necessary to detect caustics and to eliminate them. But already in 2D this is a dicult task. Often only caustics with self-intersections will be detected, but already in the case of non-constant isotropic reaction rate there may be caustics without self-intersections. intersections: There may be also other intersections of the shifted boundary grid, may be with other boundary grids, but also self-intersections. These intersections also have to be detected, but they cannot be eliminated, but require a change of the topology. This requires a very accurate classication of the intersection, because every failure leads to an 3
obsolete data structure or global errors. Especially we have to consider rounding errors. In 2D we have to nd an intersection point, in 3D we have a global intersection line which may intersect other boundary lines or caustics. The intersecting parts are obviously in the same space region, but usually "far away" in the data structure of the boundary grid. To make the test fast enough we have to use search tree structures. We see, that changing the standard geometry description is a very dicult. So, the usual method is to make small time steps to avoid caustics and not to consider topological changes. So it seams natural to look for alternatives, which allow big time steps and topological changes. An example of such an alternative is a "characteristic function" of a region which is positive in the region and negative outside. Often it is easier to dene such a characteristic function for the new region (for example a minimal distance to a boundary) as to compute the new boundary grid. Sometimes such a characteristic function may be natural for the given physical problem (some critical concentration). The natural generalization for the case of more than two dierent regions is the "region-function" | a function dening the region a point with given coordinates lies in. In some sense this seams to be the minimal geometry description | if this function cannot be dened, the geometry is not correctly dened. On the other hand, using this minimal interface we also obtain a lot of problems:
the problem of nding the correct boundary now has the grid generator. Obviously, using this geometry description you also cannot avoid topological errors, if your grid is not ne enough to detect very small regions and so on. For time dependent geometry it is necessary to nd the position of the boundary with small error. Otherwise, especially for small time steps, the error in the velocity of the boundary becomes very big. If we have only the region-function, we need an iteration process to compute the boundary position | not an ideal solution. It is not possible to transfer in this way boundary data | for example dierent boundary types, boundary charge distributions and so on. 4
special boundary lines cannot be detected | edges and vertices usually will be "rounded" by the grid generator.
3 denition of the interface
The interface we consider here is based on the second concept | a function dening the region for given coordinates. But we have modied the concept so that it is possible to transfer more information through the interface. In some sense this interface can be considered as the maximal interface that can be approximated by the minimal interface. So, if you have a geometry dened by the minimal interface, you can use it as the input of the grid generator. But if you have more detailed information about the geometry, you have the possibility to transfer some of this information through the interface or to use faster and more accurate algorithms for the interface. Since in this concept it is possible to use the minimal interface as input, the main problem of this interface | the possibility of topological errors | cannot be completely solved. But we obtain new possibilities to avoid such errors. So the starting point of the interface we want to dene is the regionfunction: returning a descriptor of the region the point with the coordinates x lies in. This descriptor in principle may be a pointer to some data, but we prefer to use an integer value enumerating the regions.
int Region(float x[DIM]);
3.1 the node data type
The rst question we have to consider is the allocation and computation of other nodal data (function values) which are necessary in the process of grid generation. For example, assume, that some function values are necessary for renement criteria. Than it is necessary to compute and store them in the process of grid generation. The storage allocation has to be done by the grid generator ("knowing" the current node number), but the data have to be computed by the geometry description ("knowing" the data on the old grid). 5
So, node functions can (and have to) be interpolated in the process of grid generation, and the natural way is to use the region-function to interpolate these data. Then they have to be transferred to the renement criteria. The grid generator only has to allocate storage for the data and to transfer the data. We have realized this using an abstract node data type containing all information about the node:
typedef structf float x[DIM], region dom, ... g Node;
So the region-function now has the form
int Region(Node *node);
3.2 data for the grid search algorithm
Another change in the interface was done to allow a fast implementation of the region-function for a geometry dened by a grid (for example the grid of the previous time step). Let us consider now algorithms nding a point with given coordinates in a grid. There are two interesting algorithms for this search:
The quadtree/octree method. A special search data structure (called quadtree in 2D, octree in 3D) has to be created to search elements near the given point. For every node we have operation for a grid with nodes. The neighbourhood search. Beginning with some start element we test if the point is inside the element. If this test fails because the nodes lies behind a side of the element we go to the neighbour element at this side.
O(log n) n
The eciency of the neighbourhood search depends of the start element we use. It seams the best to start with the element containing the nearest point we have found before. Then the eciency of the algorithm depends on the point order. As an example consider a regular rened grid with the same density for the old and the new grid. Then the nodes of the highest renement level usually have a direct neighbour of a lower level which was created before. The "distance of one neighbour" needs a xed time (This will be true also for local rened grids if the density of the old and the 6
new grid are nearly the same). For lower level nodes we have a two times bigger distance to the nearest node created before. For a regular, isotropically rened rectangluar grid we can explicitely compute this time. For an 1D grid we have nest level nodes with distance 1, nodes with distance 2, nodes with distance 4 and so on. The result is an log time behaviour for the full grid. For a 2D or 3D grid we can have a better behaviour if we use a good renement order. Then we will have (for 2D) 3n=4 nodes of the nest level with distance 1, 3n=16 nodes of the rst coarser level with distance 2 and so on with resulting linear behaviour. That means, we have a chance to get a linear behaviour using the neighbourhood search algorithm with a good start element. But to realize this we have to change the interface denition to organize the transfer of the information about the nearest previously created node (which is available only for the grid generator) and the element it was found in (which is an object of the old grid and hidden from the grid generator by the region-function). To do this we include the pointer to the data of the nearest previously found point into the region-function call:
n=2 n=4 n=8 n n
int Region(Node *nnew, Node *nold);
Into the Node data structure we include the descriptor of the element where we have found the node (usually an integer number). The regionfunction now has to write the result of the search into the Node structure and can read it from nearest. For the rst node nold will be NULL.
3.3 the position of boundary faces
Consider now the problem of exact denition of the boundary position. For time dependent problems we need the boundary position with high accuracy. That's why it is not enough to dene only the region for regular grid nodes and to consider grid nodes near the boundary as the boundary nodes. If we have found an edge with dierent region numbers at the ends we have to compute the intersection of this edge with a boundary. In principle it is possible to do this using the bisection algorithm and the region-function. Depending of the accuracy we need we have a xed number of region-function calls. The number of edges intersecting boundaries is small compared with the number of nodes especially if the grid is ne, so this is time consuming compared with other parts of the algorithm only for 7
coarse grids. But usually there are also better algorithms possible depending of the realization of the geometry description. Especially we can get more information about the intersection | special boundary types, function values dened on the boundary and so on. Therefore it is useful to include a special function into the interface dening the intersection with boundary faces: Input are the node data for the ends of the edge 1 and 2. The output data are parts of the Node structure of the "intersection node" int, especially the coordinates and the boundary face identier. int is dened as the rst intersection with a boundary with the line from 1 to 2. We have a default implementation using the bisection algorithm with iterative calls of the region-function. So it is not necessary to implement this function. But this implementation is not ideal. For example, we may not nd the rst intersection point if there are many, we obtain only an approximation and cannot transfer boundary data. So for some interesting input data we have special implementations: For the search in a grid there is an analogon of the neighbourhood search | the neighbourhood search along the edge. We nd here the exact position, really the rst intersection with the boundary face, the algorithm will be faster, and we have a possibility to interpolate nontrivial boundary data from the grid. For a change of the geometry with a characteristic function we use linear interpolation to nd the boundary face intersection | the zero value of the characteristic function. This algorithm will be usually faster and allows higher accuracy, especially if the function is approximately linear. Therefore we recommend to use characteristic functions with nearly linear behaviour near the border.
n n n n n n
int Face(Node *nint, Node *n1, Node *n2);
3.4 the position of boundary lines
The same problem as for the boundary face position we have also for the position of boundary lines. For many applications a good approximation for the position of edges is necessary. It is possible to compute the position of an boundary line by an iteration process of calls of the region- and face-function, but for special geometry data there are algorithms which will be faster, more 8
accurate and allow the transfer of other data. Therefore it seams useful to include another function dening the position of boundary lines into the interface:
int Line( Node *nint, Node *nface, Node *n1, Node *n2, Node *n3);
The input data are the three nodes of a triangle 1 3 and an intersection f ace of the line from 1 to 2 with a boundary face. This intersection is assumed to be the result of a previous call of the face- or line-function. There are two possible results:
n ; :::; n n n n
There may be other intersections of this boundary face with the border of the triangle. Than int is the rst of these intersections if we go along the intersection curve of the boundary face with the triangle. There may be no other intersections of this boundary face with the border of the triangle. That means, the intersection curve of the border with the triangle ends inside the triangle. Than int is this end | the intersection of a boundary line with the triangle.
n n
We also have implemented a default algorithm (an analogon of the bisection algorithm for the boundary faces) using calls of the region- and the face-function (So using a fast face-function makes also this algorithm faster). In the following we describe this algorithm: At rst we "go around" the triangle from the given intersection f ace to the next intersections with a boundary using face-calls. If we nd an intersection with the same face number we return this intersection. This is the analogon of "comparing the sign" of the bisection algorithm. Here we may return an incorrect result if there are dierent intersections. If the triangle is smaller than the necessary accuracy we return the middle point of the triangle. Else we divide the triangle into four parts and call recursively the line-function for the small triangle containing the given intersection f ace. If the output is an intersection of the border with the outside or an inner intersection with an boundary line we return this intersection. Else we go to the neighbour and call the line-function for this neighbour with the output of the previous call as input. It is possible to get an innite cycle (since the line-function can make errors). So we have to break down this loop after a given number of steps. Than we return the last (inner) border intersection
n n
9
we have found as an inner intersection with an edge. That means we also can get erroneous boundary lines dividing the border into dierent parts using this algorithm. Because of the possible errors of this algorithm it seams useful to have special implementations, especially for grids, there it is possible to nd the exact solution.
3.5 the position of boundary vertices
The same has to be done to dene the position of boundary vertices. We have included the following call into the interface:
int Vertex (Node *nint, Node *nline, Node *n1, Node *n2, Node *n3, Node *n4);
Input are the four nodes of a tetrahedron and an intersection line of the rst side (n1 2 3 4) with an edge. Output int are the data of the next intersection of this boundary line with the border of the tetrahedron or the vertex | the end of the boundary line | inside the tetrahedron. An equivalent default implementation calling region-, face- and line-function is available. A special (fast and exact) realization for a grid is also possible (but at the current state not implemented).
n ;n ;n ;n n
3.6 degenerate cases
For all of the previous function we have dened their behaviour only for the case of "general position". But what we have to do in degenerate cases? There aren't any special output conventions for degenerate cases. If the answer is not unique, you have to give one of the possible answers. For example if the input point of the region-function lies on a boundary, every region containing this boundary can be used. But you have to be sure that you obtain in this way a consistent picture. So if you have returned dierent region numbers for two points, you must be able to return a boundary face intersection for this edge. Considering such degenerate cases and problems caused by rounding errors is the most dicult part of the implementation of such an interface.
10
3.7 coincidence functions
Using this algorithm we cannot avoid topological errors, if the grid density dened by the renement criteria is to small to detect all regions. But it is possible to do very much to avoid such errors testing the coincidence of dierent geometrical objects. Assume, for example, we have a very thin layer which intersects a grid edge. Using the face-function we have found one of the two intersections of the face of the layer with the edge. How can we nd the other? If we have some test of the coincidence of the boundary face number and the region numbers of the ends of the grid edge we can easily detect that there is a problem and use further renement to solve it. There will be a lot of other possible errors which can be detected using such coincidence information. So it seams useful to include such tests into the interface:
int TestRF(int region, int face); int TestRL(int region, int line); int TestRV(int region, int vertex); int TestFF(int face, int face); int TestFL(int face, int line); int TestFV(int face, int vertex); int TestLL(int line, int line); int TestLV(int line, int vertex);
with the possible answers 0 - no, 1 - yes, 2 - unknown. In the case of the minimal interface we have a priori no boundary face numbers. But for every boundary face point we have found by the bisection algorithm we know the two (dierent) region numbers of the nearest points we have considered. We can use these region numbers to create a default boundary face number. Using this number we can get a nontrivial implementation of the TestRF- and TestFF-functions. Remark that for a boundary face there must not be dened a unique left and right region. 11
int TestVV(int vertex, int vertex);
3.8 making recursive calls possible
We have considered possibilities to describe a geometry, for example given by a grid or by a characteristic function. But for time-dependent problems we need a possibility to change a given geometry. This also can be done using a characteristic function for the part which has to be changed. Assume that we have dened such a function. To compute the result of the regionfunction we have to compute the result of the function and the result of the old region-function. So we can have recursive calls of the region-function. There may be recursive calls also for the other functions. It seams not good to transfer the parameters of the geometry using global variables. Therefore we include an abstract pointer to "user data" as the rst parameter into each function call.
3.9 conclusion
So we have dened a new interface for the geometry description. It consists of the region-function, three boundary intersection functions and nine test functions. In our implementation we have included also service functions, for example a free-function to release the storage occupied by the parameters if the geometry is no longer necessary: Other useful functions may be save and load (currently not implemented). There is a lot of information usually available for the grid generator which is not included into the interface and therefore not available for our grid generator:
void Free(void *user-data);
the number and lists of identiers of the regions, faces, lines and vertices. explicit discretizations of the faces and lines, coordinates of the vertices. other derived information, for example curvature and so on.
12
4 realizations of the interface
Let's consider now possibilities to realize this interface for dierent input data. Some of these realizations we have already considered in the previous considerations:
We have considered the approximation of the interface if only a regionfunction is given. We have considered the case of a complete grid. In this case the neighbourhood search algorithm seams to be the best way to realize the region function. There are also fast and exact realizations of the three intersection functions. We have considered the usage of characteristic functions, especially as a method to change a given geometry. Characteristic functions can also be used to subdivide boundary faces and lines. This method can be used to avoid the "rounding o" of lines and corners, because after the subdivision of a face or edge the exact position of the edge or corner dividing the face resp. edge can be computed by the grid generator using the edge resp. corner function. A trivial geometry (only one region, no boundary) can be used as a start for further changes. It is easy to dene "boolean operations" (union, intersection) not only for regions, but also for boundary faces and lines. It is possible to use graphical input data (pixel maps) to dene a regionfunction. For time-dependent problems the following algorithm can be used: Start with the grid of the previous time step. Then compute the characteristic function (or necessary data for the computation of this function) on this grid. Then use the region-function for a geometry dened by a grid to interpolate the function characteristic for the nodes of the new grid. At last make the change-operation using this characteristic function. 13
To guarantee that this algorithm works the old region-function will be called in the change-operation before the characteristic function will be evaluated. It is also possible to realize this interface for the "standard" geometry descriptions. If the geometry is described by a complex boundary grid fast search algorithms have to be used to make the realization ecient. For example, search trees (quadtree/octree) can be created. Another possibility is to create the Delaunay grid containing the boundary nodes and use the neighbourhood search algorithm described before.
We see, that it is possible to realize the interface for all usual types of input data for a geometry description. With this interface it is possible to create even innite geometries (for example Julia sets)!
4.1 How to use this interface
Because a lot of possible information about the geometry is not visible for the grid generator it is not possible to guarantee the topological correctness of the resulting grid. This can be easily proved considering the case of a very small enclave. Since there is no information about the existence of this enclave and there may be no grid point inside also in the nest grid we cannot even detect that we have made an error. But there is a class of geometries there it is possible to construct the geometry without errors. These are geometries where all objects are convex:
Theorem 1 Assume we have a geometry with a nite number of objects
(regions, faces, lines and vertices) so that all these objects are closed convex sets, no interior point of an object 1 lies in another object of the same or higher dimension, and there is an exact implementation of the interface for this geometry. Then there is an algorithm which can reconstruct the geometry in a given bounded region exactly using only the rst four functions.
1
Interior point for an object means interior in the subset topology of the union of all objects with the same dimension
14
The principal way the algorithm works is clear: In the rst step we dene the region number of all nodes of some grid containing the region of interest using the region-function. Then there can be dened lines with dierent region numbers for the end points. For these lines we use the face-function to dene the rst boundary face intersection. Using the inverse order for the end points we can proof that there is only one intersection. Else we make further renement of the edge. This process stops because of the nite number of faces and the convexity (so there cannot be two intersections of the edge with the same face). Now we consider the sides. For every intersection of their border with a boundary face we use the line-function to nd the continuation. If there are dierent inner intersections with boundary lines we use further renement of the side. Then for every inner intersection point we subdivide the side connecting the point all boundary face nodes on the border of the side by a line. If we have two boundary face nodes of the same face on the border we also subdivide the side connecting these points by a line. At last we do the same for elements and the vertex-function. 2 This theorem shows a way we can avoid errors: We can subdivide objects so that the resulting topology is more convex. For example, if we have a thin layer inside a region, this layer may not be found. Subdividing the region into two part (over and under the layer) makes it possible to detect the layer. Subdividing a face along a sharp edge or an edge at a corner makes it possible to nd the exact position of the edge or the corner. In general it is not possible to avoid geometrical and topological errors. But the topological errors are always local, and they can be reduced by further renement. For many applications this may be not good enough, because also a little topological error can have global consequences. But in the applications we have in mind | the simulation of processes which may change the topology | it is clear that the topology is dened only modulo the accuracy of the simulation itself. So for these applications our concept seams to be even saver than the standard method, where every error can have global consequences.
For an accurate proof we need a more accurate denition of the geometry and the interface, especially we have to consider degenerate cases in detail.
2
15
5 Grid Quality Criteria
There are a lot of dierent criteria to measure the quality of a grid. This seams to be natural, because there are dierent applications and discretizations, and a grid may be good for one but bad for another case. So for dierent equations and discretizations there have to be dierent quality measures for the grid. From this point of view it makes no sense to introduce a criterion independent of an equation or discretization. Consider, for example, the case of the diusion equation with nonconstant diusion coecient. The function and the diusion coecient are dened on the nodes of the grid. Now we have dierent possibilities for the discretization. Let's consider the following: 1. the FEM discretization (with lumping of the mass matrix). 2. the FV (nite volume or box) discretization based on the Voronoi boxes, the diusion coecient will be considered as constant on elements, it's value computed by averaging over the nodes of the element. 3. the same FV method, but now the diusion coecient will be considered as constant on the inuence region of the edges and computed by averaging over the two nodes of the edge. For the diusion problem we have to avoid negative concentrations. So a stable discretization is necessary. In this case it is possible to nd conditions so that there will be no negative concentrations (they are sucient for the M-matrix property of the resulting discretization matrix). Consider at rst the 2D case. For our three schemes we obtain the following conditions: 1. a grid without obtuse angles. 2. a grid without obtuse angles. 3. a Delaunay grid without obtuse angles opposite to a boundary edge. So, in 2D for a given node set the Delaunay grid will be optimal for all of these methods. In the 3D case we have the follwoing picture: 16
1. For the FEM discretization we obtain a criterion which may be fullled for non-Delaunay grids and not fullled for the correspondent Delaunay grid. 2. For the FV method we need a Delaunay grid so that the centre of the minimal ball containing the element must be inside the element. 3. In this variant we need a Delaunay grid so that the centre of the minimal ball containing the element must be inside the region. From point of view of grid generation it is easy to guarantee only the Delaunay property, because there is a standard algorithm to construct a Delaunay grid for a given point set. There is no algorithm known to nd the optimal grid for the FEM method for a given point set in 3D. The standard method is to start with the Delaunay grid and to use local transformation and additional point inclusion. An extreme example for the dierence between these criteria is the so-called sliver | a highly degenerated tetrahedron, but without obtuse angles on it's sides. Such slivers are very bad from point of view of the FEM criterion, but they may be optimal from the Delaunay criterion. In our grid generation algorithm we create a Delaunay grid and try to avoid situations which are bad especially from point of view of the FV method, because the chance to get a good grid from this point of view is much higher. Slivers can occur in our grids, and we recommend to use FV methods instead of destroying the slivers.
6 Anisotropy
Anisotropic grid generation techniques are useful to minimize the number of nodes. There are applications where anisotropic grid generation is obviously necessary because of the dierent length order in dierent directions. For example in environmental processes we have often dierent scales in horizontal and vertical direction. Another example are thin boundary layers in viscous uids. But anisotropic renement is useful not only for such special applications. In principle, anisotropic renement can be used every time. In the worst case there may be no reduction of the node number. But this worst case 17
not often appears. Usually in almost every application you have locally a "gradient direction" which has to be rened better than the orthogonal directions. The other point is that the reduction eect is much greater in 3D, because we have usually two "orthogonal" directions which must not be rened. If we obtain a factor 2 in a 2D problem, we may obtain a factor 4 in the similar 3D problem. So anisotropic techniques seams to be useful for every 3D application.
grid quality criteria for anisotropic grids: As in general, in the case of
anisotropic renement we have to consider our problem and our discretization to formulate the correct criterion. For example, consider a problem with anisotropic diusion coecients. Then the optimal grid may be obtained by replacing the Delaunay criterion by an anisotropic analogon | it depends on the discretization we want to use. But usually the anisotropy is not caused by the equation and discretization, but by the special situation (geometry, initial values, boundary conditions) we have to consider. So the quality measure we have to use coincide with isotropic case. So we create a Delaunay grid and try to fulll the additional conditions for the FV criteria also in the anisotropic situation.
6.1 problems with anisotropic Delaunay grids
There are a lot of specic problems in anisotropic Delaunay grid generation. The reason is the following alignment property: In a highly anisotropic Delaunay grid the nodes have to be aligned in the direction of "greatest renement". This is easy to see. In a highly degenerated triangle the only possibility to avoid big obtuse angles is to have two nearly right angles. But this means that the shortest edge lies in the direction of high renement. The same consideration can be made also in 3D. This eect creates a lot of problems, because we cannot use standard methods which do not preserve the alignment property:
point insertion: inserting a point into a Delaunay grid is a standard
method to solve dierent problems:
local renement boundary correction (divide edges going through the boundary) 18
destroy slivers.
But if the new point is not "aligned" the resulting grid is very bad. We obtain elements which are too big or non-Delaunay elements. node shifting: this method often will be used near the boundary | an inner node near the boundary can easily create bad elements near the boundary, and so it may be better to shift this node to the boundary. But here we also have to consider the alignment. So we can shift the node only in alignment direction. grid smoothing: standard grid smoothing procedures cannot be used, because nodes may be shifted away from the "alignment curve". We see, that our possibilities to manage an anisotropic Delaunay grid are restricted. We have to modify these methods so that they preserve the alignment or we cannot use them.
7 the algorithm
The main steps of the algorithm are
An octree-like anisotropic renement controlled by application dependent criteria. A boundary shift procedure. Intersections of the boundary with the octree will be computed, nodes will be shifted to these intersections. A Delaunay step computes the resulting Delaunay grid for the point set created by the previous steps. A boundary correction step detects and corrects incompatibilities of the resulting grid with the geometry description.
Let us consider now these steps in detail.
19
7.1 renement
We start with a modication of the anisotropic quadtree/octree renement procedure. This is a standard algorithm starting with one element (a quader in 2D, a cuboid in 3D). Then external criteria will be used to rene the elements. A tree of elements will be created. In the isotropical variant, a quader will be splitted into four and a cuboid into eight parts. In the anisotropic variant an anisotropic criterion is necessary, telling in which direction the element has to be splitted. The element will be splitted into two parts. Octree techniques are considered for example in [1], [2], the anisotropic variant in [5]. The main dierence between the standard quadtree/octree method and our method is that we rene not the elements but their edges. This is a little bit more exible. So in the standard octree method in the highest level we have at minimum four nodes, in our method there may be one. We also use a dierent data structure. Our structure is based on nodes. For every node we have the six (four in 2D) neighbour nodes in the orthogonal directions (minimum two of them must be dened). The cuboids play only a secondary role, and we don't use a tree structure. We have included them into our data structure only to make some search operations faster. In 2D we don't use the quaders. The renement criteria we use are also based not on elements, but on nodes and edges. We use the following criteria:
a application dependent function dening the maximal length of an edge around a given point (isotropical criterion). a application dependent function dening if a given edge has to be rened (anisotropic criterion). parameter-controlled regularization criteria.
Some minimal regularization property is necessary | direct neighbour edges in the same direction must have the same or the neighbour renement level. This is automatically guaranteed by the edge renement procedure | a neighbour edge of a coarser level will be rened by a recursive call of the renement procedure before the edge itself will be rened. After this renement the other parameter-controlled regularization criteria will be 20
tested, and some of the neighbours with the same renement level as the rened edge may be also rened. The time complexity of the renement algorithm is Because we rene only edges in coordinate directions we obtain automatically an alignment in these directions. If one of these directions approximately coincides with the "highest renement" direction we can obtain very anisotropic grids. For a skew highest renement direction the resulting grid will be more isotropic. So, if the highest renement directions are equally distributed you will have parts with isotropical renement and parts with high anisotropic renement. But also in this case we can obtain approximately a factor 2 in the 2D case and factor 3 in 3D for the number of nodes compared with pure isotropical renement.
O(n).
7.2 boundary shift
In the rst step only inner nodes of the regions have been created. Now we have to create the boundary nodes. One possibility to do this (for example [3]) is to compute the intersections of the edges with the boundary faces and of the rectangle sides with boundary lines and the cuboids containing boundary vertices and to include these points into the grid. 3 We use the following modication of this method: instead of including a new node at the intersection we shift an existing grid node to this place, and we don't consider all of the intersections. This gives us more freedom for manipulation to o...