1 Abstract *
2 Introduction *
2.2 DXF Files *
3.2 Visualisation of Results *
3.3 Automatic Meshing *
4.2 Split Overlapping Objects *
4.3 Linearise Objects *
4.4 Form Closed Boundaries *
4.5 Check Boundary *
4.6 Adjust Order *
4.7 Hierarchy Determination *
4.8 Output Boundaries *
4.9 Viewer Module *
4.10 CommandLine Options *
4.11 Other Developments *
4.11.1 Debugging Information *
4.11.2 Boundary Viewer *
5.2 Adding Layer Sections *
5.3 Smoothing *
5.3.1 Boundary Node Smoothing *
5.3.2 Internal Node Smoothing *
5.3.3 Constrained Smoothing *
5.4 Seaming *
5.5 Joining Intersecting Boundaries *
5.6 Row Adjustment *
5.7 Closing Boundaries *
5.8 Cleanup Mesh *
5.9 Writing Meshes *
5.10 Command Line Options *
5.11 Other Developments *
5.11.1 Debugging Information *
5.11.2 Connection Viewer *
5.11.3 Mesh Viewer *
5.11.4 Settings File *
6.2 Asymmetry *
6.3 Sensitivity to Segment Lengths *
6.4 Performance *
8 References *
9 Glossary *
This report describes algorithms for converting 2D DXF files into a hierarchy of linear segment object descriptions and producing from these an unstructured all quadrilateral mesh, suitable for finite element analysis applications.
Some simple examples of the meshes automatically generated are presented and areas where the methods used could perform undesirably are outlined, with suggestions on how to avoid or mitigate these circumstances.
The brief given for the project was:
"AutoCad is widely used to produce engineering drawings of objects. The first aim of this project is to develop a program that can interpret geometrical objects described in AutoCad format and reconstruct them on screen. The second aim is to develop an algorithm which can break down the objects into simpler shaped regions automatically (quadrilaterals and triangles) so that they can be used as the input to a finite elements package at a later stage."
The essence of this is a requirement to automatically generate meshes from a design, minimising the need for human input. This is an area of considerable academic and commercial research. The final target is a completely automatic system that will allow a user to design a product and model it’s characteristics, without requiring any knowledge of the modelling processes.
Some brief background information on relevant topics, namely finite element analysis and DXF files are useful before proceeding further.
The method of finite elements is used extensively for modelling the behaviour of designs before manufacture or testing. It essentially works by dividing the model up into an arbitrary number of discrete domains, the finite elements. A numerical model is then used to approximate the behaviour of each element to meet a set of boundary conditions. If a good enough model and sufficient elements are used the simulated system approaches the real answer.
The elements are either triangular or quadrilateral in 2D or tetrahedral and brick shaped in 3D. The triangular elements provide faster and simpler modelling but quad elements produce more accurate answers. When modelling a system the division of the modelled object into elements is very important in facilitating accurate results and rapid calculations. The points in a model and their connectivity are often called a mesh, the connections being the threads between the knots (nodes) of an imaginary section of netting.
Meshing can be done by hand and an experienced eye can produce good meshes by taking many factors into account. Depending on the complexity of the model and the accuracy desired this could take a long time, up to manyears. To reduce this expense several automatic techniques have been proposed to allow computers to undertake most of the necessary work and cut down the time required to prepare a model.
DXF (Drawing eXchange Format) files were originally introduced in the AutoCAD series of computer aided design tools. They are now widely supported by products in the CAD/graphics sector.
DXF files are able to represent many different objects used in technical drawing and allow objects to have a wide range of properties such as layer, colour, line width, line style etc. DXF files are available in two forms, ASCII text that can be read and edited using many text editors and a more compact binary form. The ASCII form is more widely supported and easier to read so only this form will be considered. The objects that can be represented range from simple straight lines to complex nonuniform rational Bspline (NURBS), text and dimension objects.
The DXF file format stores information in pairs of values, first an integer to describe the data type of the next value, followed by the data value. Data types supported include integer, floating point, text strings and hexadecimal values.
The DXF file is divided up into sections delimited by control strings. Each section has a set range of content but need not contain any information. The range of content has increased over the years with successive additions to the file format, but in a consistent manner which allows good backward compatibility.
For simplicity it was decided to create two separate programs, one to allow viewing and linearisation of the DXF files and another to perform the meshing. Discussions with the project supervisor established that the programs would be written in ANSI compliant C, and should be as platform independent as possible. The source code for all algorithms and tools used in this project are included in Appendices B to F.
After a brief perusal of the DXF file format [1] it became clear that supporting a complete set of objects would be complicated and unnecessary. To reduce the number of structures to support the following assumptions were made:
Only 2D objects will be read (even if a 3D meshing algorithm is used it will be easily adapted to 2D).
Only objects that can be considered as lines will be read (text is not needed and shapes are likely to be formed from line primitives).
This leaves a relatively small set of objects, namely lines and curves, to interpret. Of these the most complex is the SPLINE which is a version of a NURBS curve. A little research in the standard references material [2], illustrated that this would be complicated to implement and there is little need as AutoCAD, and most other CAD packages, support converting such curves into linear segments. The objects to interpret are now narrowed down to LINE, POLYLINE, ARC, CIRCLE and ELLIPSE.
To form a mesh it is necessary to have linear edges to the shape being meshed (triangles and quadrilaterals do not have curved sides). This will require approximating curves with linear segments.
The next consideration is how to link up the objects that form the outline of the shape to model. It is reasonable to assume that the endpoints of objects forming a shape will match exactly. However if more than two endpoints match more information is required. The simplest way to achieve this is to use the layer property of the objects read in, only objects with matching endpoints and the same layer will form a shape. As the layer property is a string it can also be used to convey information about the shape, such as it's material properties.
Having produced a shape it must be written out in a consistent manner, which will arbitrarily set as anticlockwise ordering of points. To express the hierarchy of other shapes inside the current shape only the shapes directly within the current shape will be written out clockwise, to form "holes", any shape inside a shape inside the current shape will not be included in the current shape's output. By doing this for every shape a full description of all shapes in a model are expressed in terms of an exterior and, possibly some, interior boundaries. An example of this is given later and should clarify the concept.
To allow viewing of the DXF objects read in, a graphical user interface (GUI) of some kind is required. As graphic devices and interfaces are system dependent there is no way to ensure that graphics can be shown on any platform. However there are some standardised graphics packages that are available for a range of platforms. These include SPHIGS, OpenGL and SRGP. Of these OpenGL is the most widely supported and likely to be so in the future, and free implementations (both proprietary and open source [3]) are available, so this is the chosen format for display. By itself OpenGL is capable of rendering an image to a buffer but does not support display of the image generated. This can be accomplished using the OpenGL Utility Toolkit (GLUT) which supports a simple interface for displaying windowed graphics and handling user interaction on Win32 and XWindows platforms.
A wide variety of automatic meshing algorithms have been described in professional literature, a good summary of some techniques being given in Reference [4]. Many of the algorithms are combined and implemented in commercial and research [5] meshing packages.
Automatic mesh generation techniques can be split roughly into two categories, structured and unstructured meshes. In structured meshes every interior point must be connected to a fixed number, sometimes called the valency, of other nodes, in the case of 2D quadrilaterals this is four. In unstructured meshes there is no restriction, and are easier to create, as valency constraints require complex iterative smoothing. Structured meshes are mainly used in computed fluid dynamics. As the particular interest of this project is to model electromagnetic systems unstructured meshes are acceptable.
To allow maximum utility of the resulting program it is best to implement a mesher producing quadrilateral elements, as these can be easily split into triangles for FE packages that require simpler elements.
The meshing algorithm implemented will be in 2D as this is much less complex than 3D and can be useful in some 3D situations, namely in electrical machines which are axially homogenous.
There are a variety of approaches to 2D unstructured quad mesh generation [6], including:
Moving front algorithms work by projecting elements from the interior and exterior boundaries and using these to form the next active boundary to add elements to. In this way layers of elements are added which tend to have the best elements at the surfaces of the object. This is a desirable characteristic when using finite elements to model electromagnetic problems, as the surfaces of objects are generally the areas of interest.
The input file to the boundary creator will be a 2D DXF format drawing (XY plane only, Z components will be ignored), and the output will be a boundary description file in a simple text format, described later. A boundary in the input file is formed from a closed loop of entities with the same layer string. Boundaries may not overlap but may be identical in sections i.e. two boundaries may be exactly adjacent with a common edge. Boundaries must also be simple i.e. no selfintersections. These characteristics are summarised in Fig. 1.
A brief overview of the processes carried out is as follows:
4.1 Reading
Data from DXF files
To read DXF files usefully it is necessary to translate string information into a format in which decisions can be made, and to read in numeric values correctly. If the line contains a string a comparison is made with a set of predefined control strings i.e. all the strings the program is expected to understand, and if a match found the corresponding flag returned, if no match is found a flag indicating an unknown string is returned. Another function can then be used to retrieve the string for further analysis.
A description of all supported entities and their component values are given in Appendix A.
The DXF file is subdivided into several sections but the only ones used by the boundary converter are the HEADER and ENITIES sections.
The HEADERS section describes general attributes of the file and of the way in which it is to be displayed. For the purposes of this project the only relevant information is on the angle description (direction and degrees/radians/grads format used) and precision of angles, as these are needed to read ARC entities correctly.
The ENTITIES section of the file is contains all the geometrical objects of interest, in this case ARC, CIRCLE, LINE, POLYLINE and ELLIPSE entities. Each entity can have a variety of associated properties but the only one used here is the layer string, which gives the name of layer the entity belongs to.
When loading a DXF file the relevant angle information is extracted from the HEADERS section. The ENITITIES is searched for any coordinate entries and the number of decimal places used in their representation is noted, used later for correct truncation of points.
The ENTITIES section is scanned for recognised objects, which are read into memory and any necessary angle conversions performed. The layer string of each entity is also read in and stored in a list for later use.
The object list is then searched for objects with the same type and
overlaying extrema. If an overlay is found the objects are split so that
the end points of an object are only coincident with the end of another
object. This is necessary to avoid cases where the linearised approximations
of geometrically similar objects are not identical, as illustrated in Fig.
2. If such a case were to occur the result would be overlapping meshes.
Figure 2 Splitting Objects
The objects are then converted to linearised approximations of their real shape. This is initially performed on the basis of longest length of linear segment to be produced. This can be specified by the user in the command line on invoking the DXF to boundary converter or a value calculated internally. This is translated into the number of segments to approximate the objects with. The parametric representation of the object is then used to find points on the object to form the endpoints of linear segments. The linear approximations are stored as a list of lists of points, an N segment approximation needing N+1 points.
The end points of each linear approximation are truncated to the same
precision as in the original DXF file. This is done to allow the correct
linking of objects by linear approximation, Fig. 3b, not a problem for
LINE entities but very necessary for ARC and ELLIPSE entities. If endpoints
do not match no boundary can be formed, as illustrated in Fig. 3a.
Figure 3 Truncation of EndPoints
It is worth noting that the precision of angular units can effect the
truncation, if too coarse an angular scale is used a point on a curve may
be truncated wrongly. Also note that the area of approximated regions is
not the same as the original where curved edges are used. It would be possible
to produce equal area circular arc approximations by altering the radius
of the circle approximated. Elliptical arcs would be much more difficult
to handle so no area equalisation was implemented.
Firstly a check is made that there are not more than two instances of an end point in the same layer, if there are then a correct boundary may not be formed.
A boundary is formed by searching for an approximation in the same layer that has an end point coincident with the end point of the current approximation. The matching approximation and the end of it that is coincident with the current approximation recorded. The matching approximation becomes the current approximation and this is iterated until the original approximation is reached. The boundary can now be read out by going through the approximations in the order matched and copying the points in the order of matching (forwards or backwards).
This is procedure is repeated until all approximations have been accounted for. A search for duplicate boundaries is made, as duplicates would prevent the correct calculation of boundary hierarchy.
Each boundary is checked for selfintersections by checking that no segment intersects any other segment, using a parametric method, derived from the CyrusBeck lineclipping algorithm [7].
Two intersecting vectors are shown in Fig. 4 and illustrate the method.
Where two lines (P_{1} & P_{2}) intersect the dot product
of the normal to one line (N_{2}) and the vector formed by subtracting
a point on the line (P_{22}) from it’s intersector (P_{1})
is zero, see Equation 1. Substituting the parametric form of the lines,
Equations 2 & 3, give the equations of the parameters as Equations
4 & 5.
Figure 4 Parametric Intersection
If 1<u<0 and 1<v<0 the intersection is within both segments. If the segments are parallel the lower term in Equations 4 & 5 will be zero which would cause an error so this case is checked for and rejected.
If no selfintersections are found the boundary is simple, a nonsimple boundary could escape the next stage erroneously.
To perform the later stages it is necessary to ensure that all boundaries are described in the same order.
The external angle of a boundary can be found by summing the product of the angle of each vertex, found using Equation 6, by the sign of the signed area, Equation 7. The sign of the signed area is equivalent to checking whether the node following a vertex is to the right (negative) or left (positive) of the vector to the vertex from the node prior to it. This is briefly shown in Fig. 5.
If the sum is negative the boundary was scanned in the wrong order i.e.
clockwise, so the boundary is reversed to give an anticlockwise ordering,
necessary for working out the hierarchy of the boundaries.
Figure 5 Order Calculation
a, b & c are the vertices of a triangle associated with angles
A, B & C.
The hierarchy of the boundaries is defined as which boundaries are inside which. An example is Fig.6 where B is inside A, and C is inside both A and B.
The hierarchy of the boundaries is then found by taking pairs of boundaries
and checking whether one is completely inside the other. By keeping track
of which boundaries are inside which it is possible to avoid checking every
pair of boundaries (N² complexity). First of all a check is made
for whether the boundaries overlap by searching for intersections. If they
do overlap they are invalid and the program exits. However if no intersections
are found but sections of identical boundary detected the special case
of an insert has occurred. This is illustrated in Fig. 7a. and requires
splitting the boundaries so that no boundary inside another has identical
boundary sections, as this would be impossible to mesh, there being no
area between the boundaries to fit elements in.
Figure 6 Hierarchy
Figure 7 Insert Handling
For checking whether a boundary is inside another a raycasting technique is used, illustrated in Fig. 8. From each segment an outward normal vector is projected and intersections of this ray with segments of both boundaries are checked for using the parametric method described above. The boundary intersected and intersect distance are stored, identical distances and boundary numbers removed, then the list is sorted by intersect distance.
If an even number of intersections have occurred and there is only one intersection on the origin of the ray it is assumed that one boundary is inside the other and a consistency check performed. This works by checking that the intersections with the assumed inner boundary only occur inside the outer boundary. If this check is passed and splitting is necessary the boundaries are searched for the start and end points of identical edge sections and the points between these copied to new boundaries, as in Fig 7b.
Once this splitting has been performed the hierarchy determination process is restarted so that the new boundaries are taken into account.
Once the hierarchy is known it is necessary to eliminate redundant information,
any boundaries inside a boundary inside the current boundary must be removed
from the current boundary’s list of inside boundaries (they are masked
by the boundary immediately enclosing them). This is demonstrated in Fig.
9. The three boundaries A, B & C will form three distinct objects to
mesh, with the following exterior and interior boundaries:
Figure 6 Boundary Hierarchy
Object to Mesh
Exterior Boundary
Interior Boundary(s)
1
A
B
2
B
C
3
C

The hierarchy is then used to write out the boundaries in an ASCII text file formatted as follows:
"BOUNDARY" nHoles nPoints nHoles & nPoints are integers
"INFO" STRING STRING is the LAYER string of the exterior boundary
X Y nPoints sets of floating point coordinates, ordered anticlockwise around boundary
Followed by nHoles entries as follows:
"HOLE" nPoints nPoints
X Y nPoints points ordered clockwise
The output file may contain more than one boundary and will be named as the input file except with a BND file extension. Note that if any holes are present each hole will have it's own entry as a BOUNDARY.
The viewer module may be started if the relevant preprocessor defines and libraries were included at compile time and the user specifies it in the command line. The viewer module makes use of the OpenGL API, GL Utilities (GLU) and the GL Utility Toolkit (GLUT). The first two are available with most operating systems and GLUT is available for Win32 and XWindows systems, making the viewer available to the majority of desktop systems i.e. Windows and Linux.
The viewer displays the linear approximations and for curved objects a dashed line approximation to the curve (it is possible for the linear approximation to be a better approximation to the curve than the dashed version) as in Fig 10.
Figure 10 Viewer/Editor Screenshot
The viewing and interaction operations include panning, zooming, selection, addition and removal of segments and points and the moving of individual points.
Pan and Zoom
The user may pan using the cursor keys and zoom in and out using the ‘I’ an ‘O’ keys or the appropriate options from the right click menu.
Selecting an Approximation
Clicking on a linear approximation will select it and highlight it, showing the points that the user may adjust, see Fig. 10. The layer information for that object and any exactly coincident object will be displayed in the title bar.
Adding/Removing Segments
The user can increase or decrease the number of segments in the selected linear approximation by right clicking and selecting "Add/Remove Points" the appropriate number of points from the submenu. The linear approximation will then be updated with the points defining it's segments spread evenly (linearly for line and angularly for curved objects).
Bisecting Segments
By clicking on a segment in the selected approximation and selecting "Insert Point on line" from the right click menu the user may add a point the centre of a segment to decrease segment size locally.
Moving Points
The user may also select a point on the highlighted approximation and move it by clicking and holding down the left mouse button on the point and dragging the mouse. The new point and adjoining segments will be displayed in dashed form while the user drags the point. During point movement the position of the point will be displayed in the title bar. Releasing the left mouse button will place the point.
The movement of the last point can be undone by right click menu selecting "Cancel Move".
When selecting and moving a point holding down the shift key will constrain the point's motion to be exactly on the object and to be between the points on either side. This should make it easier to place points exactly on the boundary.
Deleting Points
Individual points can be removed by selecting a point and then selecting "Delete Point" from the right click menu.
Where approximations are exactly coincident any addition, removal or movement of points is performed on both objects.
When the user is satisfied with the linearisation of boundaries right click menu "Generate Boundary" will generate the boundary, as described above.
When compiling the DXF to boundary conversion tool all the files in Appendix B must be included.
To make use of the viewer module VIEWGL must be defined (commonly a DVIEWGL command line argument to the compiler) and the OpenGL, GLU and GLUT libraries included during the linking stage.
When invoking the tool one, two or three command line arguments are used. The first argument is mandatory and is the path and file name of the input DXF file. Note that the output file will have the same name but with a BND extension. The optional second argument is the floatingpoint maximum length of segment to allow in the linear approximations, if not included or 0.0 is specified a value is estimated. The final argument is optional, if "V" is specified and the viewer module was included at compile time the viewer is invoked after loading the DXF file and the user may edit the linear approximations.
During development it became necessary to know the results of the major stages of processing and to this end an additional set of functions were implemented (see Appendix A debug.c). If this is included in the compilation with the preprocessor define DEBUG the files objs.txt and las.txt will be created in the working directory listing respectively entities after splitting and the endpoints of the linear approximations before linking into boundaries.
4.11.2 Boundary Viewer
To permit visual checking of the resulting boundaries a boundary viewer
has been created (see Appendix C for source code). This displays differing
layers in different colours to help in interpreting the boundaries produced,
an example output is given in Fig. 11.
Figure 11 Boundary Viewer
The meshing algorithm used is a variation of the "Paving" approach first proposed by Blacker & Stephenson [8]. This algorithm produces unstructured meshes by layering quadrilaterals from the edges of the object.
The principal variations from the algorithm originally described are:
Each operation is described below
Figure 12 Paving Flow Diagram
Objects to meshed are loaded from an input file, using the BND file format, one at a time. Each boundary (external and internal) is checked to ensure that it has an even number of points. If an odd number is found an additional point inserted in the middle of the longest boundary segment.
Vertices on a boundary can be classified by the angle (q ) at the node N_{i} as row end, row side, row corner or row reversal nodes roughly as follows:
Row End 0<=q <=Pi/2
Row Side Pi/2<=q <=Pi
Row Corner Pi<=q <=3Pi/2
Row Reversal 3Pi/2<=q <=2Pi
Quadrilaterals are added to the mesh by projecting nodes from a row end node on the current boundary and using these to form new elements. Row end nodes can usually be found somewhere on the boundary to start from but if none are found one can be created by projecting two vectors to form a new quad, as in Fig. 13. Once a row has been started the projection necessary is selected primarily on the basis of the angle at the boundary node. These projections are row side, row corner and row reversal projections. These are illustrated in Figures 14 to 16 respectively and are all formed by rotating the vector from the boundary point to it's predecessor by a fraction of the vertex angle, and adjusting the magnitude.
If another row end node is encountered the layer section is completed by forming a quad as in Fig. 17.
Figure 13 Row End Creation
Figure 14 Row Side Projection
Figure 15 Row Corner Projection
Figure 16 Row Reversal Projection
Figure 17 Row End
The original algorithm proposes selecting the projections used by classifying a portion to add as a primitive and finding which projections produce the least number of irregular nodes. This is potentially quite computationally costly so a simplified method has been implemented. This is based on the angle of the node and the number of other nodes connected to that node, hopefully minimising irregular (i.e. valency not equal to four) nodes.
Two different forms of smoothing are recognised, boundary and internal node smoothing. Boundary node smoothing is based on an isoparametric smooth with adjustments for desired length and angular smoothness. Internal node smoothing is a length weighted Laplacian smooth with angular smoothness adjustment near boundaries.
Boundary smoothing is performed on nodes on a boundary that are not fixed i.e. not in the original object description. After a boundary smooth a smooth of internal nodes connected up to a depth of N is carried out. N equal to 3 is considered sufficient and allows any major changes at the boundary to be spread out through the mesh.
Internal smoothing is used on all nonfixed nodes that are not part of a boundary.
The isoparametric smooth is found by taking a quad that the point being
smoothed (N_{i}) is in. The vectors from the origin to the vertices
either side of the point (N_{j} & N_{l}) are summed
and then the vector from the origin to the vertex opposite (N_{k})
and the vector to the point being smoothed subtracted. This is summed for
all M quads the point is in to give the movement vector (V_{a})
from N_{i} to the isoparametric point. This is summarised in Fig.
18.
Figure 18 Isoparametric Smooth
The length adjustment is made by scaling the vector from the isoparametric
smoothed point to the node from which the node being smoothed was projected
(N_{i}N_{j}V_{a}), to the length from N_{i}
to N_{j} (referred to as the desired distance l_{D}), to
give a new vector V_{b}, as in Fig. 19. This preserves the size
and aspect ratios of elements.
Figure 19 Length Adjustment
Following this an adjustment is made for the angular smoothness. The
vector bisecting the angle between N_{j} and the opposing vertices
of the two adjoining quads is found (V_{bis1}), as illustrated
in Fig. 20. This is then used to find the vector bisecting it and the vector
from N_{j} to N_{i} (V_{bis2}). The length of this
vector (l_{Q}) is the distance from N_{j} to it's intersection
with a line joining the opposing vertices in the adjoining quads (Q). The
length is adjusted by comparing l_{Q} and l_{D}, if l_{D}
is longer than l_{Q} the magnitude is adjusted to the average of
l_{D} and l_{Q}, otherwise it is the same as l_{D}.
Figure 20 Angular Smoothness
This gives an angular adjustment vector of:
V_{c} = V_{bis2}  V_{i}
These are then averaged to give a boundary smoothing vector of:
V_{BND} = (V_{b} + V_{c})/2
Internal node smoothing is simpler and is demonstrated in Fig. 21. For
each of the N nodes connected to the node being smoothed (N_{1}N_{4})
a vector is formed of the vector from the node being smoothed (Ni) to the
node in question (V_{1}V_{4}). The sum of these vectors
would be the Laplacian, which attempts to move a node to the centre of
all its connected neighbours. If the node in question (N_{?})is
on boundary (permanent or temporary) the angular adjustment vector described
above (V_{a}) is added to the contributing vector, to give a modified
contribution V’_{2}. The length weighted sum over all connected
nodes is the adjustment vector, V_{L} as in Equation 8.
Figure 21 Laplacian Smooth
The above methods are relatively simple and efficient to implement, converging in 3 or 4 iterations. However these smoothing methods can reduce the quality of quadrilaterals so a constrained smooth can be enabled. This only allows nodes to be moved if the movement produces an improvement in the quadrilaterals attached to the node. To achieve this the concept of a quality metric for triangles and quadrilaterals is introduced [9].
The metric used breaks down the quadrilateral into the four possible
triangles, by using the two diagonals of the quad as shown in Fig. 22,
and evaluates the "equilateralness" of each of those triangles, using Equation
9. The quality metric of each triangle is in the range 1 to 1, lower values
indicating a more distorted element, negative values indicating inverted
elements.
Figure 22 Component Triangles
A, B & C are the vertices of an anticlockwise triangle.
In a quad the metric is the minimum metric of the four triangles in that can be made from the quad with an adjustment for two or more inverted triangles, given in Equation 10.
N is the number of negative q s minus one.
When constraining a smooth all quads attached to a node are examined and the movement is allowed only if the metric improves for all quads, improves significantly for some quads without significantly decreasing it for any others or improving the average metric. If this is not the case the movement vector is bisected and the metrics reevaluated, an arbitrary number of times.
This is the operation by which small angle and inverted node vertices are removed by shifting the points on the boundary either side of the vertex to a common point, as illustrated in Fig. 23. The threshold angle at which this occurs is larger for irregular nodes than regular nodes, to encourage the propagation of irregular nodes away from the mesh surface.
Figure 23 Seaming
When large differences in segment length on either side of the vertex are found the longer segment is split into three and a quad inserted, as illustrated in Fig. 24, to help the transition of element sizes. This can be applied iterativley to allow large transitions in element size ([8] states that length differences of up to 30 can be accommodated).
Figure 24 Transition Seam
5.5 Joining Intersecting Boundaries
When boundaries overlap it is necessary to connect them. The overlap is detected by using the same parametric method described earlier for the DXF to boundary converter. Once an overlap is discovered all overlaps between the boundaries involved are listed. The pairs of segments intersecting and the pairs of segments either side are evaluated to find the best way to join the boundaries. The evaluation checks that the intersecting boundaries will leave an even number of segments on either side of the join. It also attempts to reduce any disparities in intersecting segment lengths by allowing for the insertion of extra quads. The goodness of the join is ascertained by taking into account the separation of the points to join and the angle between the segments, as in Equation 11. This metric can be any positive value from 0 up, the lower the metric the better the join. The threshold (a ) prevents nearly parallel segments from producing excessively low (i.e. very good) metrics even when the separation of points is large.
P_{11}P_{12} and P_{21}P_{22} are
the endpoints of the intersecting segments, none of which are coincident.
f is the angle between the segments.
a is the minimum allowable sin(f
)
Once the best place to join has been found any necessary transition quads are inserted by adding points to a segment 1/3 and 2/3 along and projecting two points from these, parallel to the sides of the quad. The original quad is removed and replaced by four quads, as in Fig. 25. This is iterated and allows segments of very different lengths to be joined.
Figure 25 Transition Quad
Now the boundaries are joined by merging the end points of the segments, shown in Fig. 26. If a boundary is intersecting itself two boundaries will be formed, if two different boundaries intersect once joined there will be only one resulting boundary. If a boundary is interior and self intersects one of the resulting boundaries will be interior so a check is made by calculating the total external angle of each boundary, if this is negative the boundary is marked as interior. If two boundaries intersect and both are interior the resulting boundary will be interior and is marked as such.
Figure 26 Joining Boundaries
Once a layer has been added to a mesh it may be necessary to adjust the mesh in areas of change to maintain as regular a mesh as possible. These are primarily the insertion of quads in expanding section (wedging) and the deletion of quads from contracting sections (tucking). Areas of the layer boundary where these are needed are detected by looking for runs of vertices with concave or convex angles. Once such a run has been detected elements are added/removed as close to the centre of each Pi/2 section of external angle as possible.
The wedging operation is only performed if the angle at the node is greater than a threshold (must be greater than Pi) and the adjoining boundary quads expand more than a ratio, i.e. d_{2}C×d_{1} and d_{4}C×d_{3.} The wedge is inserted by adding two points on existing segments and projecting a third to produce a quad to form a corner, demonstrated in Fig. 27.
Figure 27 Wegding
Tucking may only occur where two succeeding vertex angles are less than a threshold (at most Pi) and where the quad between them is contracting more than a ratio i.e. d_{2}<C×d_{1}. A tuck is formed by removing the two points and their connecting quad from the boundary and adjusting the adjoining quads, Q1 and Q2, to fill in the gap in the boundary, as shown in Fig. 28.
Figure 28 Tucking
When a boundary has few enough vertices some meshing operations can produce undesirable results e.g. seaming a two vertex boundary does nothing. To prevent such cases it is necessary to remove small boundaries and complete meshing using special cases, rather than the general algorithm. After any operation that reduces the number of vertices in a boundary a check is made whether the boundary operated on has six or less nodes and is not an interior boundary.
If a boundary has zero or two nodes no action is necessary, the mesh is complete. If there are four nodes a quad formed from those nodes is added. In the case of six nodes the original algorithm [8] proposes matching the boundary to one of a set of primitives to insert 2, 3 or 4 quads. However as closing a boundary is likely to occur far from the edges of a meshed object the quads inserted do not have to be well formed, as they have little effect on the simulation result at the surfaces. A simpler method of closing a six vertex boundary is used in which all possible bisectors that produce two quads are searched for that which will produce the maximum (best) average distortion metric of the quads produced. This is illustrated in Fig. 29.
Figure 29 Closure
Automatically generated meshes often contain elements that have undesirable characteristics and a wide variety of improvements are possible, some being described in Ref. [10]. The only mesh improvements implemented are node improvement by quad removal and quad improvement by quad merging.
In the first case quads are removed if the regularity of the resulting nodes improves upon that of the original nodes. This is illustrated in Fig. 30 and works by merging two opposing vertices of a quad.
Figure 30 Quad Removal
The second case searches for interior nodes that are connected to only two quads and also the quads share two edges, then combines the quads into one, as in Fig. 31.
Figure 31 Combining Quads
In some case this may lead to the formation of an undesirable chevron shaped quad, Fig. 32b, in which case combination is not performed. Instead the node connected only to those quads is moved so that it is on the line joining the nodes either side of it in the quad. This will produce "triangular" quads (quad metric can be at most 0), Fig 32c, but will avoid the case where one quad (Q1) has a bad metric (less than 0), Fig 32a.
Figure 32 Chevron Quad
The process of cleanup not only improves the mesh but also reduces the number of nodes, reducing the processing required in finite element analysis.
Once an object has been paved the resultant mesh is added to a global mesh, only new points are added to this mesh, as identical points will occur at the interface between meshes, and the quadrilaterals copied. The information string is also added to a list, if not already in use, and each quad in the added mesh elements given an index to the string. This allows sets of boundaries describing complete systems to be formed into one mesh.
When all objects from the file has been meshed the global mesh is written out as an ASCII text file using the following format:
"MESH" nInfos nPoints nQuads three integers.
...
String
String of property description. nInfos entries.
...
X Y Floating point coordinates. nPoints entries.
...
Q1 Q2 Q3 Q4 Qstring Q1Q4 Zero based indices of points forming quad, in anticlockwise order. Qstring is a zerobased index of property string, if no property string is applicable this should be –1. nQuads entries must occur.
Only one argument needs to be supplied when invoking the paving program, the path and name of the *.BND file containing the objects to mesh. The mesh file will be created in the same directory with the same name but a .MSH extension.
To assist in debugging the algorithm a set of functions, listed in debug.c in Appendix D, for outputting the state of the meshing task were created. A full description of the mesh quads, points and active boundaries can be written to a text file for in depth analysis. The current mesh of quads is also written out in the same format as the global mesh output.
For a rapid judgement of the state of the meshing in progress the connections between points may be extracted from the text file and displayed using a graphical viewer, shown in Figure 33 and listed in Appendix E.
Figure 33 Connection Viewer
To view the mesh format another viewer (Appendix F) has been implemented. This draws the quad perimeters and fills them with a colour based on their property string. An example is given in Fig. 34.
Figure 34 Mesh Viewer
1.25 Wedge ratio
0.0628 Wedge angle threshold (PI + val) (approx 3.6 degrees)
0.8 Tuck ratio
0.0628 Tuck angle (PI  val) (approx 3.6 degrees)
0.7853981634 Seam angle for nodes with valency of 4 (approx. 45 degrees)
0.9299310886 Seam angle for nodes with valency not 4 (approx. 60 degrees)
2.5 Seam wedge ratio
0.25 Minimum allowable sin(q ) used in joining metric
3 Boundary smooth depth
4 Smoothing iterations (should converge in 3 or 4)
1 Constrain interior smooth (1 = on, 0 = off)
1 Constrain boundary smooth (1 = on, 0 = off)
10 Constraining iterations (should be less than number of bits in floating point number)
Several test cases were attempted, starting with DXF files and resulting in meshes. Only very simple designs were tested and the DXF to boundary converter behaved as expected. However the meshing of the designs highlighted some potential issues with the paving algorithm.
6.1 Breakout
"Breakout" is a well known problem with the paving algorithm and is mentioned in the literature [11]. This occurs when boundaries do not intersect correctly and overlapping meshes are produced. This is particularly likely when boundaries with short linear segments are close to boundaries with long segments as the former boundary may be entirely enclosed by a single layer projected from the latter, as shown in Fig. 35.
Figure 35 Breakout
It also occurs when joining of the intersecting boundaries is not handled correctly, as illustrated by the shaded are in Fig 36. In this instance the seaming operation after the joining of the boundaries has failed, owing to the "butterfly" quad Q1. This led to fatal errors later on and no valid mesh was produced.
Figure 36 Breakout  Failure of Seaming
One way to avoid these problems is to alter the algorithm is to check for intersections on an element by element basis [11]. However this can produce nonoptimal boundary joining and requires considerable reworking of the algorithm structure.
An alternative approach would be to detect overlaps after the addition of a layer portion and remove any overlapping quads produced. This would fit more easily into the flow of paving but would require methods for removing quads from the surface and restoring the previous boundary. This could be complex (if doing so "on the fly") or memory consuming (if keeping a copy of the mesh state prior to addition of the layer portion).
The simplest solution is to ensure that there are no segments longer than the smallest dimension of the object to mesh. This has the disadvantage of slowing down the paving algorithm, as more elements must be produced and also slowing down any finite element analysis by increasing the number of nodes in the problem.
6.2 Asymmetry
Paving a model that exhibits symmetry can result in an asymmetric mesh, as shown in Fig. 37. This is not desirable as symmetric meshes can be reduced into constituent meshes and a lower amount of modelling required, particularly true in electrical machines. To take advantage of this it is necessary to take any symmetry into account before paving the model.
Figure 37 Asymmetric Mesh
6.3 Sensitivity to Segment Lengths
The paving algorithm can cope with large variations in segment length as in Fig 38. The resulting mesh is not optimal (it has an average quad metric of 0.436286) but does have a small number of nodes (74 points forming 56 quads).
Figure 38 Different Length Segments
The same boundary described with more nodes, giving more even segment lengths, produces a much improved mesh, with an average quad metric of 0.600143, but with more nodes (176 points forming 140 quads) Fig 39.
Figure 39 Similar Length Segments
If the amount of computation required for any following finite element analysis is not critical it is preferable to produce a finer and better mesh by using short boundary segments. A good approach would be to start with a very short segment length and gradually increase it until just before the algorithm failed, or was just producing a mesh of the desired quality.
6.4 Performance
To mesh the cases shown in Figures 32 and 33 took 0.985 and 2.679 seconds respectively on a 166MHz Pentium class PC, compiled with lccwin32 [12]. This suggests that a realistically large model, thousands of segments in the boundaries, could be meshed in acceptable time (minutes) on a modern system. As the paving algorithm is sequential in nature the use of parallel computing is not likely to improve the speed of the algorithm.
One possible consideration is the memory usage, the limiting factor is not likely to be the size of available memory on current systems but the allocation of relatively large monolithic blocks of memory may be of concern. If this were a problem the implementation could be adjusted to make use of smaller memory blocks.7 Conclusions
The aims of the project brief have been met; DXF file objects can be displayed and decomposed into quadrilateral elements.
The algorithms used are relatively simple and areas where problems may occur are understood and readily identifiable. From the limited examples presented the paving algorithm appears quite robust and capable of producing good meshes, provided some consideration is made of limitations of the method.
The algorithms have been implemented in ANSI C and do not show excessive computational cost. Given suitable minimal user interaction a CAD design can be formed into an acceptable mesh.
It is unlikely that the method presented, or any other automated mesh generation scheme, will totally replace high quality hand meshing. However such methods offer a much more rapid design cycle and as computation becomes cheaper the need for good quality meshes with minimal nodes will become less pressing.8 References
[1] AutoCAD Release 14 DXF Reference at http://www.autodesk.com/techpubs/autocad/acadr14/dxf/index.htm
[2] NURBS Book (Second Edition), L. Piegl & W. Tiller, SpringerVerlag (1997), ISBN 3540615458.
[3] MESA an OpenGL conformant graphics library, source code available from http://www.mesa3d.org/
[4] A Survey of Unstructured Mesh Generation Technology, S. Owen, available at http://www.andrew.cmu.edu/user/sowen/survey/index.html.
[5] CUBIT Toolkit for robust and unattended mesh generation http://endo.sandia.gov/cubit/.
[6] Automatic mesh generation: application to finite element methods, P.L. George, Wiley (1991), ISBN 0471930970
[7] Introduction Computer Graphics, J. Foley et al, AddisonWesley (1990), ISBN 0201609215
[8] Paving: A New Approach to Automated Quadrilateral Mesh Generation, T. Blacker & M. Stephenson, International Journal for Numerical Methods in Engineering, Vol 32 811847 (1991)
[9] An Approach to Combined Laplacian and Optimisation Based Smoothing for Triangular, Quadrilateral and Quad Dominant Meshes, S. Canaan, J. Tristano & M. Staten, Presented at the 7^{th} International Meshing Roundtable (1998). Available from http://www.imr.sandia.gov/papers/imr7/canann98.ps.gz
[10] CleanUp: Improving Quadrilateral Finite Element Meshes, P. Kinney, Presented at the 6^{th} International Meshing Roundtable (1997). Available from http://www.imr.sandia.gov/papers/imr6/kinney97.ps.gz
[11] Redesign of the Paving Algorithm: Robustness Enhancements through Element by Element Meshing, D. White & P. Kinney, Presented at the 6^{th} International Meshing Roundtable (1997). Available from http://www.imr.sandia.gov/papers/imr6/white97.ps.gz
[12] LCCWin32 a free compiler system for Windows http://www.cs.virginia.edu/~lccwin32/9 Glossary
ASCII American Standard Code for Information Interchange numeric representation of alphanumeric symbols.CAD Computer Aided Design
DXF Drawing eXchange Format
GLU OpenGl Utilities
GLUT OpenGL Utility Toolkit, simple interface for windowed OpenGL graphics on Win32 and XWindows systems.
NURBS NonUniform Rational BSpline
OpenGL A standardised 3D graphics programming interface. Allows low level access to compatible graphics devices for fast graphics.
SPHIGS Simple Programmers Hierarchical Graphics Standard.
SRGP Simple Raster Graphics Package.