Some algorithms to correct a geometry in order to create a finite element mesh

==SOME ALGORITHMS TO CORRECT A GEOMETRY IN ORDER TO CREATE A   FINITE ELEMENT MESH==

'''R. Ribó, G. Bugeda,   and E. Oñate'''

{|
|-
|International Center for Numerical Methods in Engineering (CIMNE)
|-
| Universidad Politécnica de Cataluña;   Campus Norte UPC; 08034 Barcelona; Spain
|-
| web page: [http://cimne.upc.es/ http://cimne.upc.es/]
|}

{|
|-
|fax.: +34 93 4016517, e-mail: [mailto:bugeda@cimne.upc.es bugeda@cimne.upc.es]
|}

==Abstract==

One of the major difficulties in meshing 3D complex     geometries is to deal with non proper geometrical definitions coming     from CAD systems. Typically, CAD systems do not take care of the     proper definition of the geometries for the analysis purposes.  In     addition, the use of standard CAD files (IGES, VDA,...)  for the     transfer of geometries between different systems introduce some     additional difficulties.      In this work, a collection of algorithms to repair and/or to improve   the geometry definitions are provided. The aim of these algorithms is   to make as easy as possible the generation of a mesh over complex   geometries given some minimum requirements of quality and correctness.   The geometrical model will be considered as composed of a set of NURBS   lines and trimmed surfaces.       Some examples of application of the algorithms and of the meshes   generated from the corrected geometry are also presented in this work.

'''keywords''' Mesh generation, geometry correction, NURBS.

==1 INTRODUCTION==

Mesh generation using structured and unstructured grids is still   nowadays one of the bottlenecks in the practical application of the   finite element method (FEM). See <span id='citeF-1'></span>[[#cite-1|[1]]] and <span id='citeF-2'></span>[[#cite-2|[2]]].

Geometrical models created for finite element mesh generation purposes   are typically created by the mesh generation software or imported from   an external CAD system.  In both cases, the geometrical model must   satisfy a series of quality constraints.

Normally, geometrical models are constructed as a set of NURBS lines   and trimmed surfaces. These entities can be mathematically correct and   suitable for other uses like visualization or CAM but, in many cases,   they are not good enough for meshing operations. In these cases, it   is necessary to adapt/repair the geometrical entities by changing   their mathematical description while maintaining the same geometrical   shape.  The final entities should be better suited for the mesh   generation operations. See <span id='citeF-3'></span>[[#cite-3|[3]]] and <span id='citeF-4'></span>[[#cite-4|[4]]].

The algorithms proposed in this paper are considered as a set of   filters that will be applied to a geometry definition just after being   imported or created. The algorithms are:

* collapse of entities
* correction of the list of ''knots''
* reparametrization
* conversion to similar cubic entity
* union of curves
* reorientation of the boundary of the surfaces
* collapse of small angles

With the help of these algorithms, it is possible to automate the   process of importing geometry from CAD and mesh on it while maintaining an   acceptable level of quality criteria. This process should be performed   without the need of manual intervention.

To understand better what follows, a short description of the NURBS lines   and surfaces entities is given. A deeper description can be found in [] farin-cadg and <span id='citeF-5'></span>[[#cite-5|[5]]].

===1.1 NURBS lines===

A NURBS line is a geometrical entity, described as a parametric line in the   3D space, that is defined with a set of ''knots'', a set of ''control   points'', and a set of ''weights'' if it is rational.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-bspline-quadratic-example.png|511px|Example of a quadratic NURBS line with 4 control points evaluated   for a given value of u.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Example of a quadratic NURBS line with 4 control points evaluated   for a given value of <math>u</math>.
|}

The ''knots'' are a list of non-decreasing numbers <math display="inline">u_0,\dots ,u_{L+n}</math>   that begin in the lower range of the parameter space (usually 0.0) and   finish in the high range of it (usually 1.0). A ''knot'' has   multiplicity <math display="inline">r</math> if its value is repeated <math display="inline">r</math> times inside the list of '' knots''. The initial and the end knots must have multiplicity <math display="inline">r=n+1</math>, where <math display="inline"> n</math> is the degree of the curve.

The ''control points'' <math display="inline">\vec{P}_1,\dots ,\vec{P}_L</math> are a list of points   in the 3D space that are part of the NURBS definition. The line will   interpolate the first and the last point and will smoothly approximate the   other points.

The ''weights'' <math display="inline">\omega _1,\dots ,\omega _L</math> are a set of non-negative real   numbers, one for every control point. They allow the shape of the   curve to change and also to have an exact representation of conic curves like circles   or ellipses.

The number of ''knots'', <math display="inline">N_k</math> must be equal to <math display="inline">N_k=N_p+n+1</math>, where <math display="inline"> N_p</math> is the number of control points.

The evaluation of a NURBS curve for a given value of the parameter <math display="inline">u</math> can be   done in a recursive manner. First, it is necessary to identify the interval <math display="inline"> [u_I,u_{I+1}]</math> in the list of knots that contains the value of <math display="inline">u</math> (<math display="inline">u\in [ u_I,u_{I+1}]</math>). Next, the following recursive expression can be applied:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vec{P}_i^k(u)=\frac{(1-\alpha _i^k) \omega _{i-1}^{k-1}\vec{P}_{i-1}^{k-1}(u)+ \alpha _i^k\omega _{i}^{k-1}\vec{P}_{i}^{k-1}(u)}{ \omega _{i}^{k}} \quad \hbox{with: }\left|   \begin{array}{l}k = 1,\dots ,n-r \\    i = I-n+k+1,\dots ,I+1 \\    \vec{P}_{i}^{0}(u) = \vec{P}_{i} \\    \omega _i^0=\omega _i   \end{array}   \right.  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|-
| style="text-align: center;" | <math> \alpha _i^k=\frac{u-u_{i-1}}{u_{i+n-k}-u_{i-1}}  </math>
|-
| style="text-align: center;" | <math> \omega _i^k=(1-\alpha _i^k)\omega _{i-1}^{k-1}+ \alpha _i^k\omega _{i}^{k-1}    </math>
|-
| style="text-align: center;" |

then:

<math>   \vec{s}(u)=\vec{P}_{I+1}^{n-r}(u)   </math>
|}
|}

===1.2 NURBS surfaces===

A NURBS surface is the extension of the NURBS line to one additional   dimension in the parametric space. Most of the properties of the NURBS   curves applies here. There is a ''list of knots'' for every parametric   direction <math display="inline">u</math> and <math display="inline">v</math>. The ''control points'' are a set of <math display="inline">P_{ij}</math>   points with <math display="inline">i\in \lbrack 1,\dots ,L_{u}]</math> and <math display="inline">j\in \lbrack 1,\dots ,L_{v}]</math> .

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-nurbs-surface.png|493px|NURBS surface with the corresponding control polygon.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' NURBS surface with the corresponding control polygon.
|}

The evaluation of the surface can be made in different ways:

*  To evaluate first the NURBS curve corresponding to one of the   parametric directions <math display="inline">u</math> (maintaining <math display="inline">v</math> constant) and obtain a NURBS   line. Then, to evaluate the resulting NURBS curve in the second direction <math display="inline">v</math> .

*  Recursive evaluation. Similar to the expression (1) given for the   NURBS curve.

*  Evaluation by means of the ''b-spline'' base. This technique is   also usable for curves.

The evaluation with the ''b-spline'' base is done by defining a base of    ''b-splines'', <math display="inline">N_i</math>, as a recursive function in the following recursive   way:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>N_i^0(u)=\left\{   \begin{array}{l}1\qquad \hbox{if }u_{i-1}\leq u < u_i \\    0\qquad \hbox{if not.}   \end{array}   \right. </math>
|-
| style="text-align: center;" | <math> N_i^n(u)=\frac{u-u_{i-1}}{u_{i-n-1}-u_{i-1}}N_i^{n-1}(u)+ \frac{u_{i+n}-u}{ u_{i+n}-u_{i}}N_{i+1}^{n-1}(u)  </math>
|-
| style="text-align: center;" |

Then, the following expression is applied:

<math>   \vec{s}(u,v)=\frac{ \sum _i\sum _j\omega _{ij}\vec{P}_{ij}N_i^m(u)N_j^n(v) }{ \sum _i\sum _j\omega _{ij}N_i^m(u)N_j^n(v)}     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

==2 COLLAPSE OF ENTITIES==

The standard CAD exchange formats (IGES, VDA,...) do not contain the   topological connection between the different surface patches defining a   closed geometry. These files just contain the mathematical description of   the patches. Nevertheless, in order to proceed with a correct mesh   generation it is necessary to check the correctness of the 3D geometry by   ensuring, for example, that it is defined by a completely closed surface. In addition,   many mesh generation algorithms need to know the neighboring relation between patches   in order to advance in their work. For this reason, it is necessary to obtain the topological connection   between the different surface patches.

An important additional difficulty is that, very often, the boundary curves   defining neighboring patches are very similar, but not identical. Normally,   the corresponding pairs of curves are close enough for the different   visualization or CAM operations, but not for the necessary operations for   the analysis. This makes  the task of identifying neighbor   patches very difficult.

In order to solve this difficulty it is necessary to define a geometrical   tolerance. When two different geometrical entities are separated by a   distance smaller than the prescribed tolerance they are considered to be   identical and they are collapsed. This tolerance can be specified by the   user or can be obtained in an automatic way as a certain percentage of the   total size of the geometrical model.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-picture-94d646.png|600px|The collapse procedure reduces the initial four curves to only   three.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' The collapse procedure reduces the initial four curves to only   three.
|}

Given a geometrical tolerance <math display="inline">\epsilon </math>, the collapse criteria for the   different geometrical entities are the following:

*  Two points are collapsed into a single one if:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\Vert \vec{P}_{1}-\vec{P}_{2}\Vert <\epsilon    

</math>
|}
|}

where <math display="inline">\vec{P}_i</math> is the vector of coordinates of point <math display="inline">i</math>.

*  Two curves are collapsed if they share their end points and the   maximum distance between them is smaller than <math display="inline">\epsilon </math>. In addition, they   are also collapsed if a curve belongs to the interior of another   one. This last property means that, in addition to the first criterion, its   end points <math display="inline">\vec{P}_{i}</math> accomplish

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathrm{dist}(\vec{P}_{i},L)<\epsilon    

</math>
|}
|}

where <math display="inline">\vec{L}(t)</math> is the second curve. In the latter case, the final result   of the collapse operation is not a single curve but a group of curves   as shown in Figure 3.

*  Two patch surfaces are collapsed if they share their boundary curves   and a similar maximum distance criterion is accomplished.

*  Two volumes are considered as equivalent if they share their   boundary surfaces.

The computation of the maximum distance between different curves or   surfaces is made in an approximated way as follows: a fixed number of   points <math display="inline">\vec{P}_i</math> are chosen in the interior of the curve or surface   depending on their geometrical characteristics. The distance   point-curve or point-surface <math display="inline">d_i</math> is computed for each of the   selected points. Then, <math display="inline">d_{max}</math> is approximated as <math display="inline">d_{max}=max(d_i)</math>.

In general, a big number of control points and a high degree of the line or surface   formulations will imply a bigger geometrical complexity. For this reason,   the number of selected points will be related to the two mentioned values.

In the computation of the <math display="inline">d_i</math> values the <math display="inline">\vec{P}_i^{\prime }</math> point   is obtained as the ''mapping'' of the original point over the   geometrical entity and <math display="inline">d_i=\| \vec{P}_i-\vec{P}_i^{\prime }\| </math>. The   ''mapping'' operation is the transformation of an original point to   the closest point contained into a geometrical entity like a curve or   surface. This procedure is described in&nbsp;<span id='citeF-5'></span>[[#cite-5|[5]]].

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-picture-9f224e.png|600px|Saving of elements when unnecessary details are eliminated.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Saving of elements when unnecessary details are eliminated.
|}

The collapse operation is also useful to eliminate small details of the   geometry that can be of no relevance for the analysis. If some lines or surface   patches are smaller than the tolerance, they can disappear and their   neighbor entities are reconnected. This operation is normally called '' feature reduction''. From the practical point of view, many details that are   crucial for the design task, can be totally irrelevant for the analysis.   This operation can save a lot of elements and degrees of freedom in the   final finite element mesh. In Figure&nbsp;[[#img-4|4]] a graphical   example of this possibility is shown.

==3  CORRECTION OF THE LIST OF    ''KNOTS''==

The list of ''knots'' is a set of non decreasing real values that   belong to the parametric space and are used for the definition of a   NURBS. A more detailed definition of ''knot'' can be obtained   in&nbsp;<span id='citeF-6'></span>[[#cite-6|[6]]].

In the formulation of a NURBS curve or surface the list of ''knots''   cannot contain any ''knot'' with a multiplicity higher than the degree of   the mathematical entity. Nevertheless, in many cases the standard exchange   files contain this type of error. In these cases, it is necessary to correct   the list of ''knots'' by eliminating the overdefined ''knot'' and the   associated weights and points. In this way, if we have a <math display="inline">n</math> degree NURBS   curve defined by <math display="inline">L</math> points with the following list of ''knots'':

<span id="eq-7"></span>
<span id="eq-8"></span>
<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_{0},\dots ,k_{i},k_{i+1},\dots ,k_{i+n+p},\dots ,k_{L+n}\qquad \hbox{with:}\quad k_{i}=k_{i+1}=\dots =k_{i+n+p}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|-
| style="text-align: center;" | <math> \vec{P}_{1},\dots ,\vec{P}_{L}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|-
| style="text-align: center;" | <math> \omega _{1},\dots ,\omega _{L}     </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}
|}

the list ([[#eq-7|7]]) will be reduced to:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_{0},\dots ,k_{i},k_{i+1},\dots ,k_{i+n},k_{i+n+p+1},\dots ,k_{L+n}   </math>
|}
|}

and the list of control points ([[#eq-8|8]]) and weights ([[#eq-9|9]]) to:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vec{P}_{1},\dots ,\vec{P}_{i+1},\vec{P}_{i+p+2},\dots ,\vec{P}_{L} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|-
| style="text-align: center;" | <math> \omega _{1},\dots ,\omega _{i+1},\omega _{i+p+2},\dots ,\omega _{L}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}
|}

Note that the new list will have <math display="inline">L-p</math> control points and that the   resulting curve will be an approximation to the original one, being   the original one a mathematically incorrect representation of a curve.

==4 REPARAMETRIZATION==

A very desiderable property for the mesh generation algorithm is   to use the arc length as the defining parameter over the curve or surface. This   property is not normally available in the NURBS curves or surfaces.   This can be partially corrected by using a constant distance between   each pair of consecutive control points. Nevertheless, in order to   obtain the best possible quality during the mesh generation process it   is convenient to improve the parametrization from the very beginning.   The corresponding changes that are produced in the curves or surfaces   during the improvement of their parametric definition are named   reparametrization.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-picture-84767d.png|600px|Comparison between the derivatives of the curve before and after   the reparametrization.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Comparison between the derivatives of the curve before and after   the reparametrization.
|}

A major problem arises when there are big discrepancies between the   spacing of distances between the control points and the spacing between the orresponding   associated points of the curves or surfaces. This problem can be detected by   a comparison between the modulus of the derivatives in the following way:

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{1}{F}< \frac{\left\|\left.\frac{d\vec{L}}{dt}\right|_{t=k_{i-}}\right\| }{\left\|\left.\frac{d\vec{L}}{dt}\right|_{t=k_{i+}}\right\|} < F      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where the <math display="inline">k_i</math> in ([[#eq-12|12]]) represents an interior ''knot''    and <math display="inline">F</math> is a maximum acceptable parameter that,    for acceptable parametrizations, can have values between 4 and 5. In   this formula, <math display="inline">\vec{L}(t)</math> represent the curve.

The first step of the correcting process consists of decomposing the   curve into the addition of the set of the equivalent Bézier curves.   This can be done by inserting multiple ''knots'' until a   multiplicity equal to the order of the mathematical entityi is reached. The process of   inserting ''knots'' is described in <span id='citeF-6'></span>[[#cite-6|[6]]].

From the geometrical and the parametrization points of view, the defined curve is completely equivalent to the original one.   If we have a <math display="inline">n</math> degree curve and <math display="inline">L</math> defining points with the following values:

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_0,\dots ,k_i,\dots ,k_{L+n}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|-
| style="text-align: center;" | <math>   \vec{P}_1,\dots ,\vec{P}_L  </math>
|-
| style="text-align: center;" |

new ''knots'' are inserted with the same value of the       already existing ''knots'' into ([[#eq-13|13]]) until they all       have a multiplicity of the same order as the curve

<math>     \underbrace{k_0,\dots ,k_0}_{n+1},\dots ,\underbrace{k_i,\dots ,k_i}_{n+1},     \dots ,\underbrace{k_{L+n},\dots ,k_{L+n}}_{n+1}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}
|}

After inserting these knots, the number of control points of the curve   will increase. The new number of points will Bézier equal to the sum   of the difference between the original multiplicities and   the order of each of the original ''knots''.

Taking the list of points and grouping them in sets of <math display="inline">n+1</math> points,   it is possible to buid successive Bézier curves with the same degree   than the original curve and the same continuity between successive   Bézier curves. These new curves will not depend on the list of   ''knots'' and, therefore, their shape will not change if any of the   ''knots'' is moved. Consequently, the increments between each   groups of ''knots'' can be recomputed in order to obtain a better   parametrization. A possible modification based in the ''chord   length'' parametrization consists in creating the following new list:

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_{0},\dots ,k_{n},\dots ,k_{i},\dots ,k_{i+n},\dots ,k_{L^{\prime }},\dots  ,k_{L^{\prime }+n} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|-
| style="text-align: center;" | <math> k_{0}=\dots =k_{n}\qquad k_{i}=\dots =k_{i+n}\qquad k_{L^{\prime }}=\dots  =k_{L^{\prime }+n} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|-
| style="text-align: center;" | <math> k_{i+n+1}-k_{i+n}=\frac{l_{c}}{l_{tot}}     </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}
|}

where <math display="inline">l_{c}</math> in ([[#eq-17|17]]) is the length of the corresponding Bezier   curve and <math display="inline">l_{tot}</math> is the total length.

The new curve parametrized in this way will have the same geometrical shape   than the original one but with a different parametric definition. Figure&nbsp; fig:conv-to-bezier shows the typical improvement than can be obtained with   the described algorithm.

==5 CONVERSION TO SIMILAR CUBIC ENTITY==

In some cases, the correction described in the last section is not   enough for ensuring a good parametrization. Typically, these cases   arise when one or more control points are repeated, or else when the   orders of magnitude of the distances between points have big   discrepancies between them.

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-picture-0f2f3a.png|600px|]]
|[[Image:draft_Ribó_247186029-picture-e9d3f9.png|400px|Conversion of a NURBS to an approximated interpolant cubic NURBS.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 6:''' Conversion of a NURBS to an approximated interpolant cubic NURBS.
|}

In these cases it is acceptable to use a new curve or surface not   identical to the original one but with a good enough approximation. We   should keep in mind that the exchange of CAD files is always carried   out with different geometrical tolerances. Hence, it seem logical to   allow geometrical changes below the limits of this   tolerance.

The process consists of computing a set of points <math display="inline">\vec{P}_{i}</math>   belonging to the interior of the curve. These points will be the base   of an interpolation algorithm that will be used for the definition of   the new curve with the required approximation to the original one. The   criteria for the selection of the number of points <math display="inline">L_{I}</math> is obtained   by correlating the number of control points <math display="inline">L</math> with the accepted   tolerance. A new cubic curve can be interpolated from the new set of   points using different procedures like the one described in   <span id='citeF-6'></span>[[#cite-6|[6]]]. Figure&nbsp;[[#img-6|6]] shows an example of   this type of conversion. In this case, the discrepancy between both   curves is high due to a point with <math display="inline">C^{0}</math> continuity.

==6 UNION OF CURVES==

It is possible to convert a set of connected NURBS curves into a single   curve.  Normally, it is advantageous to keep the number of curves as   small as possible because it can make the mesh generation task easier.   Typically, the presence of extremely small segment lines can increase   the complexity of the mesh generation.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-picture-e2d0fb.png|600px|Different criteria for accepting (or not) the union of curves.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Different criteria for accepting (or not) the union of curves.
|}

The criteria for joining different curves are the following:

<ol>

<li> Enough degree of continuity between the different segments must   exist. Typically, a <math display="inline">C^{1}</math> continuity is considered as enough.     </li>

<li> None of the segments can support any individual surface patch. The   set of segments to be joined must belong to the boundaries of the same   surface patches.   </li>

</ol>

Figure&nbsp;[[#img-7|7]] shows different examples of these situations.

Before the union can proceed, the different curves must have the   same degree. Hence, the degree of all the segments will be increased   until they reach the maximum value <math display="inline">n</math>. Next, the different control   points will be joined and a new list of ''knots'' will be computed   starting from the original ones in the following way:

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_{A,0},\dots ,k_{A,L_A+n}\qquad \hbox{curve A} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|-
| style="text-align: center;" | <math> k_{B,0},\dots ,k_{B,L_B+n}\qquad \hbox{curve B} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|-
| style="text-align: center;" |

The new list will be:

<math>   \alpha k_{A,0},\dots ,\alpha k_{A,L_A+n},\alpha k_{A,L_A+n}+ \beta  k_{B,0},\dots ,\alpha k_{A,L_A+n}+\beta k_{B,L_B+n} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|-
| style="text-align: center;" | <math> \alpha =\frac{l_A}{l_A+l_B}\qquad \beta =\frac{l_B}{l_A+l_B}     </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}
|}

where <math display="inline">l_A</math> and <math display="inline">l_B</math> in ([[#eq-21|21]]) are the respective lengths of the   curves.

==7 SUBDIVISION OF CURVES==

Some of the algorithms related with the mesh generation task solve   small non linear systems of equations that involve the derivatives of   the analytical expression defining the curves. The presence of   strong changes of these derivatives inside the curve can produce   convergence difficulties in the solution of the mentioned systems. Due   to this reason, it is very convenient to maintain a <math display="inline">C^{1}</math> continuity   over the whole curve (In practice, small angle discontinuities can   also be allowed). Hence, it is convenient to subdivide (cut) the   original curves at points not satisfying the required   continuity.

This subdivision is made through the insertion of additional ''knots'' at   the required cutting points until a multiplicity equal to the order   of the curve is reached. The new curves will be:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_{0},\dots ,\underbrace{k_{i},\dots ,k_{i}}_{n+1},\dots ,k_{L+n}    </math>
|-
| style="text-align: center;" | <math> \frac{k_{0}}{k_{i}},\dots ,\frac{k_{i}}{k_{i}}\qquad  \hbox{ curve A}   </math>
|-
| style="text-align: center;" | <math> \frac{k_{i+n+1}-k_{i+n+1}}{k_{L+n}-k_{i}},\dots ,\frac{k_{L+n}-k_{i+n+1}}{ k_{L+n}-k_{i}}\qquad  \hbox{ curve B}   </math>
|}
|}

The control points will be simply spread depending on the number of '' knots'' corresponding to each curve.

==8 REORIENTATION OF THE BOUNDARY OF THE SURFACES==

Some of the mesh generation algorithms require a proper orientation of all   the boundary curves. A curve <math display="inline">\vec{L}(t)</math> that belongs to the boundary of a   surface with a normal vector <math display="inline">\vec{N}</math> is considered as well oriented if the   vector <math display="inline">\vec{V}</math> defined by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vec{V}=\vec{N}\times \frac{d\vec{L}}{dt}    </math>
|}
|}

always points towards the interior of the surface.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-picture-983f81.png|600px|Criterion of signs and orientations of the trimming lines of a   trimmed surface.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Criterion of signs and orientations of the trimming lines of a   trimmed surface.
|}

In some cases, the information imported from a CAD system does not   satisfy the mentioned criterion. It is then convenient to check this   possibility and to change the corresponding orientation when necessary.

Nevertheless, some times it is not easy to identify the interior of   a surface. In this work, we use the fact that for non trimmed surfaces the   curves must belong to the boundary of a NURBS surface. For   trimmed surfaces, it is necessary to look for any of the trimming curves   placed on the surface.

A boundary line of a NURBS surface patch is well oriented if:

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{d\vec{L}(t)}{dt}\cdot \frac{\delta \vec{S}(u,v)}{\delta u}>0 \qquad    \vec{L}\in v=0   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|-
| style="text-align: center;" | <math> \frac{d\vec{L}(t)}{dt}\cdot \frac{\delta \vec{S}(u,v)}{\delta u}<0 \qquad    \vec{L}\in v=1   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}
|}

where <math display="inline">\vec{S}(u,v)</math> is the parametric surface. If the curve lays on <math display="inline">u=0</math> or <math display="inline">u=1</math> the corresponding similar expressions will be used.

An additional convenient check consists of computing the normal   vector to the surface in its center and to compute the approximated   normal vector to the boundary curves. The approximate normal vector   can be computed as:

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vec{N}_{L}=\sum _{i}\vec{L}(t_{i})-\vec{C}\times \vec{L}(t_{i}+\Delta t_{i})- \vec{C}     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

where <math display="inline">\vec{C}</math> is the center of the surface. The boundary curves will be   well oriented if:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vec{N}\cdot \vec{N}_{L}>0   </math>
|}
|}

==9 COLLAPSE OF SMALL ANGLES==

A very common problem in files imported from CAD systems is that some   of the surfaces contain almost null angles that make impossible the   generation of a proper mesh. In those cases, it is convenient to   collapse part of the lines that define these angles until converting   the angles into bigger ones allowing to place acceptable mesh elements.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-picture-cb47e3.png|600px|Collapse of a too small angle. The collapse will be accepted if L   is bigger than a minimum value.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Collapse of a too small angle. The collapse will be accepted if <math>L</math>   is bigger than a minimum value.
|}

The collapsing angle will be computed depending on the given tolerance <math display="inline"> \epsilon </math>. In general, the resulting criterion should guarantee that   the size of the resulting curve is in agreement with the rest of the   contiguous curves (see Figure&nbsp;[[#img-9|9]]).

==10 EXAMPLES==

The algorithms presented in previous sections have been implemented   into the GiD pre/postprocessing system<span id='citeF-7'></span>[[#cite-7|[7]]] developed at   CIMNE. The following examples are a set of representative applications   of these algorithms for the preparation of different geometries in   order to be meshed. Typically, the geometry has been defined using a   CAD system and the corresponding information has been introducind into   GiD using standard CAD files (IGES, VDA, ...). In all cases, the use   of the presented algorithms has carried out the   improvements/reparations needed to allow meshing of the geometry   without the necessity of any additional operation.

====10.1 Geometry of a set of casting teeth====

This example (see Figures&nbsp;[[#img-10|10]] and&nbsp;[[#img-11|11]]),   corresponds with the generation of a finite element mesh over the solid   model of a casting set of teeth. The number of surface patches used for the geometry   definition is around 400. In this particular case, before the mesh generation it has been necessary   to correc the geometry in order to create a closed volume (without gaps). The final obtained mesh has   around 40,000 elements.

====10.2 Tarazona cathedral====

Figures&nbsp;[[#img-12|12]] and&nbsp;[[#img-13|13]] represent a structural analysis model   with a non-linear damage material model of the Tarazona cathedral.    In this case, it has been necessary to correct many of the surfaces used   to describe the fine details of the shape in some parts of the   geometry. The final mesh has around 40,000 elements.

====10.3 Sheet stamping processes====

This is a typical case corresponding with a sheet stamping analysis (see Figure&nbsp;[[#img-14|14]]). The    geometry has been modeled using traditional   CAD systems and before proceeding to the classical mesh generation operations   it has been necessary to use all the   correction algorithms presented in these pages. Some of the problems   that usually appear in this type of geometries are:

* Some surfaces are too small and must disappear.
* Other surfaces have a bad mathematical description.
* Some of them are correctly defined but with some properties that     are not suitable for meshing on them.
* Some of the geometrical details are very small and produce extremely complex geometries.

Figure&nbsp;[[#img-14|14]] shows a typical mesh for this type of cases with around 400,000 finite elements.

====10.4 Aluminium casting====

These type of geometrical models have the same problems as the previous one.   However, they are more complicated to improve due to the   inherent complexity of their surface patches. See   Figures&nbsp;[[#img-15|15]] and&nbsp;[[#img-16|16]].

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-dents-geom.png|520px|Geometry of a set of casting teeth for construction       machines. This set is the one used in the casting process.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Geometry of a set of casting teeth for construction       machines. This set is the one used in the casting process.
|}
<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-dents-render.png|536px|Fotorrealistic render of the previous geometry.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Fotorrealistic render of the previous geometry.
|}

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-geom_catt.png|783px|Geometry of the Tarazona cathedral (Spain).]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Geometry of the Tarazona cathedral (Spain).
|}

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-malla_catt.png|783px|Mesh of the cathedral.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Mesh of the cathedral.
|}

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-estampacio.png|518px|Analysis of sheet stamping processes.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Analysis of sheet stamping processes.
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-injeccio-plastics-s.png|480px|Mesh of an aluminium casting model.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Mesh of an aluminium casting model.
|}

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Ribó_247186029-injeccio-plastics-2s.png|480px|Detail of the previous mesh.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Detail of the previous mesh.
|}

==11 CONCLUSIONS==

CAD models are not always best suited as starting points for mesh generation.   CAD models and the corresponding exchange files are more   suited for design and CAM operations than for the analysis task. For this   reason, it is normally necessary to repair the CAD models in order to make   them suitable for mesh generation.

The proposed algorithms have shown to be efficient in dealing with three   different problems:

*  The elimination of some small features that are necessary for the   design but are not important for the numerical analysis.

*  The correction of some mistakes and accuracies that are found in many   geometrical interchange files.

*  The change in the mathematical definition of some entities in order   to make them better suited for mesh generation.

All these algorithms can be executed automatically as a filter to the   imported or created geometry reducing the human interaction to a minimum and   facilitating a lot the mesh generation task.

An academic version of the GiD pre/postprocessing system, where the   algorithms presented in this paper have been implemented, can be   freely downloaded at

http://gid.cimne.upc.es

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' R. Löhner and P. Parikh. (1988) "Generation of three-dimensional unstructured grids bythe advancing front method". Int. j. numer. methods fluids. 8

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' P. L. George. (1991) "Automatic Mesh Generation. Application to Finite Element Methods". Wiley

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' J M Le Gouez and M A Boheas and C J Miles. (1995) "Data requirements for analysis modelling". NAFEMS

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' NAFEMS. (1996) "Integrate CAD and Analysis"

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' Ramon Ribó. (2000) "Development of an integrated system for the geometry dealing, mesh generationand data for the finite element analysis"

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Gerald Farin. (1973c1993) "Curves and surfaces for CADG". ACADEMIC PRESS, INC., Third Edition

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Ramon Ribó and others. (1998) "GiD reference manual". CIMNE

<div id="cite-8"></div>
'''[8]''' D.N.Knuth. (1973) "The art of computer programming, Vol. 3". Addison-Wesley

<div id="cite-9"></div>
'''[9]''' Hoschek, Josef and Dieter Lasser. (1993) "Fundamentals of computer aided geometric design". Addison-Wesley

<div id="cite-10"></div>
'''[10]''' Foley and van Dam and Feiner and Hughes. (1993) "Computer Graphics. Principles and practice". Addison-Wesley, second Edition

<div id="cite-11"></div>
'''[11]''' P. Lancaster and K. Salkauskas. (1988) "Curve and Surface Fitting". Academic press, first Edition

<div id="cite-12"></div>
'''[12]''' O. C. Zienkiewicz and R. Taylor. (2000) "The Finite Element Method", Volume I. John Wiley and Sons, 5 Edition

<div id="cite-13"></div>
'''[13]''' Adi Levin. (1999) "Interpolating nets of curves by smooth subdivision surfaces". Computer Graphics Proceedings (SIGGRAPH)

<div id="cite-14"></div>
'''[14]''' Eugenio Oñate. (1995) "Cálculo de estructuras por el método de los elementos finitos". CIMNE