The generation of triangular meshes for NURBS‐enhanced FEM

This paper presents the first method that enables the fully automatic generation of triangular meshes suitable for the so‐called non‐uniform rational B‐spline (NURBS)‐enhanced finite element method (NEFEM). The meshes generated with the proposed approach account for the computer‐aided design boundary representation of the domain given by NURBS curves. The characteristic element size is completely independent of the geometric complexity and of the presence of very small geometric features. The proposed strategy allows to circumvent the time‐consuming process of de‐featuring complex geometric models before a finite element mesh suitable for the analysis can be produced. A generalisation of the original definition of a NEFEM element is also proposed, enabling to treat more complicated elements with an edge defined by several NURBS curves or more than one edge defined by different NURBS. Three examples of increasing difficulty demonstrate the applicability of the proposed approach and illustrate the advantages compared with those of traditional finite element mesh generators. Finally, a simulation of an electromagnetic scattering problem is considered to show the applicability of the generated meshes for finite element analysis. ©2016 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons Ltd.


Introduction
High-order discretisation methods have gained an increased popularity during the last decade due to the potential of providing higher accuracy with a reduced computational cost compared to traditional low-order methods [1,2,3,4,5].The advantages of high-order methods are particularly important in wave propagation problems, where low-order methods are known to suffer from excessive dissipation and/or dispersion when propagating waves over long distances [6,7,8,9,10,11].These methods have also attracted significant attention within the computational fluid dynamics community due the ability to propagate vortices over long distances [12,13,14,15].
The use of curved elements is crucial in order to fully exploit the benefits of high-order methods [16,17,18,19,20].This has prompted a great interest in the research community on the generation of high-order curvilinear meshes and nowadays there are some approaches to automatically generate such meshes [21,22,23,24,25,26].To fully exploit the potential advantages of high-order methods, very coarse curvilinear meshes and very high-order approximations are preferred, but as the size of the elements increases, the effect of geometric inaccuracies induced by the traditional polynomial approximation of curved boundaries inherent to the isoparametric formulation becomes evident.The error induced by the isoparametric formulation can be up to one order of magnitude higher than a formulation with an exact boundary representation of the computational domain [27].In addition, complex geometries of large scale objects often contain very small geometric features (e.g.holes and fillets) making necessary to produce extremely refined meshes in some regions of the domain in order to properly represent the geometry under consideration.In many occasions, it is necessary to invest a non-negligible amount of time removing these small features present in the computer aided design (CAD) model before attempting to produce a mesh suitable for finite element analysis [28].Although the feature removal can be assisted by some semi-automatic tools [29,30,31,32], the problem is not just the huge amount of time that is invested in removing geometric features but the uncertainty that the geometric simplification can generate in the finite element solution obtained in the simplified model.The problem becomes particularly dramatic when different simulations are required in the same geometry.For instance, a geometric feature that is not relevant in a fluid mechanics problem could be extremely relevant when an acoustic problem is solved.Furthermore, features not relevant in electromagnetic simulations at low frequency can become extremely influential at higher frequencies.
Although CAD and numerical simulation are ubiquitous in modern product development, these two technologies are still far from being integrated in a seamless manner.Automatic methods have been developed in the last decade that are capable of ensuring water tightness with minimum user intervention [33].However, when dealing with complex geometries that contain multi-scale features, it is often necessary to manually remove small geometric features that induce excessive and unnecessary mesh refinement.NURBSenhanced finite element method (NEFEM) was proposed to bridge the gap between CAD and finite element analysis [34].The main idea is to define the curved elements in contact with the CAD boundary in terms of the exact boundary description and not using the traditional isoparametric approach.The advantages of NEFEM compared to other curved finite element methods have been studied in detail from a theoretical and practical points of view both in two and three dimensions [35,27,36,37].
Despite all the benefits reported in the literature, the widespread application of NEFEM has been hampered by the lack of an automatic mesh generator able to generate the meshes that will allow to fully exploit its potential.In particular, a mesh generator able to generate elements which size is independent on the complexity of the boundary and that can be used for finite element analysis is not available.In fact, this lack of such an automatic mesh generator for NEFEM has motivated the development of methods that do not require the generation of fitted NEFEM meshes but still use the NE-FEM rationale [38,39,40].
This paper proposes a novel mesh generation technique that allows to produce triangular meshes where the elements account for the exact boundary representation of the domain irrespectively of the desired element size and the geometrical complexity.In particular the produced meshes contain elements with edges defined by more than one curve, avoiding small elements when very small NURBS curves are used to represent the boundary of the domain.Furthermore, the meshes generated with this technique can produce elements where an edge contains corners of the boundary representation of the domain, enabling to encapsulate small and complex geometric features within coarse triangular elements.The produced meshes allow to extend the NEFEM formalism introduced in [34], where it was assumed that an edge of an element must be defined by at most one curve.This paper also proposes a technique to extend the meshes to higher order by using a solid mechanics analogy only for the elements in contact with the NURBS bound-ary.Several examples demonstrate not only the applicability and potential of the proposed mesh generation technique but also demonstrate the applicability of the produced meshes by presenting a simulation of the scattering of electromagnetic waves by complex geometric objects using NEFEM.
The outline of the paper is as follows.In Section 2 the NEFEM formulation is recalled and the concept of a NEFEM element is extended with respect to its original definition.Section 3 summarises the mesh requirements and presents the proposed technique to generate linear meshes suitable for NEFEM.The generation of high-order meshes for NEFEM in presented in Section 4. Three numerical examples of increasing difficulty are used in Section 5 to illustrate the potential of the proposed technique and the proposed meshes are used to perform a simulation with NEFEM.Finally, Section 6 summarises the main conclusions of the work that has been presented.

NURBS-Enhanced FEM
This section introduces the fundamental concepts of NEFEM in two dimensional domains and generalises the original definition of NEFEM elements [37].

NURBS curves
A qth-degree non-uniform rational B-spline (NURBS) curve is a piecewise rational function defined in parametric form as where {B i } are the coordinates of the n cp + 1 control points (forming the control polygon), {ν i } are the control weights, and {C q i (λ)} are the normalized B-spline basis functions of degree q, which are defined recursively by for k = 1 . . .q, where λ i , for i = 0, . . ., n k , are the knots or breakpoints, which are assumed ordered 0 ≤ λ i ≤ λ i+1 ≤ 1.They form the so-called knot vector which uniquely describes the B-spline basis functions.The multiplicity of a knot, when it is larger than one, determines the decrease in the number of continuous derivatives.The number of control points, n cp +1, and knots, n k +1, are related to the degree of the parametrization, q, by the relation

NEFEM triangular elements
An open bounded domain Ω ∈ R 2 is considered, where its boundary ∂Ω is described using NURBS curves C k for k = 1, . . ., M , with M the total number of boundary curves A regular partition of the domain Ω = e Ω e in triangles is assumed, such that Ω i Ω j = ∅, for i = j.

Original definition of NEFEM triangular elements
The original definition of a NEFEM triangular element [34] assumes that • curved elements have, at most, one edge on the boundary, • every curved edge belongs to a unique NURBS curve and • internal edges are straight.
If Γ represents the edge of the element Ω e on the NURBS boundary parametrized by C, and x 1 , x 2 ∈ ∂Ω the two vertices on the NURBS boundary, see Figure 1, the curved edge is defined by where C(λ 1 ) = x 1 and C(λ 2 ) = x 2 .Then, a curved triangular element is defined as a convex linear combination of the curved edge and the interior node as as illustrated in Figure 1.Remark 1 (Visibility condition).The parametrisation (1) used to define curved elements with one edge defined by one NURBS curve implicitly assumes that the straight segments connecting the interior node x 3 with the points in Γ are contained in the domain Ω.This means that where n denotes the outward normal vector to Γ.

Extending the definition of NEFEM triangular elements
In this work, the original definition of a NEFEM triangular element is extended to include elements with more than one boundary edge, elements with boundary edges defined by more than one NURBS curve and elements that do not fulfil the so-called visibility condition detailed in Remark 1.As it will be shown, this extension facilitates the possibility of generating meshes with a characteristic element size completely independent of the complexity of the NURBS boundary.To simplify the presentation and without loss of generality, it is still assumed that internal edges are straight.
An edge Γ on the boundary with vertices x 1 and x 2 is defined by where n c denotes the total number of curves defining the edge Γ, The boundary of the element ∂Ω e is assumed to be a closed piecewise curve formed by trimmed NURBS curves and straight segments.
Formally, an element Ω e can be defined using the Jordan's curve theorem [42] as the subdomain where ind ∂Ωe is the index or winding number of a point x with respect to the closed curve ∂Ω e formed by the edges of the element.Figure 2 shows a triangular element that satisfies the new definition of a NEFEM element but does not satisfy the original definition.It can be observed that the triangle has two edges defined by NURBS curves and one edge is defined by two NURBS curves.In addition, the element does not satisfy the visibility condition in Remark 1.

NEFEM rationale
A distinctive feature of NEFEM is the definition of the polynomial approximation in the physical space, with Cartesian coordinates, rather than in the reference element with local coordinates.This choice ensures the reproducibility of polynomials in the physical space and has been shown to be advantageous compared to the traditional approximation with local coordinates [27].For simplicity, Lagrange nodal basis are considered to define the polynomial approximation although other basis (e.g.hierarchical) can be considered.
NEFEM performs the computation of the integrals appearing in a weak variational formulation using the exact boundary description of the computational domain.This section describes the strategy to compute the boundary and element integrals for the new definition of NEFEM elements introduced in Section 2.2.2.

Preprint of
Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247

Boundary integrals
A boundary integral to be computed along an edge on the boundary given by (3) can be written as where f is a generic function (usually a polynomial).It is worth noting that the integral has been split as a summation of integrals over different NURBS curves as the new definition of NEFEM curved elements allow the possibility to have one edge defined by several NURBS.A detailed comparison of different numerical integration techniques to evaluate the integrals of polynomial functions along NURBS curves [43] showed that one dimensional Gaussian quadratures defined over the parametric space of the NURBS provide the most efficient option.Therefore, the boundary integrals are approximated as where λ j i and ω j i are the coordinates and weights of the n j ip Gaussian integration points in [λ j 1 , λ j 2 ].It is worth recalling that a NURBS parametrisation is a piecewise rational function whose definition changes at the breakpoints so it is mandatory to define composite quadratures that account for the discontinuous nature of the parametrisation [43].

Interior integrals
The proposed strategy to compute interior integrals according to the new definition of a NEFEM element is to build a composite two dimensional quadrature in elements affected by the NURBS boundary representation of the domain.
Elements satisfying the hypothesis stated in Section 2.2.1 and the visibility condition detailed in Remark 1 are parametrised using the mapping given in Equation (1), as usual in a NEFEM context.An element Ω e not satisfying any of the previous hypothesis, see an example in Figure 2, can be split into K cells with only one face described by a NURBS curve and L straight-sided cells, namely where Ω c k denotes a cell with only one edge on defined by one NURBS curve and satisfying the condition (2) and Ω s l denotes a straight-sided cell.The partition satisfies that Ω c i Ω c j = ∅ and Ω s i Ω s j = ∅ for i = j and Ω c i Ω s j = ∅ for all i and j.
Figure 3 shows a partition of the NEFEM element shown in Figure 2 with three integration cells that have only one edge defined by one NURBS curve.An interior integral in an element Ω e that has been partitioned according to (4) is computed as

Integrals on cells Ω c
k , with at most one edge described by a NURBS curve and satisfying the visibility condition, are computed using the mapping given in Equation (1) whereas integrals in straight-sided cells are computed using standard triangle quadratures [44,45,46,47].
The details of the strategy to perform the partition of elements into integration cells with at most one edge described by a NURBS curve and satisfying the visibility condition are presented in Section 3.5.It is important to emphasise that the subdivision proposed to compute interior integrals is aimed at simplifying the implementation (i.e., to avoid the definition of a different mapping for each type of element).It is also worth remarking that this subdivision does not introduce new degrees of freedom, as the cells are only used to perform the numerical integration.

Generation of low-order NEFEM meshes
A mesh is usually generated in a hierarchical manner.For a two dimensional domain, vertices are meshed with nodes, curves forming the boundary are discretised with edges and, finally, the domain is meshed using elements.
The technique presented here follows a different rationale and allows to discretise the boundary with a required element size, independently of small geometric features.Finally, a variant of the advancing front method is proposed whereby mesh fronts coincident with the NURBS boundary match the exact geometric definition from CAD independently on the element size.

Mesh requirements
The following requirements for a NEFEM element are considered: 1. the characteristic element size is not restricted by small geometric features 2. element edges on the boundary are exactly described by the CAD boundary representation of the domain using NURBS 3. the boundary of a NEFEM element is a simple closed or Jordan curve (i.e.there is an injective map from a circle to the boundary of the element) 4. element edges interior to the domain are straight segments The first two requirements constitute the basis of NEFEM and ensure that de-featuring of complex geometries is not required before a finite element mesh suitable for a NEFEM analysis can be generated.The third requirement ensures that there are no self-intersections between the edges of an element.Finally, the last requirement is optional and is only aimed at improving the performance of the solver.For instance, a discontinuous Galerkin algorithm for solving Maxwell's equations in the time domain can be accelerated 10 times for low-order approximations and up to 100 times for high-order approximation by considering straight-sided triangles compared to curved triangles [48].

Boundary sampling points
The first step consists on defining a set of sampling points for each NURBS curve on the boundary.The approach considered uses a distribution function [49] defined through a background mesh and a set of points and line sources.The sampling points are generated at a distance α times the minimum value of the spacing distribution function, where α is a user defined parameter that by default is taken as 0.1.In practice, the user can also specify the regions where an exact boundary representation is of interest and the set of sampling points will only be built on those boundaries.This avoid the generation of a large set of sampling points in problems such as external flow computations where the far field shape and geometric representation is not relevant.
The initial set of sampling points induces a polygonal description of curved boundaries.Extra sampling points are added in regions where the distance from the straight sided segment connecting two sampling points and the true CAD boundary is more than a specified tolerance.In the proposed implementation, extra points are added when the distance is larger than 10% of the minimum spacing.Next, the desired mesh spacing at the sampling points is computed from the spacing distribution function.Finally, the spacing of the sampling points is checked and corrected if it induces an element on the boundary that is ε times smaller than the desired element size.Typically, a value of ε ∈ [0.5, 0.75] is considered in the implementation of the proposed algorithm.
It is worth noting that the parameter α should be chosen to ensure that all geometric features are captured.Due to the marginal cost associated to the generation of sampling points, the use of small values of α is advocated here, but other techniques based on the generation of sampling points over the Bézier curves forming the NURBS with curvature control could be devised.

Boundary discretisation
The first mesh requirement described in Section 3.1 implies that the traditional, hierarchical, meshing paradigm cannot be adopted.It is clear that using a hierarchic approach, where the first step consists on adopting the vertices of the CAD model as mesh nodes, implies that small geometric features might dictate the element size.As an example, let us consider a domain Ω = [−25, 25] 2 \[−5, 5] × [−0.25, 0.25], representing a square plate with a rectangular inclusion of dimension 10 × 0.5 whose centre of gravity is located in the origin, as represented in Figure 4(a).If the desired element size is, for Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247instance, 3.5 and the traditional hierarchical approach is utilised, the minimum element size will be dictated by the thickness of the inclusion (i.e., 0.5) as illustrated in Figure 4(b).The proposed approach starts by combining the boundary curves into loops.A loop is defined as a collection of boundary curves L := {C k } k=1,...,N L , where N L is the total number of curves forming the loop, The curves forming the loop are ordered such that the first curve, and the total length of the loop is denoted by ).The proposed procedure to discretise boundary loops with edges of a desired size is detailed in Algorithm 1.

Preprint of
Compute the candidate mesh node, if not isValid then Get the sampling points within x i−1 and x i , λ l i−1 l=1,...,ns , and the associated index of the curves where the points belong, s l l=1,...,ns A candidate boundary vertex x i at a distance h(x i−1 ) of a given vertex x i−1 is first identified.As the proposed methodology enables the creation of element edges across different boundary curves, it is necessary to identify the curve where the candidate vertex belongs, C s , and the parametric coordinate of the candidate vertex, λ i .The strategy to identify both the curve C s and the parametric coordinate λ i is detailed in Algorithm 2. It is worth noting that Algorithm 2 requires the solution of a one-dimensional non-linear problem.In practice this is achieved using a simple bisection approach as the curves are initially sampled and convergence is achieved in very few iterations.
Next, the mesh front formed by the given vertex x i−1 and the candidate vertex The objective is to guarantee that, in the second stage of the proposed meshing technique described in Section 3.4 (i.e., domain discretisation), an element satisfying the third mesh requirement stated in Section 3.1 exists.The proposed strategy to check the validity of a mesh front is summarised in Algorithm 3.
For a given edge Γ, the horizon of an edge vertex x k of Γ is defined as the most distant point x k ∈ Γ that can be connected to x k without intersecting , for all x ∈ Γ and x k x k ⊂ Ω, where d denotes the distance function and xy is the straight segment connecting x and y end for 2. Check if a third node can be found to form a valid element Find P as the intersection between the line connecting x i−1 and x i−1 and the line connecting Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Γ, as illustrated in Figure 5.The intersection of the lines connecting the boundary vertices and the corresponding horizons, denoted by P in Figure 5, is used to decide the validity of the mesh front.If the distance from P to the boundary vertices is less than the desired spacing h this implies that an interior node to form a triangle with the required spacing can be found.An example of a non-valid mesh front is shown in Figure 6.As detailed in Algorithm 1, if the mesh front is valid, the second boundary vertex x i is accepted.Otherwise, a more restrictive criteria is applied in order to find the boundary vertex x i that will guarantee a valid mesh front.The procedure, summarised in Algorithm 4, uses the set of sampling points defined in Section 3.2 to find a boundary edge of maximum length that forms a valid mesh front.Compute the candidate mesh node, Starting with the position of x i that produces a non-valid front (marked with a red dot in Figure 7), the sampling points, x l i−1 l=1,...,ns , are checked sequentially starting with the sampling point closest to the discarded vertex x i .For each sampling point, the proposed algorithm checks the validity of the mesh front formed by x i−1 and the current sampling point until a valid front is found.Figure 7 shows the sampling points that produced invalid mesh front with red crosses.The position of x i is adjusted to be the position of the first sampling point producing a valid front.
Remark 2. The check for the validity of a mesh front proposed in Algorithm 3 is only relevant when the produced mesh has to verify the fourth mesh requirement described in Section 3.1.When curved internal edges are allowed, it is possible to eliminate this step and build the boundary discretisation only based on the desired spacing.
Remark 3. In the implementation of the proposed method, the user can specify parts of the boundary where the curves will be grouped in loops and other parts of the domain where the standard, hierarchical, approach will be utilised.For instance, in the domain shown in Figure 4, only the curves defining the inclusion can be selected to form an inner loop and maintain the curves defining the outer part of the domain to be treated in a hierarchical manner.This is of particular importance when problems in an exterior domain are of interest (e.g., external flows, electromagnetic scattering) as the Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247 Figure 7: Adjustment of the vertex x i using the sampling points, denoted with a cross, as detailed in Algorithm 4 to produce a valid front of maximum length after the non-valid mesh front shown in Figure 6 is found.
so-called far field boundary is introduced for computational purposes and its shape and geometric representation is irrelevant.

Domain discretisation
The advancing front method (AFM) is one of the most popular unstructured mesh generation procedures [50].The main idea is to use the concept of a generation front, initially a sequence of segments that connects the mesh nodes obtained during the discretisation of the boundary.The procedure used to build a triangular element consists on selecting a segment from the front and looking for the ideal point to form the element.The method is designed to try to form equilateral triangles whenever is possible and it is usually finalised by a mesh quality enhancement procedure [50].
The standard AFM is not applicable to generate meshes suitable for NE-FEM because boundary edges are considered to be straight segments.Therefore, if the AFM is applied to the discretisation of the boundary obtained after following the strategy described in Section 3.3, the third mesh requirement stated in Section 3.1 is only guaranteed if the straight segment is inside the domain Ω.This is illustrated in Figure 8.A zoom of the domain represented in Figure 4 is shown in Figure 8(a).the elements generated with the standard AFM might be invalid if the exact boundary representation is considered.
This section proposes a variation of the standard AFM that accounts for the exact boundary representation of the domain and meets the mesh requirements stated in Section 3.1.The procedure is described in Algorithm 5.
The first step consists on computing the horizon of the boundary vertices, x 1 and x 2 , as proposed in Algorithm 3. The intersection of the lines connecting the boundary vertices and their respective horizons is denoted by P as shown in Figure 9.A line of search is defined as the angle bisector of the two lines connecting the boundary vertices and their respective horizons, denoted by a dashed red line in Figure 9.The third vertex of the element is defined as where t takes the maximum value that guarantees that the edges x 1 x 3 and x 2 x 3 do not intersect the boundary and have length not greater than the desired element size h.
, for all x ∈ Γ and x k x k ⊂ Ω, where d denotes the distance function and xy is the straight segment connecting x and y end for 2. Compute the third node using the modified AFM Find P as the intersection between the line connecting x 1 and x 1 and the line connecting x 2 and x 2 Define n P = (x 1 − x 1 ) + (x 2 − x 2 ) Find x 3 ∈ Ω along the line defined by P and n P such that d(x i , x 3 ) ≤ h and x i x 3 ∈ Ω for i = 1, 2. In some situations, using the bisector to obtain the position of the third node induces the generation of small internal edges, as illustrated in Figure 10(a).In such cases, it is advantageous to introduce a rotation of the line where the third vertex will be positioned.The following expression for the vertex x 3 is used when the expression in Equation ( 5) produces an internal edge much smaller than the desired spacing and the angle α is selected to minimise the difference between the length of both internal edges and the desired spacing.Figure 10 illustrates the benefit of using this definition for the position of the third vertex.In this example the definition in Equation ( 5) produces an internal edge of size 0.6h, being h the desired element size, whereas the definition in Equation ( 6) produces an internal edge of size 0.9h.In this figure, the shaded area in blue shows the region of possible positions for the third vertex of the element.
It is worth noting that the proposed modification of the AFM described in Algorithm 5 reduces to the standard AFM if the mesh front is a straight-sided segment.

Numerical integration cells
As mentioned earlier, in Section 2.3.2, the strategy to perform the interior integrals appearing in the weak formulation is to use the convex linear 21 Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247 Figure 10: Illustration of the modified AFM described in Algorithm 5 to build a triangular element with the definition of the internal point using (a) Equation ( 5) and (b) Equation (6).
mapping of Equation ( 1) [34].The new definition of curved elements introduced in this work requires a partition of the elements with more than one edge on the NURBS boundary or with an edge defined by more than one NURBS curve.This section presents the proposed strategy to build the integration cells for those elements.It is worth emphasising that the cells are only used to perform the numerical integration, so no new degrees of freedom are introduced.
A recursive technique is proposed in Algorithm 6 to build integration cells that have, at most, one edge defined by a NURBS curve.To simplify the notation, it is assumed that when the element has two edges on the boundary, the internal edge connects the vertices x 1 and x 3 , otherwise the node not on the boundary is x 3 .
The proposed algorithm uses a point interior to the domain, Q to check the visibility of the sampling boundary points.The point Q is selected to be the third node if the element has only one edge defined by NURBS as illustrated in the example of Figure 11.Otherwise, if the element has two Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247

Algorithm 6 integrationCellsOfAnElement
Compute the integration cells for a boundary element Input: Vertices forming the element x 1 , x 2 and x 3 .Input: NURBS curves, {C i } i=1,...,nc , describing the boundary edge/s Γ Input: Sampling points x l i−1 l=1,...,ns on the NURBS boundary Output: edges on the NURBS boundary, the point Q is selected to be the mid point of the only internal edge.The visibility of the sampling points is checked starting with the first sampling point x l 1 = x 1 until a sampling point that is not visible is found.In the example of Figure 11(a), the first non-visible sampling point is x l 3 .The first integration cell, in blue in Figure 11(b), is formed using the last visible Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247point, x l 2 .The algorithm flags the last visible point as x1 and continues until a new visible sampling point is found, x2 in Figure 11(b).This allows to form a second integration cell that has no edges on the NURBS boundary, the cell in red in Figure 11(b).The definition of the point Q is now changed to be the mid-point between x1 and x2 and the algorithm is applied to the triangle formed by x1 , x2 and Q, producing two new cells as illustrated in Figure 11(c).The final partition in integration cells with at most one edge defined by a NURBS curve is shown in Figure 11(d).

Extension to higher order elements
This section proposes a procedure to build high-order meshes from the low-order meshes produced using the method described in Section 3. The proposed strategy is based on the linear elastic analogy [23] but it is simpler because an element-by-element solution is proposed in order to maintain internal edges as straight-sided segments, following the mesh requirements introduced in Section 3.1.

High-order boundary nodal distribution
An edge Γ on the boundary with vertices x 1 and x p+1 , defined in (3), is considered where C 1 (λ 1 1 ) = x 1 and C nc (λ nc 2 ) = x p+1 .Algorithm 7 details the proposed procedure to build a high-order nodal distribution on Γ.
First, the desired distance between high-order boundary nodes is computed.This distance can vary from node to node and depends upon the desired nodal distribution that is specified by the user, namely {ξ k } k=1,...,p+1 ∈ [0, 1].Equally-spaced nodal distributions in [0, 1] can be easily constructed although other alternatives, with better approximation properties, are available [51].The second step follows the same rationale of the procedure described in Algorithm 1, in Section 3.3, to discretise boundary curves forming loops, but changes a desired element size h by a desired distance between boundary nodes l k .
For interior edges, considered as straight-sided segments, a high-order nodal distribution of the desired degree p is placed by mapping the nodal distribution {ξ k } k=1,...,p+1 ∈ [0, 1] to the edge using the linear mapping uniquely defined by the vertices of the edge.

High-order elemental nodal distribution
For an element Ω e with closed boundary ∂Ω e formed by three edges {Γ i } i=1,2,3 and with at least one edge on the NURBS boundary, the triangle formed by the vertices of Ω e is denoted by T e .A high-order nodal distribution of the desired degree p is placed in T e .This can be easily performed by mapping the nodal distribution from a reference element to T e using the linear mapping uniquely defined by the vertices of T e .
Then, the following linear elastic problem is defined The constitutive relation relates the stress tensor σ to the deformation tensor ε introducing the Young's modulus E and Poisson's ratio ν for the elastic medium [52,53].In the above expression, I is the identity tensor, tr denotes the trace operator and the deformation tensor is defined in terms of the displacement field u as The problem is closed with the following discrete Dirichlet boundary condition where x 0 k denotes the position of a high-order edge node in T e and x k denotes the position of an edge node obtained from Algorithm 7 described in the previous section.It is worth noting that this represents an homogeneous Dirichlet boundary condition if x k belongs to an interior edge.
The procedure is completed by solving the linear elastic problem (7).The triangle T e represents the initial, undeformed, configuration.After solving the elastic problem, a high-order nodal distribution is obtained in Ω e .It is worth remarking that due to the nature of the imposed boundary conditions and the solution of the problem element-by-element, if T e satisfies the mesh requirements specified in Section 3.1, then the element Ω e produced with this strategy also satisfies the same mesh requirements.
Figure 12 illustrates the proposed procedure to build a high-order nodal distribution in a NEFEM element.Figure 12(a) shows a fourth-order nodal distribution with equally-spaced nodes in T e .The boundary nodal distribution obtained from Algorithm 7 is used to define the Dirichlet boundary condition illustrated in Figure 12(b) with arrows.The resulting nodal distribution in Ω e is also represented in Figure 12(b).
Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Remark 4. As mentioned earlier in Section 2.3, NEFEM defines the approximation directly in the physical space, with Cartesian coordinates.Lagrangian polynomials are directly built in Ω e based upon the position of the high-order nodes obtained with the elastic analogy [34,27].This gives complete freedom on the position of the nodes in Ω e , contrary to the isoparametric FE formulation, where the position of the nodes in the physical element must be carefully selected in order to guarantee the optimality of the approximation [54,55].It is therefore possible to introduce source terms in the elastic problem (7) to obtain the high-order nodes in Ω e following a particular spatial distribution.In fact, it is even possible to specify the position of the interior nodes without solving an elastic problem.This could represent an advantage if the solution is known to be complex in certain regions as the nodes can be clustered to better represent the physics of the problem under consideration.
Remark 5.The strategy proposed in this section can be adapted to produce meshes conforming to the CAD boundary representation of the domain and with interior curved edges.This can be achieved by solving the elastic problem in the whole domain Ω instead of performing an element-by-element solution.
In fact, a more efficient layer-by-layer approach can be implemented [23].The first step involves the solution of the elastic problem in a layer of elements surrounding the boundary represented by NURBS.Dirichlet boundary conditions are used for the edges on the boundary and homogeneous Neumann boundary conditions are used in the outer boundary of the layer.The second step consists on solving another elastic problem in a second layer of elements surrounding the first layer.Dirichlet boundary conditions correspond to the displacement field from the outer boundary of the first layer.The procedure continues until the last layer of elements is deformed or until the displacement field of one layer is small enough.

Examples
This section presents three examples that illustrate the applicability of the proposed approach for the generation of meshes suitable for NEFEM.The meshes are analysed in terms of the desired element size and the actual minimum element size of the produced meshes.An application example is also considered where the produced mesh of a complex object is utilised to perform an electromagnetic scattering simulation with NEFEM.
In the implementation of the proposed methodology the input format of the CAD geometry is considered an IGES.Other input formats could Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247be easily considered because the implementation only uses the control points and associated weights and the knot vector for each NURBS boundary entity.

Low-order mesh around an aerofoil with a blunt trailing edge
The first example involves the generation of a NEFEM mesh around an aerofoil of unit chord length with a blunt trailing edge of dimension 0.001.Figure 13 shows a NEFEM mesh generated with the technique proposed in this paper.The colour field represents, in logarithmic scale, the spacing dis- tribution function that is induced by two point sources placed at the leading and trailing edges and several line sources to control the element size on the upper and on the lower parts of the aerofoil.The desired element size is 0.005 on the leading edge, 0.01 on the trailing edge and 0.025 on the upper and lower parts of the aerofoil.
The geometry of the aerofoil is described using three NURBS curves and the outer (far-field) boundary is treated in the usual hierarchical manner.The total number of elements with an edge and a node on the NURBS boundary is 89 and 92 respectively.Figure 14 shows two histograms, summarising the size of the edges on the NURBS boundary and the edges of elements with a node on the NURBS boundary.The length is normalised dividing by the desired spacing that, for each edge, is computed as the average of the desired spacing at the nodes.
The minimum and maximum generated edge length on the boundary is approximately 0.00532 and 0.0253 respectively, illustrating the ability of the proposed approach to discretise the boundary with a desired element size, irrespectively of small geometric features.In addition, the minimum and maximum edge length for an element with a node on the boundary is 0. and 0.0261 respectively, showing that the required spacing is maintained for those elements in contact with the boundary elements, generated using the modification of the advancing front method.It is worth noting that all the other elements (i.e., elements not in contact with the NURBS boundary) are generated using the standard advancing front method and, therefore, the edge length for those elements is not reported here.
Figure 15 shows a detailed view near the blunt trailing edge for two meshes, generated with a standard mesh generator and with the proposed technique.Figure 15(a) shows that using a standard mesh generator the minimum element size is given by the dimension of the blunt trailing edge, 0.001 in this example.The mesh produced with the proposed technique, shown in Figure 15(b), maintains the desired element spacing even in regions where small geometric features are present.In this case, the size of the elements near the trailing edge is five times bigger than the dimension of the blunt trailing edge.
The main advantage of the proposed meshing technique in this example is not to reduce the number of elements compared to a traditional FE mesh, but to avoid the small elements introduced by a standard FE mesh generator to capture the blunt trailing edge.The standard FE mesh shown in Figure 15 has a minimum element size five times smaller than the minimum element size of the NEFEM mesh.This difference in element size might have an important impact in the total CPU time of a simulation when an explicit time marching algorithm is employed.In addition, the proposed technique avoids abrupt transitions of the element size, inducing better approximation properties, specially if the mesh is used for a low order computation.

Low and high-order meshes around an aircraft profile
The second example involves the generation of a NEFEM mesh around the aircraft profile represented in Figure 16.A detailed view of the front part is also shown, revealing a number of very small features compared to the size of the aircraft.The aircraft profile is described by 53 NURBS curves and the minimum and maximum length of these NURBS curves is 2×10 −6 and 2.5 respectively, differing in more than six orders of magnitude.A mesh generated with the proposed approach and with a uniform desired element size of 0.04 (i.e., 2000 times larger than the smallest geometric feature) is shown in Figure 17.The mesh has 121 elements with one or more edges on the boundary described by NURBS and 119 elements with one node on the NURBS boundary.Figure 18 shows a detailed view of the mesh near the most critical zone containing the NURBS of length 2 × 10 −6 and the integration cells in three of the elements with an edge defined by four or even five NURBS curves.It is worth noting that the desired element size is larger than the length of 33 of the boundary curves.
The minimum and maximum generated edge length on the NURBS boundary are 0.0319 and 0.0497 respectively, clearly demonstrating the ability to generate meshes of the desired element size even when geometric features of significant small size are present.For the elements with one node on the NURBS boundary, the minimum and maximum generated edge length are 0.0300 and 0.0689 respectively.Figure 19 shows a histogram summarising the normalised length of the edges of elements with at least one node on the NURBS boundary.
To further illustrate the benefits of the proposed meshing technique, Figure 20 shows a detailed view of a mesh produced using a standard finite element mesh generator with the same desired element size.The resulting mesh has 311 edges on the NURBS boundary but, more importantly, the  minimum element size is equal to the minimum size of the geometric features in the domain, 2 × 10 −6 in this example.The mesh generator used here is the two dimensional version of the FLITE system, a mature in-house mesh generator.It worth mentioning that point sources have been added near the regions with very small geometric features to guarantee that the produced mesh does not encompass abrupt changes of the element size and, therefore, it is a suitable mesh for finite element analysis.
Next, the mesh shown in Figure 17 is extended to high-order by using the technique proposed in Section 4. For a degree of approximation p = 5, a reference element with the desired nodal distribution is considered.In this example, a Fekette nodal distribution [56] is employed instead of the more traditional equally-spaced nodal distribution, due to the well known superior approximation properties [51,57].The nodal distribution is mapped to each element in the mesh using the affine mapping that is uniquely defined by the vertices of each element.The result is a nodal distribution not conforming to the boundary as illustrated in Figure 21.Using the solid mechanics analogy on each element with at least one edge on the NURBS boundary, the nodal distribution is adapted to the exact geometry and the resulting high-order NEFEM mesh is obtained as represented in Figure 22.
When adapting the nodal distribution to the true geometry using the solid mechanics analogy, it is necessary to specify the material parameters (i.e., Poisson's ratio and Young modulus) for the elastic medium.When the elastic analogy is used to generate high-order curvilinear meshes for isoparametric finite elements [23] the effect of the materials parameters on the quality of the produced meshes is relevant, but this is not the case for NEFEM.First, as NEFEM defines the approximation on the physical space, the scaled Jacobian, typically used by the FE community to evaluate the quality of a FE mesh, is irrelevant.Second, this work proposes the solution of an elementby-element elastic problem instead of a global problem, to guarantee the last mesh requirement in Section 3.1.The effect of the Poisson's ratio has been found to be minimum on the final nodal distribution obtained.As an example, Figure 23 shows two detailed views of the nodal distributions obtained with ν = 0 and ν = 0.49 for an approximation with degree p = 3.Despite the large deformation of the elements shown, the variation of the position of the inner point is very small, usually more than 10 times smaller than the element size.For higher degrees of approximation this difference becomes even less noticeable.
As expected, the value of the Young modulus has no effect on the solution of the elastic problem because imposed displacements are considered on the whole boundary.

High-order mesh around a satellite profile
The last example considers a more complex geometry corresponding to a profile of a satellite as shown in Figure 24.In this example there are multiple regions with very small geometric features and an attempt to remove geometric features not relevant for the finite element analysis stage would clearly result in an extremely tedious and time consuming process.The satellite profile is described by 139 NURBS curves and the minimum and maximum length of the NURBS curves are 0.01 and 2.4 respectively.

Ω
Figure 24: Domain Ω representing the exterior of a complex satellite profile with many small geometric features.
A mesh generated with the proposed approach and with a desired element size of 0.06 is shown in Figure 25.It is important to notice that, in this example, almost 70% of the curves have a length smaller than the desired element size.The mesh has 2,519 elements, 1,410 nodes and 143 edges on the boundary described by NURBS.The potential of the proposed approach is clearly illustrated by comparing the number of boundary edges and the total number of NURBS describing the boundary.The similarity between these numbers indicates that the majority of the elements will contain edges defined by more than one NURBS curve.
Figure 26 shows a detailed view of the mesh near two zones containing small geometric features.The nodal distribution corresponds to a polynomial approximation with p = 4 generated using the solid mechanics analogy presented in Section 4.
The minimum and maximum generated edge length on the NURBS boundary are 0.0393 and 0.0837 respectively whereas the minimum and maximum generated edge length for elements with one node on the NURBS boundary is 0.0296 and 0.1045 respectively.Figure 27 shows a histogram summarising the size of the edges on the NURBS boundary normalised with the desired spacing.
Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Next, the effect of the parameter α, introduced in Section 3.2 to produce the distribution of sampling points, is studied.Table 1 reports the total number of elements (n el ), the number of elements with at least one edge on the NURBS boundary (n e el ), the total number of nodes (n no ) and the minimum and maximum edge size (h min and h max respectively), for the meshes generated using different values of the sampling parameter α.The results clearly show that the proposed methodology is not sensitive to this parameter, suggesting that a more sophisticated and computationally expensive approach to produce the distribution of sampling points is not required.
As described in Section 3, the overhead of the proposed approach and a standard mesh generator based on the advancing front method is restricted to the generation of elements with at least one edge defined by NURBS curves.In this example, the time spent on generating the elements with at least one edge on the boundary is 8% of the total mesh generation time.On average, generating a NEFEM element is 2.5 times more expensive than generating a standard triangular element using the advancing front method.It is worth remarking that this time includes the generation of the integration cells as described in Section 3.5.
Finally, to illustrate the applicability of the proposed technique, an electromagnetic scattering problem [58,8] is solved using the mesh shown in Figure 25.The problem is governed by the transient Maxwell's equations and consists of an electromagnetic wave travelling from the left to the right and scattered by the satellite profile.The frequency of the wave is such that the satellite, assumed to be a perfect electric conductor, occupies 14 wavelengths.The problem is solved using a high-order discontinuous Galerkin formulation [6] with NEFEM and a standard fourth order Runge-Kutta explicit time integrator.The three components of the scattered field are shown in Figure 28 for the transverse magnetic mode [58].A reference solution is computed on a standard FE mesh that contains 208,843 elements and 105,716 nodes.The element size has been selected to ensure that at least 20 elements per wavelength are present when a degree of approximation p = 1 is employed.For this level of mesh refinement, all the geometric details are captured by a standard mesh generator without any extra refinement.Therefore, this example shows the potential of the proposed technique in a high-order context.
A comparison of the radar cross section (RCS), an engineering quantity of interest in electromagnetic scattering problems [58], is shown in NEFEM and the mesh shown in Figure 17 is less than 2 × 10 −2 in the L 2 norm, demonstrating the possibility of computing accurate solutions with extremely coarse meshes generated with the proposed approach.
In this example, the required cost per cycle is reduced by a factor of 87 by using NEFEM with p=4 compared to the reference solution computed with linear elements in a finer mesh.In addition, the linear computation required 79 cycles in order to converge the solution to a harmonic steady state with a tolerance of 10 −6 in the RCS, whereas the NEFEM computation required 49 cycles.This means that the total cost reduction factor with NEFEM and the coarse mesh is approximately 140 (79×87/49); an improvement of more than two orders of magnitude.
It is worth remarking that the use of standard high-order isoparametric finite elements is not competitive in this example.Due to the small geometric features, the element size is always restricted by the length of the smallest curve.In fact, a high-order solution will introduce a dramatic restriction in the time step, resulting in a more costly solution compared to linear elements.Therefore, the speed up of the proposed NEFEM approach will be even higher than 140 if the comparison was performed with respect to high-order isoparametric finite elements.
Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247

Concluding remarks
The first fully automatic mesh generation technique for NEFEM has been presented.It allows to generate meshes that account for the CAD boundary representation of the domain given by NURBS curves and with an element size independent on the geometric complexity of the computational domain.
The proposed technique allows to generalise the original definition of a NEFEM element, ensuring that elements with more than one edge on the NURBS boundary or edges defined by more than one NURBS can be considered.A technique to discretise the boundary with a desired spacing and independent on small geometric features is introduced and a modification of the advancing front method is proposed.A method to extend the produced meshes to higher order approximation has also been presented.
Three examples demonstrate the applicability and potential of the proposed approach by generating meshes of complex objects with small geometric features.An analysis of the produced meshes has been presented in terms of the desired element size.Finally, a numerical example involving the simulation of the scattering of electromagnetic waves demonstrates the applicability of the proposed meshes in a NEFEM context.
It is in the scope of a future publication to present the extension of the methodology presented in this paper to three dimensional domains.

Figure 1 :
Figure 1: Parametrization of a curved triangular element with an edge defined by a NURBS curve.

Figure 2 :
Figure 2: A valid NEFEM element with the new definition that does not fulfil the original requirements stated in Section 2.2.1.

Figure 3 :
Figure 3: Partition of the NEFEM element shown in Figure 2 with three integration cells that have only one edge defined by one NURBS curve.

ΩFigure 4 :
Figure 4: (a) Domain Ω representing a square plate with a rectangular inclusion and (b) a standard FE mesh showing that the minimum element size is dictated by the thickness of the inclusion.

Algorithm 1
Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Discretisation of a closed loop formed by NURBS curves Input: Spacing distribution function h(x) Input: NURBS curves forming a closed loop L := {C k } k=1,...,N L Output: Mesh nodes {x i } i=1,...,M L discretising the loop L Set i = 1, λ i = 0, r = 1 and s = 1 Define the first node as Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Algorithm 2 boundaryNodeAtGivenDistance Find boundary vertex x i at a distance δ from another boundary point x i−1 Input: Loop L = {C k } k=1,...,N L Input: Index r of the curve to which the point x i−1 belongs Input: Parametric coordinate λ i−1 of the point x i−1 Input: Distance δ Output: Parametric coordinate λ i of node x i and index s of the curve to which the point x i belongs 1. Identify the curve to which the vertex

Figure 5 :
Figure 5: Illustration of the procedure described in Algorithm 3 to check the validity of a mesh front.
Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Algorithm 4 boundaryNodeAtMaxDistance Find the second vertex of a boundary edge Γ using sampling points Input: Loop L = {C k } k=1,...,N L Input: Boundary vertex x i−1 Input: Parametric coordinates of the sampling points λ l i−1 l=1,...,ns Input: Index of the curve to which each sampling point belongs s l l=1,...,ns Output: Parametric coordinate λ i of node x i and index s of the curve to which the point x i belongs Set isValid = false, l = n s , λ i = λ l i−1 and s = s l while isValid = false do

Figure 8 (ΩFigure 8 :
Figure 8: Generation of a coarse mesh around a rectangular inclusion with the standard AFM leading to an invalid NEFEM element: (a) domain Ω with an exact boundary representation, (b) altered domain Ω with an approximation boundary representation, (c) one element generated with the AFM in Ω, and (d) non-valid NEFEM element in Ω, showing a self-intersection of one internal edge and the exact boundary.

Figure 9 :
Figure 9: Illustration of the modified AFM described in Algorithm 5 to build a triangular element.

1
and visible=false for k = 1, . . ., n s do if Qx l k ⊂ Ω and visible= true then Set x2 = x l k−1 and visible=false Add cell with vertices x1 , x2 and Q to the list of integration cells x1 = x2 else if Qx l k ⊂ Ω and visible= false then Set x2 = x l k and visible=true Add cell with vertices x1 , x2 and Q to the list of integration cells xm = 0.5(x 1 + x2 ) subMeshOfAnElement x1 , x2 , xm , {C i } , Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247

Figure 11 :
Figure 11: Illustration of the Algorithm 6 to build the integration cells of an element.(a) Checking the visibility of sampling points, (b) formation of the first two cells, (c) recursive application of the algorithm and (d) final partition of integration cells.

Figure 12 :
Figure 12: Generation of a high-order nodal distribution in a NEFEM element: (a) domain Ω and high-order nodal distribution in the straight-sided triangle T e , (b) high-order nodal distribution adapted to Ω e obtained by solving a linear elastic problem in the element.

Figure 13 :
Figure 13: NEFEM mesh around an aerofoil with a blunt trailing edge and spacing distribution function in logarithmic scale.

Figure 14 :
Figure 14: Histogram of the normalised length of (a) the edges on the NURBS boundary and (b) the edges of elements with a node on the NURBS boundary.The histograms correspond to the mesh of Figure 13.

Figure 15 :
Figure 15: Detailed view of (a) a standard FEM mesh and (b) a NEFEM mesh near the blunt trailing edge of an aerofoil.

Ω
a) Domain Ω representing the exterior of an aircraft profile and (b) detailed view of the domain showing a number of small geometric features.

Figure 17 :
Figure 17: NEFEM mesh around a complex aircraft profile.

Figure 18 :
Figure 18: Details of the NEFEM mesh of Figure 17 near the regions containing the two shortest NURBS curves.The black dots denote the vertices of the generated triangular elements and the red lines are the subcells used for numerical integration.

Figure 19 :
Figure 19: Histogram of the normalised length of the element edges on the NURBS boundary corresponding to the mesh of Figure 17.

Figure 20 :
Figure 20: Detailed views of a standard FE mesh near the regions containing the two shortest NURBS curves.

Figure 21 :Figure 22 :Figure 23 :
Figure 21: Detail of the NEFEM mesh and the fifth order Fekette nodal distribution mapped from the reference element to the imaginary straight-sided element defined by the vertices of each element represented with discontinuous lines.

E 3 Figure 28 :
Figure 28: Components of the scattered field by a satellite profile computed with NEFEM and p = 4 on the mesh shown in Figure 17.

Figure 29 .Figure 29 :
Figure29: RCS computed with a reference FE mesh and with a NEFEM mesh and p=4.
Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Algorithm 3 checkValidityOfMeshFront Check the validity of a mesh front with nodes x i−1 and x i Input: Edge nodes of a mesh front x i−1 and x i Input: NURBS curves, {C i } i=1,...,nc , describing the boundary edge Γ Input: Desired element size h Output: Boolean isValid denoting the validity of mesh front with nodes x i−1 and x i 1. Compute the horizon of the boundary vertices for k Preprint of Ruben Sevilla, Luke Rees and Oubay Hassa, "The generation of triangular meshes for NURBS-enhanced FEM" International Journal for Numerical Methods in Engineering, 2016, DOI: 10.1002/nme.5247Algorithm 5 Build a NEFEM element from a mesh front with nodes x 1 and x 2 Input: Edge nodes of a mesh front x 1 and x 2 Input: NURBS curves, {C i } i=1,...,nc , describing the boundary edge Γ Input: Desired element size h Output: Third node x 3 to form a triangular NEFEM element Ω e 1. Compute the horizon of the boundary vertices for

Table 1 :
Effect of the parameter α used to build the distribution of boundary sampling points.