Domain Decomposition of Geometrical Objects - Semi-Automatic Mesh Generation Contents

1 Abstract *

2 Introduction *

2.1 Finite Elements *

2.2 DXF Files *

3 Algorithm Requirements * 3.1 DXF files *

3.2 Visualisation of Results *

3.3 Automatic Meshing *

4 DXF to Boundary Converter * 4.1 Reading Data from DXF files *

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 Command-Line Options *

4.11 Other Developments *

4.11.1 Debugging Information *

4.11.2 Boundary Viewer *

5 Paving of Boundaries * 5.1 Loading Objects *

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 Results * 6.1 Breakout *

6.2 Asymmetry *

6.3 Sensitivity to Segment Lengths *

6.4 Performance *

7 Conclusions *

8 References *

9 Glossary *

1 Abstract

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.

2 Introduction

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.

2.1 Finite Elements

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 man-years. 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.

2.2 DXF Files

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 non-uniform rational B-spline (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.

3 Algorithm Requirements

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.

3.1 DXF files

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.

3.2 Visualisation of Results

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 X-Windows platforms.

3.3 Automatic Meshing

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:

Of these algorithms the Paving algorithm seems to be the simplest to implement, as it only requires vector operations and will directly produce the quad mesh.

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.

4 DXF to Boundary Converter

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 self-intersections. These characteristics are summarised in Fig. 1.

A brief overview of the processes carried out is as follows:

Figure 1 Acceptable Boundaries

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 sub-divided 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 co-ordinate 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.

4.2 Split Overlapping Objects

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

4.3 Linearise 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 End-Points

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.

4.4 Form Closed Boundaries

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.

4.5 Check Boundary

Each boundary is checked for self-intersections by checking that no segment intersects any other segment, using a parametric method, derived from the Cyrus-Beck line-clipping algorithm [7].

Two intersecting vectors are shown in Fig. 4 and illustrate the method. Where two lines (P1 & P2) intersect the dot product of the normal to one line (N2) and the vector formed by subtracting a point on the line (P22) from it’s intersector (P1) 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 self-intersections are found the boundary is simple, a non-simple boundary could escape the next stage erroneously.

4.6 Adjust Order

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.

4.7 Hierarchy Determination

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&sup2; 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 ray-casting 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.

4.8 Output Boundaries

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 co-ordinates, 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.

4.9 Viewer Module

The viewer module may be started if the relevant pre-processor 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 X-Windows 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 sub-menu. 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.

4.10 Command-Line Options

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 floating-point 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.

4.11 Other Developments

4.11.1 Debugging Information

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 pre-processor 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

5 Paving of Boundaries

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:

As described in the original paper the paving algorithm makes use of the program flow in Fig. 12. This gives the outline of the process, it is a little more complicated in practice, as control may have to jump from one section to another.

Each operation is described below

Repeat Load an object
Repeat Select a boundary to add a layer to Repeat Add a layer section
Smooth layer section
Seam layer
Check for intersecting layers Join intersecting layers
Seam layers
Until full layer added
Adjust layer
Check for intersecting layers
Join intersecting layers
Seam layers
Closure of suitable boundaries
While boundaries left in shape
Cleanup mesh
While objects left in file

Figure 12 Paving Flow Diagram

5.1 Loading Objects

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.

5.2 Adding Layer Sections

Vertices on a boundary can be classified by the angle (q ) at the node Ni 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.

5.3 Smoothing

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 non-fixed nodes that are not part of a boundary.

5.3.1 Boundary Node Smoothing

The isoparametric smooth is found by taking a quad that the point being smoothed (Ni) is in. The vectors from the origin to the vertices either side of the point (Nj & Nl) are summed and then the vector from the origin to the vertex opposite (Nk) 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 (Va) from Ni 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 (Ni-Nj-Va), to the length from Ni to Nj (referred to as the desired distance lD), to give a new vector Vb, 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 Nj and the opposing vertices of the two adjoining quads is found (Vbis1), as illustrated in Fig. 20. This is then used to find the vector bisecting it and the vector from Nj to Ni (Vbis2). The length of this vector (lQ) is the distance from Nj to it's intersection with a line joining the opposing vertices in the adjoining quads (Q). The length is adjusted by comparing lQ and lD, if lD is longer than lQ the magnitude is adjusted to the average of lD and lQ, otherwise it is the same as lD.

Figure 20 Angular Smoothness

This gives an angular adjustment vector of:

Vc = Vbis2 - Vi

These are then averaged to give a boundary smoothing vector of:

VBND = (Vb + Vc)/2

5.3.2 Internal Node Smoothing

Internal node smoothing is simpler and is demonstrated in Fig. 21. For each of the N nodes connected to the node being smoothed (N1-N4) a vector is formed of the vector from the node being smoothed (Ni) to the node in question (V1-V4). 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 (Va) 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, VL as in Equation 8.

Figure 21 Laplacian Smooth

5.3.3 Constrained Smoothing

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 re-evaluated, an arbitrary number of times.

5.4 Seaming

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.

P11-P12 and P21-P22 are the endpoints of the intersecting segments, none of which are co-incident.
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

5.6 Row Adjustment

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. d2C&times;d1 and d4C&times;d3. 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. d2<C&times;d1. 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

5.7 Closing Boundaries

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

5.8 Cleanup Mesh

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.

5.9 Writing Meshes

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 co-ordinates. nPoints entries.


Q1 Q2 Q3 Q4 Qstring                                 Q1-Q4 Zero based indices of points forming quad, in anticlockwise order. Qstring is a zero-based index of property string, if no property string is applicable this should be –1. nQuads entries must occur.

5.10 Command Line Options

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.

5.11 Other Developments

5.11.1 Debugging Information

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.

5.11.2 Connection Viewer

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

5.11.3 Mesh 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

5.11.4 Settings File

During initialisation a settings file may also be specified, by a second command line argument, to allow variation of some algorithm parameters, otherwise default settings are used. The settings file contains information used in wedging, tucking, seaming, joining and smoothing operations. The details are contained in an ordered ASCII text file. The following is a commented example of the default settings:
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)
6 Results
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 non-optimal 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
[2] NURBS Book (Second Edition), L. Piegl & W. Tiller, Springer-Verlag (1997), ISBN 3-540-61545-8.
[3] MESA an OpenGL conformant graphics library, source code available from
[4] A Survey of Unstructured Mesh Generation Technology, S. Owen, available at
[5] CUBIT Toolkit for robust and unattended mesh generation
[6] Automatic mesh generation: application to finite element methods, P.L. George, Wiley (1991), ISBN 0-471-93097-0
[7] Introduction Computer Graphics, J. Foley et al, Addison-Wesley (1990), ISBN 0-201-60921-5
[8] Paving: A New Approach to Automated Quadrilateral Mesh Generation, T. Blacker & M. Stephenson, International Journal for Numerical Methods in Engineering, Vol 32 811-847 (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 7th International Meshing Roundtable (1998). Available from
[10] CleanUp: Improving Quadrilateral Finite Element Meshes, P. Kinney, Presented at the 6th International Meshing Roundtable (1997). Available from
[11] Redesign of the Paving Algorithm: Robustness Enhancements through Element by Element Meshing, D. White & P. Kinney, Presented at the 6th International Meshing Roundtable (1997). Available from
[12] LCCWin32 a free compiler system for Windows
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 X-Windows systems.

NURBS Non-Uniform Rational B-Spline

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.