|
| Two-dimensional finite-element mesh generation using regular triangular grids and hierarchical model definitions |
W. Robert J. Funnell
Departments of Biomedical Engineering & Otolaryngology
McGill University
3775, rue University
Montréeal, Quéebec H3A 2B4, Canada
| This page is a conversion to HTML of a manuscript originally written by me in 1987/1988. It was revised in 1990, but never finalized and never published. |
SummaryThis paper describes a method for generating a two-dimensional mesh of triangles over an irregular polygonal region. Using a specified mesh resolution, a regular grid of equilateral triangles is laid down over the interior of the region. The remaining strip around the periphery of the region is then analyzed and triangulated. The method can handle complex regions with holes and very narrow necks. A relaxation technique is used to make the resulting triangle shapes more uniform. A hierarchical modelling procedure is also described for creating complete finite-element models of complex shell structures. Model boundaries are defined using straight-line segments, circular arcs, and splines, defined by control points. The distinct surfaces of the model are defined in terms of boundary segments, with associated material properties. Three-dimensional curvatures may be specified. |
Many techniques have been described for dividing an irregularly shaped two-dimensional polygonal region into elementary (generally triangular or quadrilateral) subregions, whether for use as a finite-element mesh or for the purposes of interpolating or representing data defined at points in the region. Several reviews have been published{6807,3992,6811,7151,7764}. In this section I shall briefly summarize previous work and then introduce a method which we have been using to generate finite-element meshes for models of the middle ear{7791}. This method is particulary useful for generating meshes consisting of multiple complex regions with shared boundaries.
Methods may be divided according to whether they (1) do not involve internal nodes at all; (2) use prespecified internal nodes; (3) generate new internal nodes as a side effect of subdividing the region into elements; or (4) generate new internal nodes first and then generate elements using both the original boundary nodes and the new internal nodes. In this last group of methods there may be more or less information available to guide the generation of elements, depending on how the internal nodes are generated. The methods may be subdivided according to whether the elements are generated (4a) completely regularly, (4b) completely irregularly, or (4c) regularly for the bulk of the elements but irregularly for those at boundaries or in transition regions. In group 4a, the internal nodes are generated in such a regular pattern that the element definitions follow directly; in effect, the nodes and elements are defined simultaneously. The internal nodes are generated from the boundary specification either by solving partial differential equations or by algebraic interpolation. In group 4b, internal nodes are generated using some strategy to obtain a reasonably uniform distribution, and then nodes are joined together to form elements using little or no a priori information about the way the nodes were generated. Once the internal nodes have been generated, these methods face essentially the same task as those in Class 2 above. (Within this group, note that the paper by Suhara & Fukuda{6806} has often{6808,7810,7779} been cited with the authors' names in the wrong order. Lo{7779} also cites Gordon & Hall{7870} incorrectly.)
In group 4c, internal nodes are generated in some regular pattern and then that regularity is used as much as possible in joining them to form a regular pattern of elements. The regularity will in general break down at boundaries or in transition regions, and then some sort of method must be used for finishing off the mesh. The regular grid has usually been taken to be rectangular, either geometrically uniform{7867,7151,7056}, topologically uniform{7872}, or recursive (quadtrees){6800,7834,6903}. (Note that Ho-Le{7764} states that Kela et al.{7834} use a uniform rectangular mesh, whereas they actually use quadtrees. Recursive fractal subdivision using the Pascal-Sierpinski gasket has also been suggested{7865}, but it does not provide conformity between elements{8331}.)
Thacker et al. (1980){4292}, on the other hand, use a triangular regular grid. They start by superimposing a uniform mesh of equilateral triangles over the region to be meshed. They then identify a zig-zag approximation to the boundary of the region by selecting those nodes (on the regular triangular grid) which lie closest to the region boundary. (Note that Kikuchi{7151} and Ho-Le{7764} both mistakenly state that Thacker et al. select only grid points which lie inside the boundary.) Finally, the nodes on this zig-zag boundary are moved so as to lie upon the region boundary. An iterative relaxation technique is used to adjust the positions of both the internal nodes and boundary nodes so as to optimize the element shapes.
This method is very attractive, since it triangulates most of the region very quickly, and with well shaped triangles, by simply laying down a regular uniform grid. In this paper I shall describe a similar technique. The main difference is that the boundary nodes are determined at the beginning and then fixed. This is a marked advantage when generating a mesh for a structure consisting of multiple regions with shared boundaries, since the boundaries can be defined once and then the different regions be meshed independently. As discussed below, this also makes it feasible to extend the method to the generation of matching meshes for consecutive slices of a three-dimensional structure. Fixing the boundary nodes also simplifies the final relaxation phase.
In the following Sections I shall outline the main stages of this method, namely (1) fixing the boundary nodes; (2) laying down the regular mesh of equilateral triangles in the interior; (3) generating triangles linking the boundary nodes to the interior mesh; and (4) using relaxation to improve the triangle shapes. In Section VI a hierarchical modelling scheme will be described which takes advantage of the ability of this triangulation procedure to create meshes for complex structures.
The methods described in this paper have been implemented in FORTRAN. The same source files are used under the VMS™ operating system on DEC™ VAX™'s, and under the Microsoft® MS™-DOS operating system on IBM®-PC-compatible machines.
To facilitate the generation of compatible meshes for adjoining regions, it is desirable to fix the boundary nodes of the meshes using only boundary information, independent of the interior mesh. A simple way to do this is to generate nodes that divide the boundary into equal intervals, the size of the interval depending on the refinement of the desired mesh. This degree of refinement may be defined, for example, in terms of a nominal number of elements across the overall structure. If some measure of the 'diameter' of the structure is defined, say d = max(xmax-xmin, ymax-ymin), and the mesh resolution is given as m elements/diameter, then the nominal element size would be s = d/m. If the perimeter is p then the number of evenly spaced nodes to be generated is n = int(p/s) + 1, and the interval between the nodes is s' = p/n (where the function int(·) returns the largest integer ≤ the argument). The first of these new nodes is put at the beginning of the boundary. Each of the subsequent new nodes, i = 2 to n, is then located by moving along the specified boundary, line segment by line segment, until the distance traversed since the previous point is s' or greater, and then doing a linear interpolation to locate the new point within the current line segment.
This procedure results in nodes which are equal distances apart if the distances are measured along the original boundary, but which may be unequal distances apart if the distances are measured along straight lines joining the new nodes. The distances may be extremely nonuniform if the boundary has sharp corners and the mesh is coarse, as shown in Figure 1a. It would be feasible to implement a more sophisticated algorithm to deal with this case. In practice, however, it is generally not difficult to solve the problem by manually dividing the original boundary specification into segments such that sharp corners occur at the beginnings and ends of segments, as in Figure 1b.
In cases where the initial contour definition arises from some sort of tracing and digitization, manual or otherwise, it may be desirable to smooth the contour before generating the evenly spaced boundary nodes. We have used a simple three-point median filter for this purpose.
It is assumed above that the boundary of the region to be meshed is specified by a sequence of straight-line segments. As discussed in Section VI below, it may be desirable to also use arcs and splines or other curves to specify the boundary. However, these curves may always be reduced to sequences of straight-line segments to whatever resolution is desired. It would also be possible to extend the algorithm described above to handle curved arc-lengths directly, if desired.
To simplify subsequent processing, it is required at this stage that the boundary nodes be numbered in some specified direction, in our case counterclockwise. It is not difficult to reverse the node numbering of the given boundary if necessary.
The basic approach to generating the regular mesh over the
interior of the region is to establish a uniform grid of equilateral
triangles that overlays the region, and then to identify those nodes
which lie inside the boundary. The nodes lie in horizontal rows,
and in diagonal columns at an angle of 60° from the horizontal.
The left-hand and bottom edges of the mesh may be taken to be tangent
to the boundary.
Thus, the vertical position of the bottom row
is simply ymin. The horizontal position of the leftmost
diagonal column is found by effectively rotating the boundary 30°
counterclockwise, finding the coordinates of the leftmost point, and
then rotating back to the original position. Thus, the leftmost column
must pass through the point (x1,y1) where
| x1 = u cos 30 + v sin 30 |
| y1 = -u sin 30 + v cos 30, |
| u = min(xi cos 30 - yi sin 30) |
| v = min(xi sin 30 + yi cos 30). |
| x0 = x1 - (y1 - ymin)/√3 |
| y0 = ymin. |
| x2 = u cos 30 + v sin 30 |
| y2 = -u sin 30 + v cos 30, |
| u = max(xi cos 30 - yi sin 30) |
| v = max(xi sin 30 + yi cos 30). |
| x3 = x2 - (y2-ymin)/√3 |
| y3 = ymin. |
The length of the sides of the triangles
is given by s' as calculated in Section II above.
The vertical distance between
a row of nodes and the one above it is δv = s'√3/2,
and the nodes in a row are offset horizontally from the ones below
by δh = s'/2.
If the nodes are labelled (i,j), where j is the row number counting
from the bottom, and i is the number of the node within its row,
counting from the left, then
| xi,j = x0 + (j-1)δh + (i-1)s' |
| yi,j = y0 + (j-1)δv. |
| imax = int((x3-x0)/s') + 2 |
| jmax = int((ymax-ymin)/δv) + 2. |
We now wish to scan over all of the nodes in the regular grid, excluding those outside the boundary. However, it is not desirable to retain those nodes which are just barely inside the boundary, since this may later lead to very badly shaped elements. Instead we retain a node only if its distance inside the boundary exceeds some threshold (Figure 2). Once a node is determined to be within the boundary, its distance from the boundary can be determined by dropping perpendiculars from it to each line segment of the boundary. For each line segment, if the perpendicular intersects the segment within its length then we take the length of the perpendicular as the distance from the boundary; otherwise we use the distance to the closer of the two ends of the segment. The shortest such distance over all of the boundary segments is the distance of the node from the boundary.
A reasonable value for the distance threshold is s'/2. In this case the nodes retained will generally have distances from the boundary in the range of s'/2 to 3s'/2.
Once we have decided which nodes to accept as belonging to the regular interior grid, we must create triangular elements by defining edges joining the nodes. This is straightforward because of the regularity of the grid. For each node (i,j) belonging to the grid, we first consider the node above it and to the right, that is, node (i,j+1) (see Figure 3). If this node does not belong to the grid we do not generate a new triangle. If it does belong to the grid, then we check whether the nodes (i-1,j+1) and (i+1,j) belong to the grid, and generate a triangle for whichever one does, or for both or neither. It is very important that the vertices of the created triangles be numbered consistently, that is, in our case, always counterclockwise. The resulting element mesh is shown by the solid lines in Figure 4. If quadrilateral elements were preferred, it would be very easy to generate a single quadrilateral whenever all four of the nodes (i,j), (i,j+1), (i-1,j+1) and (i+1,j) belong to the grid.
Every time an equilateral triangle is created, its three edges are added to a list of edges (oriented as actually used in the triangle) and the three nodes are marked as having been used. Note that it is possible to have nodes which were retained as part of the interior mesh but which end up not being used in any of the equilateral triangles. This and other anomalies will be considered below. It may also happen that no internal nodes are generated at all. In this case we can use the same triangulation algorithm as given below for filling in the final gaps in the mesh.
Once the interior mesh has been generated, we are left with an untriangulated strip between it and the external boundary. Because of the distance criterion used above in deciding which nodes to keep as part of the interior mesh, the distance between the two boundaries is such that it is reasonable to plan to interpose a single layer of triangles between them (except possibly at sharply convex corners of the external boundary). Thus, each straight line segment on each boundary will form one edge of a triangle, with the other two edges normally being formed by joining its two ends to some one node on the other boundary. This is similar to the situation in the method of Heighway & Biddlecombe{7872}. Their approach is to use the triangle edges around the outside of the interior mesh as the bases of new triangles, using a 'best' node on the external boundary as the third vertex in each case; the newly created triangle edges are in turn used as bases for new triangles, and so on, until the triangulation is complete. The procedure may be visualized as gradually expanding all parts of the internal boundary out towards the external boundary. For every new triangle, choosing the 'best' external node requires comparing all external nodes.
The algorithm used here, on the other hand, works by marching around the gap between the internal and external boundaries, creating new triangles and filling in the gap as it goes. Once a starting point has been found by comparing all external nodes, creating each new triangle requires consideration only of the next node on each of the internal and external boundaries, which speeds things up. If the next triangle would be too badly shaped, however, a gap is left and a new starting point is chosen. A separate algorithm fills in the gaps later.
Note that reasonably smooth contours are very easy to triangulate, but exotic cases require much more care. Contours that are well behaved at high mesh resolutions will often develop problems at low resolutions. The exhaustive geometrical checks required to guarantee a valid mesh are computationally expensive, and in some applications it may be desirable to apply them only in cases where the initial quickly generated mesh is corrupted in some way{7872}.
The first step in triangulating the gap is to identify the boundary of the interior mesh. This can be done using the list, mentioned above, of all edges of the interior equilateral triangles. If we sort the list in order of the smaller of the two node numbers for each edge, then edges which appear twice will appear consecutively in the sorted list and can easily be identified. If an edge appears twice it is because it corresponds to a shared edge between two triangles, and so it cannot lie on the boundary of the internal mesh. The two instances of the edge must be oriented oppositely because of the consistent counterclockwise numbering of the triangle vertices; this provides a consistency check.
Edges which appear only once are thus identified as lying on the boundary of the interior mesh. These edges can then be linked together in sequence to form a boundary specification. We choose the first available unduplicated edge, and mark it as having been used. We then look through all the remaining unused unduplicated edges, looking for ones which start from the end of the current one. If only one is found, then it is the next edge in the boundary. If more than one is found we must choose among them in a way that will produce a counterclockwise boundary that is not self-intersecting. We can do this by choosing the edge which is encountered first when sweeping counterclockwise in an arc starting at the current edge and centred at the common node (Figure ??). The calculation of the relevant sweep angle is described in Appendix ??.
Taking the chosen edge as the current edge we then repeat the process, continuing until we return to the starting node of the first edge. Now it is quite possible that there will still remain unused unduplicated edges at this point. If so we choose the next available one and start the process again. This continues until we have some number of unconnected counterclockwise internal boundaries (Figure ??).
There may be some isolated internal nodes which were not used in any equilateral triangle, and thus are not part of an internal boundary (Figure ??). If so, we can take each of them in turn and check if there is some internal-boundary node which is one of its neighbours on the original regular grid. If found, we can insert the isolated node into that internal boundary; the new edge is added twice to the list of edges, once in each orientation. Since an isolated node may be joinable to some boundary only through the intermediary of some other isolated node or nodes, we must keep iterating until either no more isolated nodes exist at all, or no isolated node can be joined to any internal boundary. If any isolated nodes still exist after all this, one of them can be considered as a degenerate one-node internal boundary, and the whole process repeated. If there had been no internal boundaries to start with but there were isolated internal nodes, then we would have started by considering one of them as a one-node boundary.
In the simplest case, there is a single internal boundary and the gap between it and the external boundary consists of a single continuous strip. The basic approach used here involves choosing a starting pair of internal and external nodes. Two candidate triangles are then compared, corresponding to advancing by one node along either the internal or the external boundary (Figure ??). The better triangle is chosen and its edge is then used as the basis for the next choice. The process is repeated until we arrive back at the starting point. As each triangle is created its edges are added to the list of edges.
This approach is based on a common method for triangulating a three-dimensional surface between two contours lying in different planes. There are two classes of such triangulation algorithms{7792}: in one class, the algorithm passes along the boundaries once, making irrevocable decisions about which nodes to join using local information based on some sort of cost function. The other class of algorithms determines a globally optimal solution by effectively comparing all possible triangulations, again using some cost function. This cost function is designed to favour triangles which are well shaped; that is, they should be as nearly equilateral as possible. The global algorithms are more robust but considerably more expensive computationally. In the present application we have found a local algorithm to be adequate.
The cost function used is simply the length of the new edge to be created between the internal and external boundaries: the triangle with the shorter edge will be better shaped. If the length of the edge is greater than 1.5s', it is considered to be unacceptable. In addition to the length calculation, each candidate edge is checked for topological correctness. First we check that the internal end of the edge lies to the left of the neighbouring segments of the external boundary; this requirement follows from the counterclockwise orientations of both the internal and external boundaries, and is checked using the sweep angle described in Appendix ??. Next we check that the external end of the edge lies to the right of the internal boundary. ((N.B. 2pi --> 0!!). Finally, we check that the edge does not intersect either boundary, nor any of the edges that have been created so far. ((details of intersection check??))
There may arise cases where neither of the candidate triangles is acceptable, for either geometric or topological reasons. In this case the triangulation sequence is stopped, a gap is left, and another starting point is sought. Using the current internal-boundary point, we search forward along the external boundary looking for an acceptable point. If none is found, we try again with the next internal-boundary point, and so on. (We do not try re-using the current internal-boundary point if it is the original starting point for this internal boundary, to avoid a rare situation which may prematurely terminate the triangulation sequence.)
In cases where there are multiple internal boundaries, the above method is applied to each one in turn.
!!special testing; N.B. Thacker avoided coarse meshes, quit if pathological case arose; also Heighway & Biddlecombe, 1982 [7872]
The above procedure may result in untriangulated gaps between the internal and external boundaries. These gaps can be identified using the list of edges. First we add all segments of the external boundary to the list of edges, with their orientations reversed. Thus any external edge which has already been used in a triangle will be in the list twice, with opposite orientations, and those which have not been used yet will be in the list once. We now use the same procedure as described above for identifying the internal boundaries, but with the edges reversed and using a clockwise sweep. This produces a counterclockwise boundary for each gap, within which we wish to produce triangles without generating further internal nodes (corresponding to the first category of triangulation algorithms defined in the Introduction). The approach used here is to look for the boundary vertex with the most acute convex angle and to create a new triangle using that vertex and its two neighbours. The boundary is modified by removing the new triangle, and the procedure is repeated until the gap has been completely triangulated.
Once the mesh is complete, a check on its topological consistency can be made by calculating the areas of all of the triangles and adding them, and comparing the sum with the calculated area of the external boundary. Since the sign of the area calculated for a polygon depends on its orientation, this procedure confirms that all of the triangles are oriented counterclockwise as well as checking that no gaps remain in the mesh and that no spurious superimposed triangles exist.
Once the elements have been defined, it is useful to adjust the positions of the internal nodes slightly in order to make the triangles around the border as close to equilateral as possible. A simple and well known way to do this is to require that each node be positioned so that its x and y coordinates are the means of the x and y coordinates of those nodes to which it is connected. This method has been used by Kamel & Eisenstein (1971){6805} and by many others since. Figure 8 shows the mesh of Figure ?? after such a relaxation. (Alternative relaxation criteria{7799} can be used if necessary but this simple one has been found to be quite satisfactory.)
This relaxation has generally been carried out iteratively. It can be formulated as a set of linear simultaneous algebraic equations, or a matrix equation, with the unknowns being the positions of the internal nodes, the vector on the right-hand side being defined in terms of the positions of the boundary nodes, and the matrix coefficients depending on the connectivity of the mesh. This matrix equation can be solved directly. However, for a typical mesh the matrix is large and a direct solution is computationally demanding. The solution can of course be accelerated by taking advantage of the symmetry, sparseness and bandedness of the matrix, as done for the finite-element solutions themselves, but this introduces considerable complications. The iterative approach is made particularly attractive in this case by the fact that it is not necessary to position the nodes with great precision, since the aim is simply to make the element shapes reasonably uniform. In practice it generally takes only five or ten iterations to position the nodes to within 0.01s'. With the present method of mesh generation the bulk of the elements are already perfectly equilateral, and only those at or near the boundary need to be adjusted. It might be feasible to speed up the iterative relaxation by not treating nodes deep in the interior of the regular mesh, but this additional complication has not so far been found necessary.
The mesh-generation technique described above has been embedded in a modelling scheme for creating complete finite-element models of complex shell structures. The geometry of the model is defined in a hierarchical way{6805,6809}. The distinct surfaces of the model are defined in terms of named boundary segments. These boundary segments are in turn defined using straight-line segments, circular arcs, and splines, themselves defined in terms of named control points. (Note in passing that the methods of Kamel & Eisenstein (1971){6805} for boundary definition, in terms of straight lines, circular arcs and parabolic segments, were redescribed by Tsybenko et al. (1980){4774} without citing their source.) It should be emphasized that the division of the structure into different regions is not done just for the convenience of the triangulation algorithm, as is the case for some mesh-generation schemes. Rather, this division is done only as required for the specification of variations in material properties, thicknesses, loads, and so on.
Figure ? shows an example of a model consisting of two regions with different material properties. The definition file for this model is shown in Figure ?.
The control points are defined by specifying coordinates. Any
coordinate may be following by the letter R, indicating that
it is to be taken as relative to the corresponding coordinate of the
previous point. Thus, for example, point L7 in Figure ? has
a coordinate specification of (3.55,2.1R), meaning that it
has an absolute x position of 3.55 and a y position
2.1 units above that of the previous point, L6. (Note that
in Figure 8 all the points have z coordinates of 0.) The
ability to mix absolute and relative values for the coordinates of a
single point often greatly simplifies the design and modification of a
mesh definition.
The lines comprising the surface boundaries are defined
by linking control points
using the commands M, L, S, CL or CC,
standing for Move, Line, Spline, CLockwise arc
and CounterClockwise arc, respectively.
For example, in Figure 8,
line LEFT1 consists of a straight-line segment
from C2 to L2; followed by a circular arc from L2 to
L4 with its centre at L3; followed by three more straight-line
segments. Line CENTRE consists of a single spline segment defined
by the points R8, C1, C2 and R1: the actual line is
between C1 and C2, with R8 and R1 serving as control
points to define the curvature of the quasi-cubic spline{4648}.
Splines are particularly useful in defining meshes for naturally
occurring, irregularly shaped structures.
The definition for each line
is accompanied by a set of six boundary conditions which apply along
that line.
These boundary conditions correspond to the
three displacements ux, uy and uz
at each node, and
the nodal rotations about the x, y and z axes. Each
boundary condition may be either F, C, S or M, depending
on whether that degree of freedom is Free, Constrained to zero,
Slaved to a master degree of freedom, or is itself a Master.
In the example, the nodes on the outside boundaries
are all specified as fully clamped (CCCCCC)
while the nodes on the hole boundaries
are completely free (FFFFFF).
Each of the two surfaces is defined by concatenating line segments,
with a minus sign indicating that the sense of a line is to be reversed.
The outer boundaries of the two surfaces, LEFT
and RIGHT, are defined by the lines LEFT1, LEFT2, RIGHT
and CENTRE, with CENTRE being shared.
The holes are defined by the lines HOLE__L and HOLE__R;
lines JOIN__L and JOIN__R serve to join the holes to the
outside boundaries.
Boundary conditions are specified for the internal nodes that will be generated for each surface. The available boundary condition codes are the same as described above for lines.
A surface definition may also include a material-type designator
(e.g., MAT:1)
and a thickness value (e.g., THICK:0.1),
as well as control information for generating
three-dimensional curvatures as described below.
The final part of the model definition includes various non-geometric parameters. In the example, an applied pressure of 28.28 dyn cm-2 is specified, as well as two sets of material properties (including density, Young's modulus and Poisson's ratio).
The mesh generator described here has been used to generate meshes for shallow curved shells by simply generating a two-dimensional mesh, treating it as the projection of the desired three-dimensional mesh onto the x-y plane, and then determining z coordinates for each of the mesh nodes in some way. (For example, Figure 9 shows a finite-element mesh generated with the present method to model the eardrum of the cat{7791}. The z coordinates for the internal nodes of this mesh were obtained by interpolating between the specified z coordinates of the boundary nodes using circular arcs perpendicular to the x-y plane.) This approach works as long as the shell is not too deep. More general approaches to generating three-dimensional shell meshes{6801,7965} could also be used with the triangulation technique described here.
Once the model has been defined, the triangular mesh can be generated for any desired resolution. Figure 10 shows meshes generated for the same model as shown in Figure 6, for two different resolutions. The resolution is specified in terms of a nominal number of elements per diameter as described in Section II above. The mesh in Figure 6 has a resolution of 12 elements/diameter; those in Figure 10 have resolutions of 6 and 30.
The nodes of the completed mesh are renumbered in order to reduce the bandwidth of the system equation{4840}. The mesh is output, complete with material properties and so on, in a form directly usable by the SAP IV finite-element programme{7789}.
As an indication of the computer time involved in the mesh generation process, the mesh of Figure 6 with 12 elements/diameter required, when generated on a VAXstation™ 2000, roughly 10 seconds of elasped time to process the model definition and generate evenly spaced boundary nodes; another 10 seconds to do the triangulation; 3 seconds for the relaxation; and about 35 seconds for the bandwidth minimization and the production of the final finite-element problem-definition file. For the mesh of Figure 10 with 30 elements/diameter, the corresponding times were approximately 15, 80, 25 and 140 seconds. Note that these figures include times for extensive graphical display of the steps of the meshing process, as well as file reading and writing.
The meshes described here have a uniform element size over the entire model. A number of approaches are possible if nonuniform meshes are required. The nonuniformity may be introduced as an intrinsic part of the mesh generation procedure{6805,6796}. In this case the nonuniformity is generally specified in terms of desired node densities along the boundaries. Alternatively, the nonuniformity may be introduced after a uniform mesh has been generated. If this nonuniformity is specified in terms of some criterion that can be computed on an element-by-element or node-by-node basis, either as a function of the problem definition{4292} or adaptively as a function of the results of the finite-element analysis{6171}, then it would be easy to modify the meshes generated by the present method by adding or deleting elements or nodes. If, on the other hand, the desired nonuniformity were specified in terms of node densities along the boundaries, it would be necessary to perform some sort of interpolation between the boundaries in order to know how to modify the existing mesh.
A method for the automatic generation of two-dimensional meshes of triangles has been presented. The method is a modification of the method of Thacker et al. (1980){4292}, which works by laying a regular mesh of equilateral triangles over the region to be treated and then adjusting it. The attractive feature of this method is that most of the surface is triangulated very quickly and efficiently with ideally shaped elements. The present method differs from that of Thacker et al. in that the boundary nodes are fixed before the triangulation is done. This has two advantages. Firstly, the iterative relaxation scheme for adjusting the element shapes is greatly simplified if the boundary nodes do not need to be adjusted. Secondly, and more importantly, fixing the boundary nodes makes it very easy to independently create compatible meshes for separate regions which share boundaries, thus facilitating the creation of meshes for complex structures. It also leads to the possibility of using a series of separately generated two-dimensional meshes as the basis for a three-dimensional mesh of tetrahedra{7790,7957}.
A method is described for embedding the mesh-generation algorithm in a hierarchical model-definition system. Noteworthy features include the ability to mix absolute and relative coordinate values; the specification of mesh resolution independent of the model definition; the use of splines for boundary definitions; and the integration of finite-element modelling parameters such as boundary conditions, loads and material characteristics.
This work was supported by the Medical Research Council of Canada.
Given the directed line segment AB and a point C (Figure ??), we wish to calculate the angle θ swept out by rotating the line segment about point B until it reaches point C.
6171. Babuška I & Miller A, 'The post-processing approach in the finite element method – Part 3: A posteriori error estimates and adaptive mesh selection', Int. J. Numer. Methods Eng. 20(12): 2311 - 2324 (1984)
6903. Baehmann PL, Wittchen SL, Shephard MS, Grice KR & Yerry MA, 'Robust, geometrically based, automatic two-dimensional mesh generation', Int. J. Num. Meth. Eng. 24(6): 1043 - 1078 (1987)
7056. Baker BS, Grosse E & Rafferty CS, 'Nonobtuse triangulation of polygons', Discr. & Comput. Geom. 3(2): 147 - 168 (1988)
7789. Bathe K-J, Wilson EL & Peterson FE, SAP IV. A structural analysis program for static and dynamic response of linear systems, Report No. EERC 73-11 (University of California, Berkeley), vii + 59 pp. + appendices, 1974
7790. Boubez TI, Funnell WRJ & Lowther DA, 'Three-dimensional finite-element mesh generation using serial sections', Proc. 11th Can. Med. & Biol. Eng. Conf.: 39 - 40 (1985)
7957. Boubez TI, Funnell WRJ, Lowther DA, Pinchuk AR & Silvester PP, 'Mesh generation for computational analysis. Part II – Geometric and topological considerations for three-dimensional mesh generation', Comp.-Aided Eng. J. 3(5): 196 - 201 (1986)
6796. Brown PR & Hayhurst DR, 'The use of the Schwarz-Christoffel transformation in mesh generation for the solution of two-dimensional problems', Computers in Engineering, Vol. 3, ASME, New York: 1 - 7 (1982)
6807. Buell WR & Bush BA, 'Mesh generation – a survey', J. Eng. Ind. (Trans. ASME) 95(1): 332 - 338 (1973)
6808. Cavendish JC, 'Automatic triangulation of arbitrary planar domains for the finite element method', Int. J. Num. Meth. Eng. 8(4): 679 - 696 (1974)
4840. Crane HL Jr, Gibbs NE, Poole WG Jr & Stockmeyer PK, 'Algorithm 508. Matrix bandwidth and profile reduction', ACM Trans. Math. Software 2(4): 375 - 377 (1976)
7791. Funnell WRJ, 'On the undamped natural frequencies and mode shapes of a finite-element model of the cat eardrum', J. Acoust. Soc. Am. 73(5): 1657 - 1661 (1983)
7792. Funnell WRJ, 'On the choice of a cost function for the reconstruction of surfaces by triangulation between contours', Comp. & Struct. 18(1): 23 - 26 (1984)
7870. Gordon WJ & Hall CA, 'Construction of curvilinear co-ordinate systems and applications to mesh generation', Int. J. Num. Meth. Eng. 7(4): 461 - 477 (1973)
4648. Hazony Y, 'Algorithms for parallel processing: curve and surface definition with Q-splines', Comp. & Graph. 4(3-4): 165 - 176 (1979)
7872. Heighway EA & Biddlecombe CS, 'Two dimensional automatic triangular mesh generation for the finite element electromagnetics package PE2D', IEEE Trans. MAG-18(2): 594 - 598 (1982)
7764. Ho-Le K, 'Finite element mesh generation methods: a review and classification', Comp. Aided Des. 20(1): 27 - 38 (1988)
7865. Jeng JH, Varadan VV & Varadan VK, 'Fractal finite element mesh generation for vibration problems', J. Acoust. Soc. Am. 82(5): 1829 - 1833 (1987)
7799. Jones RE, 'A self-organizing mesh generation program', J. Press. Vess. Technol. (Trans. ASME) 96(3): 193 - 199 (1974)
6809. Jones RE, 'The QMESH mesh generation package', ACM SIGNUM Newsletter 10(4): 31 - 34 (1975)
7867. Kalkani EC, 'Mesh generation program for highway excavation cuts', Int. J. Num. Meth. Eng. 8(2): 369 - 394 (1974)
6805. Kamel HA & Eisenstein HK, 'Automatic mesh generation in two and three dimensional inter-connected domains', High Speed Computing of Elastic Structures, BF de Veubeke (ed.), Univ. Liège, Liège: 455 - 475 (1971)
7834. Kela A, Perucchio R & Voelcker H, 'Toward automatic finite element analysis', Comp. Mech. Eng. 5(1): 57 - 71 (1986)
7151. Kikuchi N, 'Adaptive grid-design methods for finite element analysis', Comp. Meth. Appl. Mech. Eng. 55(1-2): 129 - 160 (1986)
7779. Lo SH, 'A new mesh generation scheme for arbitrary planar domains', Int. J. Num. Meth. Eng. 21: 1403 - 1426 (1985)
7965. Lo SH, 'Finite element mesh generation over curved surfaces', Comp. & Struct. 29(5): 731 - 742 (1988)
6811. Nguyen VP, 'Review of techniques for efficient network generation for finite element analysis', Forsch. Ing.-Wes. 48(2): 33 - 44 (1982)
7810. Shaw RD & Pitchen RG, 'Modifications to the Suhara-Fukuda method of network generation', Int. J. Num. Meth. Eng. 12(1): 93 - 99 (1978)
6800. Shephard MS & Yerry MA, 'An approach to automatic finite element mesh generation', Computers in Engineering, ASME, New York: 21 - 28 (1982)
6801. Sluiter MLC & Hansen DL, 'A general purpose automatic mesh generator for shell and solid finite elements', Computers in Engineering, ASME, New York: 29 - 34 (1982)
8331. Sornette D & Roux S, 'On the validity of "fractal" finite-element methods', J. Acoust. Soc. Am. 87(1): 411 - 413 (1990)
6806. Suhara J & Fukuda J, 'Automatic mesh generation for finite element analysis', Advances in Computational Methods in Structural Mechanics and Design, JT Oden, RW Clough & Y Yamamoto (eds.), UAH Press, Huntsville (Alabama): 607 - 624 (1972)
3992. Thacker WC, 'A brief review of techniques for generating irregular computational grids', Int. J. Num. Meth. Eng. 15(9): 1335 - 1341 (1980)
4292. Thacker WC, Gonzalez A & Putland GE, 'A method for automating the construction of irregular computational grids for storm surge forecast models', J. Comp. Phys. 37(3): 371 - 387 (1980)
4774. Tsybenko AS, Vashchenko NG, Krishchuk NG & Kulakovskii VN, 'Automatic formulation of triangular-element grid for arbitrary plane regions', Probl. Prochn. 12(12): 84 - 89 (1980) (in Russian)
Figure 1. Examples of generating evenly spaced boundary nodes. The small filled circles represent the initial points defining the shape, while the large open circles represent the automatically generated mesh boundary nodes. (a) Example of problems with sharp corners. Note especially the two large circles very close together at the sharp concavity. The two lowermost nodes are also too close together, and misrepresent the shape of the region. (b) Resolution of problems. In this case the outline was defined as separate right and left halves, thus forcing boundary nodes to be positioned exactly at the sharp tips and avoiding serious non-uniformities.
Figure 2. Example of generating regularly spaced internal nodes. The filled circles represent the same evenly spaced boundary nodes as shown in Figure 1b. The uniform triangular grid consists of 12 horizontal rows and 11 slanted columns. The squares represent grid points which have been accepted as internal nodes. The crosses represent grid points which are either outside the region altogether, or are inside the region but too close to the boundary to be acceptable.
Figure 3. Illustration of notation used in describing element generation (see text).
Figure 4. Internal mesh generated from the nodes shown in Figure 2. The solid lines represent elements formed entirely of newly generated internal nodes. The dashed lines define elements formed by joining external boundary nodes to nodes of the regular internal mesh.
Figure 5. Examples of meshes requiring special handling. (a) A case in which there are no acceptable internal nodes. The mesh is formed by optimally joining external boundary nodes to form triangles. (b) A case where an acceptable internal node (the small open circle at the lower right in the lefthand !BR!B) is not joined to the internal triangular mesh. As shown at the right, simply ignoring the node still results in an acceptable mesh after relaxation (see Section V). The resulting mesh would be even better if it were made finer. (c) An example of a mesh with two regions containing acceptable internal nodes, separated by a region with none. The thick lines represent elements generated with internal nodes, the thin lines represent elements involving only external nodes. The mesh shown has been relaxed.
Figure 6. The same mesh as shown in Figure 4, after relaxation. The shaded area represents the regular internal mesh which, before relaxation, consisted entirely of equilateral triangles.
Figure 7. A sample complex mesh involving two separate surfaces, and including holes. The filled circles are the labelled control points used in the definition of the mesh, as shown in Figure 8.
Figure 8. The mesh definition corresponding to the mesh of Figure 7. See text for explanation.
Figure 9. Perspective view of a mesh representing a three-dimensional curved shell. The shell is a model of the cat eardrum.
Figure 10. Meshes for the same geometry as shown in Figure 6, with nominal mesh resolutions of 6 elements/diameter and 30 elements/diameter. The mesh in Figure 6 has a nominal resolution of 12 elements/diameter.
AudiLab home page
