A critical view on the use of Non‐Uniform Rational B‐Splines to improve geometry representation in enriched finite element methods

Enriched finite element methods have gained traction in recent years for modeling problems with material interfaces and cracks. By means of enrichment functions that incorporate a priori behavior about the solution, these methods decouple the finite element (FE) discretization from the geometric configuration of such discontinuities. Taking advantage of this greater flexibility, recent studies have proposed the adoption of Non‐Uniform Rational B‐Splines (NURBS) to preserve the interfaces' exact geometries throughout the analysis. In this article, we investigate NURBS‐based geometries in the context of the Discontinuity‐Enriched Finite Element Method (DE‐FEM) based on linear field approximations. While optimal convergence is retained for problems with weak discontinuities without singularities, representing exact geometry via NURBS does not yield noticeable improvements when extracting stress intensity factors of cracked specimens. For low‐order elements, we conclude that the benefits of exact geometry representation do not outweigh the increased complexity in formulation and implementation. The choice of linear FEs hinders the accuracy of the proposed formulation, suggesting that its full potential may only be unleashed by increasing the field representation order.

is important, for instance, in problems with material interfaces, across which not only material properties can differ, but also the governing physics. In many cases it is precisely the field near discontinuities which is of main interest and which mandates for a more accurate geometric representation. In this work, we provide a critical evaluation of the use of NURBS for the enriched FEA of C 0 -and C −1 -continuous fields.
The perspective of improving shape representation and, for certain geometries, completely eliminating the geometric discretization error, has inspired researchers in computational mechanics to investigate the use of NURBS. In Iso-Geometric Analysis (IGA), 2,4 analysis is conducted on the exact spline-based (often NURBS-based) geometry that is generated by CAD software. The solution field is modeled by means of the same spline basis, following an isoparametric approach that can compute high-order fields with high accuracy. Drawbacks of IGA lie in the noninterpolatory nature of degrees of freedom (DOFs), which does not allow to enforce Dirichlet boundary conditions (BCs) exactly, and in the computational effort to process NURBS over the whole domain. In order to exploit the more than 60 years of standard finite element (FE) technology, some authors have explored combining NURBS with FEA. The NURBS-enhanced finite element method (NEFEM) 5,6 abandons the isoparametric duality between interpolant and geometry commonly used in both FEA and IGA, and adopts a spline-based geometrical representation in the Finite Element Method (FEM) while maintaining Cartesian shape functions, following an approach inspired by the blending function method. 3,7 A NURBS-enhanced geometric mapping is applied to boundary elements, while inner elements are treated as in standard FEM. The successful application of these methods has encouraged researchers to explore the use of NURBS in the analysis of discontinuous problems, to improve the representation of discontinuities and accurately describe their geometry and kinematics.
Discontinuities usually require an ad hoc treatment in numerical analysis; in the context of FEM, the conventional approach is to build a fitted or matching FE mesh, where the edges of FEs are aligned to discontinuities. This procedure is often costly and laborious, especially in three-dimensional (3D) problems and in cases where discontinuities evolve in time or have complex shapes. Mesh generators are often not robust in such cases and high refinements or special element arrangements-for instance, using quarter-point elements 8 -are needed to correctly capture field gradients. These drawbacks have pushed the development of methods such as the eXtended or Generalized Finite Element Method (X/GFEM), 9,10 where discontinuities are resolved elegantly by decoupling them from the mesh. Inspired by the possibility of including a priori knowledge of the solution into the FE formulation, 11 X/GFEM augments the standard FEM approximation by means of enriched terms that are able to capture jumps in the solution field (strong discontinuities) or in its gradient (weak discontinuities). This is achieved by means of suitable discontinuous enrichment functions, with corresponding enriched DOFs, which are associated with the original mesh nodes. However, X/GFEM is not without issues, as certain enrichment functions give rise to spurious terms in blending elements (uncut elements with enriched DOFs) that erode the solution's accuracy. 12 Moreover, nonzero Dirichlet BCs at nodes of elements containing discontinuities are enforced weakly (on average) by means of Lagrange multipliers or other techniques. 13 Finally, X/GFEM may result in ill-conditioned system matrices when discontinuities get arbitrarily close to the nodes of the mesh. As a result, recent work has focused in the formulation of stable generalized finite element methods, which aim at ensuring that the condition number of the resulting system matrices grow at the same rate as those of standard FEM. [14][15][16] The Discontinuity-Enriched Finite Element Method (DE-FEM) was recently introduced by Aragón and Simone 17 and offers a unified enriched approach for the analysis of both weak and strong discontinuities. It differs from X/GFEM in that enriched DOFs are assigned to newly generated nodes collocated at the intersection between discontinuities and elements' edges. Because of this, DE-FEM has some advantages over X/GFEM. Enrichment functions in DE-FEM are constructed simply from Lagrange shape functions in elements that are created for quadrature of local arrays (integration elements) and are therefore local to cut elements by construction. Enrichments thus vanish at mesh nodes and are exactly zero at all edges which are not crossed by discontinuities; as a result, standard DOFs maintain their physical meaning and essential BCs can be prescribed strongly as in standard FEM. In absence of strong discontinuities, DE-FEM simplifies to the Interface-enriched Generalized Finite Element Method (IGFEM) 18 or its successor the Hierarchical Interface-enriched Finite Element Method (HIFEM), 19,20 that were devised to analyze problems with weak discontinuities alone and that exhibit the same aforementioned advantages. Since its inception, DE-FEM has been extended to the analysis of 3D problems 21 and immerse boundary (fictitious domain) problems. 22 The latter work showed that the enriched formulation is capable of recovering smooth traction profiles at Dirichlet boundaries, something that is not possible with X/GFEM even with stabilization. 23 While the vast majority of literature on X/GFEM and interface-enriched FEM approximates geometries piecewise linearly, higher-order piecewise polynomial parametrizations have been proposed to approximate curved discontinuities. 12,[24][25][26] These methods are capable of achieving optimal or close-to-optimal convergence rates, provided that the geometries of discontinuities are described with sufficient accuracy by the chosen polynomial functions. F I G U R E 1 Mathematical representation of a cracked solid composed of material phases Ω m and Ω i , which represent matrix and inclusion, respectively. The parts of the domain boundary where we prescribe the solution field (Γ u ) and tractions (Γ t ) are also illustrated. The schematic contains both weak and strong discontinuities Γ i and Γ c , respectively Alternatively, recent works have incorporated NURBS to improve the geometric representation of discontinuities. Extended IGA was proposed as a combination of XFEM and IGA 27,28 to model nonmatching discontinuities in IGA. Legrain 29 introduced the NURBS-enhanced approach of NEFEM into XFEM by applying a NURBS-based mapping to describe discontinuities within standard FEA. In the context of IGFEM, the approach followed by Tan et al. 30 and Safdari et al., 31 called NURBS-enhanced IGFEM, utilizes NURBS both for geometric mapping and for constructing enrichment functions. Soghrati and Merel, 32 instead, restrict the use of NURBS to the geometric mapping alone, while the original IGFEM enrichment functions are used for the field approximation. Both approaches to combining IGFEM and NURBS were used to study weak discontinuities in composite materials 31,32 and microvascular systems, 30 and in shape optimization. 33 In these works, comparisons between classical FEM or IGFEM and NURBS-enhanced methods are mostly limited to verification examples, including the integration of areas, convergence studies, and simple benchmark problems. In the context of low-order FEM, where linear geometric approximations tend to be inaccurate in most cases, convergence studies on curved inclusions 31 and straight microchannels 30 based on bilinear elements already show improvements, with the greatest benefits on coarse meshes. However, there is a lack of works that offer an unbiased discussion on whether accuracy improvements can effectively help to reduce the number of DOFs and outweigh the complexity of the formulation and computer implementation. In addition, within the context of discontinuity-enriched methods, no research has been published on NURBS-enhanced problems containing strong discontinuities.
In this work, and following the work of Soghrati and Alvarez, 32 we investigate the use of NURBS to handle NURBS-based discontinuities in two dimensions for both IGFEM and DE-FEM. We provide insight into the potential and the limitations of combining NURBS with DE-FEM based on linear approximations of the field. The resulting methodology is then compared with standard DE-FEM, by considering factors such as accuracy and implementation efforts. We show that, while optimal convergence rates are still recovered for weakly discontinuous problems with no singularities, the recovery of stress intensity factors (SIFs) has no major improvement over standard DE-FEM, and the NURBS-based formulation has troubles reproducing rigid body modes and does not pass an immersed patch test.

FORMULATION AND IMPLEMENTATION
In this article, we are solely concerned with elastostatic boundary value problems in two dimensions. Consider a solid body Ω ∈ R 2 with closure Ω, defined on a Cartesian coordinate system spanned by {e 1 , e 2 }. Ω is decomposed into disjoint regions such that Ω = ∪ i Ω i and Ω i ∩ Ω j = ∅, i ≠ j (see Figure 1 for matrix and inclusion phases, Ω m and Ω i , respectively). Phase boundaries, defined as Ω i ≡ Γ i = Ω i ⧵ Ω i , are given an orientation defined by their corresponding outward unit normal n i . The domain boundary Γ ⊂ Γ m , with normal n, is decomposed into a region Γ u where the solution field u is prescribed, and a region Γ t where the traction field t is applied, such that Γ = Γ u ∪ Γ t , and Γ u ∩ Γ t = ∅. Internal to the domain, ∩ i Ω i ⧵ Γ identifies a series of weak or strong discontinuities. Both types of discontinuities, and even the boundary Γ, might be described using NURBS. In this work, we restrict ourselves to traction-free cracks, although DE-FEM can also treat cohesive cracks. 17 In addition, we restrict ourselves to non-intersecting discontinuities, although it is possible to handle intersecting discontinuities by means of the hierarchical approach used in HIFEM. 19 We aim to solve the linear elastostatic boundary value problem on Ω, with prescribed displacement u ∶ Γ u → R 2 and applied traction t ∶ Γ t → R 2 boundary conditions. The body force in the jth phase is denoted by b j ∶ Ω j → R 2 . Let (1) denote the vector-valued function space, where  2 is the space of square-integrable functions and the components of v| Ω i belong to the first-order Sobolev space  1 (Ω i ); note that v ∈  also satisfies homogeneous essential BCs on Γ u . Let u = v +ũ, withũ| Γ u = u and v ∈ , be the displacement field to account for non-homogeneous essential BCs. The weak or variational form can then be stated as: where the linear and bilinear forms are given, respectively, by and In Equations (3) and (4) it is implied that w j ≡ w| Ω j , and similarly for other subscripted quantities. We assume a linearly elastic material model following Hooke's law = C ∶ , where C is a fourth-order constitutive tensor and = 1 2 ( u + u ⊺ ) the linearized strain tensor. We solve the finite-dimensional version of Equation (2) by discretizing the domain Ω into finite elements so that Ω ≈ Ω h = ∪ i e i , e i ∩ e j = ∅ ∀i ≠ j. As shown in Figure 2, the FE mesh completely disregards the placement of discontinuities. In order to properly resolve such mismatch, the standard FEM approximation space is enriched or augmented by functions that capture the jump in the displacement field, and the gradient thereof. Following a Bubnov-Galerkin procedure, we choose both trial and weight functions from the DE-FEM space  h e : 17 where the first term corresponds to the standard FEM component and the other terms constitute the enrichment. In Equation (5), h is the index set referring to all mesh nodes with associated Lagrange shape function N i and standard DOFs U i . Similarly, w ( s ) is the index set that refers to weak (strong) enriched nodes, and associated to enriched DOFs i ( i ). Their corresponding enrichment functions i ( i ) are constructed by using Lagrange shape functions that are mapped to integration elements with curved edges through a NURBS mapping, as discussed in Section 2.1. Enrichment functions, which are nonzero only in cut elements, have their maximum absolute value at the location of their corresponding enriched node and ramp linearly to zero at mesh nodes. A schematic representation of a strongly discontinuous solution field is shown in Figure 3 as the sum of standard and enriched contributions.
In elements containing crack tips, discontinuities cut element boundaries only at one point where an enriched node is created and associated with both weak and strong enrichment functions. The enriched node created at the crack tip, on the other hand, is only associated with a weak enrichment in order to guarantee continuity in the displacement. 17 Equation (2) may be augmented to include crack tip enrichment functions that model the singular stress field, which can potentially improve the accuracy of the computed solution. However, these enrichment functions are prone to ill-conditioning and they may lead to loss of accuracy in blending elements if not treated properly. 34 In addition, these singular enrichments increase the method's complexity because special techniques are required for the numerical quadrature of local arrays. 35 Therefore, singular crack tip enrichments were not considered in the original DE-FEM formulation and are also not considered in this work. The quadrature of local arrays in (integration) elements around a crack tip is therefore no different F I G U R E 2 Schematic representation of the FE discretization and its interaction with the discontinuities. The mesh is intersected with the NURBS discontinuities and integration elements (some with NURBS edges, shaded in a darker color) are created. Integration elements are then stored in an ordered tree data structure F I G U R E 3 Illustration of a discontinuous solution field and its composition from a continuous field, computed from standard shape functions, and a discontinuous field, obtained from the contribution of the enrichment. On the element, the curved discontinuity and the enriched nodes x 4 and x 5 are highlighted in red F I G U R E 4 Element containing a crack tip and the creation of integration elements. NURBS-based integration elements are shaded in a darker color from the one required in elements away from it. An example of a mesh-discontinuity interaction at a crack tip is shown in Figure 4, showing the creation of integration elements.
The creation and bookkeeping of integration elements are straightforward by using an ordered tree data structure, as shown in Figure 2. Numerical integration is conducted on the leaf elements of the tree to obtain k i and f i , the local stiffness matrix and local force vector of the ith element, respectively. We arrive to the discretized system of linear equations KU = F, where the stiffness matrix K and the force vector F are obtained from the assembly of element local arrays: with A denoting the FE assembly operator. For more details on the formulation see Aragón and Simone. 17

NURBS description of discontinuities
A DE-FEM analysis entails creating a nonmatching (usually structured) mesh and the definition of discontinuities. In this work, discontinuities are represented by NURBS curves in the parametric coordinate s: 36 where c(s) is the NURBS curve describing the discontinuity, R p i (s) are NURBS basis functions, defined as rational functions of B-spline basis functions of order p, and b i are the associated control points; there are an equal number n of NURBS basis functions and control points. A NURBS curve is described with the aid of a vector of coordinates s i , called knots; each NURBS basis function will have its support on one or more knot intervals, depending on its order p.
Intersection points between the background mesh already defined and the NURBS are found by a Newton-Raphson solver. Pseudocode for this mesh interaction is found in Appendix A. In correspondence with these intersection points, new enriched nodes are created: one single node for weak discontinuities and double nodes for strong discontinuities. Crack tips are an exception: the location for the tip node is found by checking whether the end point of the discontinuity is contained within an element where only one edge is intersected. In addition to the creation of enriched nodes, the NURBS segment contained within each element is extracted. This is accomplished by reducing the continuity of the parametrization at enriched nodes via repeated knot insertion, according to the procedure described in more detail in Appendix B, until the segment can be described with only self-contained information. 30 The extraction of the NURBS segment contained in an element involves then selecting all parameters describing that segment, notably the knots s i within the span of the enriched nodes and the control points that affect the NURBS basis functions in that range. An example to demonstrate this procedure is also included in Appendix B. This step facilitates the use of data associated with each spline segment, which is stored element-wise in a suitable data structure and later used in successive stages of analysis and postprocessing. Subsequently, every element that is crossed by the discontinuity is split into integration elements by means of Delaunay triangulation 36 and stored in an ordered tree as explained earlier.
At this point, both standard shape and enrichment functions can be constructed. For integration elements that have a NURBS discontinuity edge, such as the ones represented in a darker color in Figure 4, a NURBS-based mapping is required. In the proposed implementation, the enrichment functions are first constructed as Lagrange shape functions N □ i on a square reference domain, which requires fewer integration points than the choice of a triangular reference domain. 37 The reference domain is chosen to be the square [− 1, 1] × [− 1, 1] in reference coordinates and . The enrichment functions thus constructed can be mapped from the square reference element to a triangular element by collapsing the top edge of the square into a single vertex x 3 . Hence, N △ 3 , the shape function of the integration element associated with that node, is obtained by adding N □ 3 and N □ 4 , the shape functions associated with the top nodes of the square reference domain. 37 The triangle shape functions are therefore given by where the NURBS-based mapping between the square and the triangular domain, which is shown in Figure 5, is given by 32,37 In mapping (9), s is the parametric coordinate of the NURBS curve, while t is an auxiliary coordinate, here assuming values in [0, 1]. NURBS-based enrichment functions i and i from Equation (5) are thus created by applying mapping (9) to the Lagrangian shape functions (8). In (9), it has to be taken into account that splines are oriented curves; therefore, consistency is maintained with the orientation of the element boundary when performing numerical integration. It should also be observed that the mapping (9) only applies to elements with a single NURBS edge. For integration elements with multiple NURBS edges, for example, created by intersecting discontinuities, a dedicated mapping needs to be introduced that can account for the parametrization of multiple splines. On top of this, additional processing tasks are required to handle geometric operations and store data for different NURBS discontinuities, resulting in a significant increase in complexity compared with DE-FEM with linearly approximated discontinuities. The integration of shape and enrichment functions, needed to compute the respective contributions to the system arrays, is performed by numerical quadrature. The numerical integration of splines is an open field of research, as in general it is not possible to integrate NURBS exactly by Gauss quadrature. 31 For consistency with standard FEM elements, here we adopt high-order Gauss-Legendre quadrature rules for integrating elements containing a NURBS edge. These are widely used for integrating NURBS functions 2,30,32,37 and have been found to be appropriate in terms of both accuracy and efficiency. 6 Gauss points ( i , i ), defined over the square reference domain, are mapped onto the integration elements via the mapping given by Equation (9). The integration order in the two directions is decoupled to reduce the total number of Gauss points, based on considerations on the linearity of the mapping along . 38 We use six integration points in the direction along the NURBS edge and two points in the linear direction.
An important remark at this point involves self-intersecting integration elements 29,32 due to the subtriangulation of parent elements, as depicted in Figure 6(A). Generally, a discontinuity splits an intersected element into a triangular and a quadrilateral part, unless it passes through one of the vertices. Based on the choice to handle only triangular children elements, the quadrilateral part needs to be subdivided further, for example, by taking one of its diagonals as a splitting edge. This procedure may result in the construction of splitting edges that are secant to the curved discontinuity, which in turn could lead to self-intersecting (or tangled) integration elements with negative Jacobians. This problem can appear particularly in the presence of discontinuities with relatively high curvature and can compromise the overall accuracy of the method. Although several strategies can be adopted to prevent the appearance of tangled elements, they all require additional implementation effort, limiting the attractiveness of the method.
A first approach consists in finding an alternative triangulation that does not construct degenerate elements ( Figure 6(B)); this straightforward strategy is frequently sufficient to address the issue. However, as also highlighted by Legrain,29 such an intersection-free triangulation does not always exist. In such case, one solution is to refine the triangulation on the parent element by placing additional nodes on the discontinuity, until a triangulation is found that is free from tangled elements. This form of adaptive refinement, shown in Figure 6(C), comes at the price of additional checks for intersections and for finding the best location for new enriched nodes, eventually resulting in increased computational cost. An alternative form of adaptive refinement is proposed by Legrain,29 who refines the parent element by introducing new nodes on its edges, as depicted in Figure 6(D); the creation of more enriched nodes follows from the higher number of intersections with the discontinuity. The drawbacks of this method include an increase in the processing tasks, the introduction of hanging nodes, and the necessity for local mesh refinements, 32 which is what methods that decouple mesh from discontinuities aim to avoid.
Alternatively, Soghrati and Merel 32 propose the use of NURBS as subdivision edges: placing control points on both the discontinuity and the element edges ensures that the spline is strictly contained in the hull defined by the control points, thus preventing intersections with the discontinuity itself. The method, as shown in Figure 6(E), is robust by construction, although it increases the complexity of the geometric operations, the number of NURBS edges to be processed, and the amount of spline data that needs to be stored. In fact, integration elements with two NURBS edges require a dedicated mapping, capable of handling the parametrization of both splines.
A final approach to the issue of tangled elements consists in admitting nontriangular integration elements, implying that parent elements may be split into children elements by means of discontinuities only, as in Figure 6(F). In this way, the construction of additional subdivision edges becomes unnecessary and the possibility of intersections with discontinuities is eliminated. The downside of this strategy would be the need to account for both triangular and nontriangular integration elements in the mapping; a formulation for quadrilateral integration elements can be derived from the blending mapping proposed by Szabó et al. 7

Note on rigid body displacements
The formulation described above is superparametric, that is, the description for the internal NURBS edges is of higher order than that of the field interpolation. It is well known that superparametric elements suffer from artificial stiffness, and therefore, they cannot properly describe rigid body displacements. 3 In order to investigate the accuracy of DE-FEM with NURBS discontinuities for rigid body displacements, we analyze the lowest energy modes of a single background element with a NURBS interface, by performing an eigenvalue analysis on the stiffness matrix.  Figure 7(A)-(C) illustrate the eigenmodes for the unconstrained problem with different material properties on either side of the interface: E 1 = 10 and 1 = 0.3 on the left of the interface, and E 2 = 1 and 2 = 0.3 on the right. The three modes are linear combinations of the two rigid body translation and a rigid body rotation. It is found that the corresponding eigenvalues are indeed zero to numerical precision. However, in this situation, the enriched DOFs are zero, as there is no jump in the gradient of the field in a rigid body displacement of the full element. For more rigorous testing, we place the NURBS interface in an immersed setting, 22 where a material with properties E 1 = 10 and 1 = 0.3 is placed on one side of the interface, and the other side is a void, and the original nodes of the void domain are fixed. This is shown in Figure 7(D)-(I), where we expect, again, to recover three zero-energy modes. However, the three lowest eigenvalues are no longer exactly zero, indicating that the NURBS enrichment introduces artificial stiffness to the formulation ( Table 1).
The test shown in Figure 8 confirms once more that DE-FEM with NURBS-based enrichments indeed suffers from this numerical defect. In this test, a square mesh consisting of 3 × 3 × 2 triangular elements is defined between coordinates The superparametric "variational crime" severely limits the attractiveness of NURBS-enrichments in DE-FEM with low-order (linear) background elements. However, hp-refinement could be employed to alleviate this problem, 7 by increasing the order of the field interpolation, while simultaneously reducing the ratio between the curvature of the F I G U R E 7 Lowest-energy modes of a single background element with a NURBS interface. (A)-(C) show the rigid body modes of an element with a NURBS interface between two different materials, where the undeformed configuration is shown with dashed lines. In (D)-(F), the domain on the left of the NURBS is void, and therefore the nodes on the left are fixed. The lowest-energy modes correspond to a nonzero eigenvalue, that is, the displacements are not rigid. In (G)-(I), the region on the right is void, so a single node is fixed. The mode in (G) is a rigid rotation of the full background element, and therefore the energy associated with this mode is zero. However, the modes in (H) and (I) again correspond to nonzero eigenvalues NURBS and the element size. Alternatively, the adoption of an isoparametric formulation could be investigated, to ensure that the parametrization of both solution field and geometric representation have the same order.

Convergence study
The convergence of the method is first studied with the Eshelby inclusion problem, as shown in Figure 9. The problem consists of a circular plate (radius r 2 = 2) containing a circular inclusion (radius r 1 = 0.4). The material elastic constants are chosen as E 1 = 1 and 1 = 0.25 for the inclusion, and E 2 = 10 and 2 = 0.3 for the matrix. The plate is subjected to a prescribed displacement The exact displacement field for this problem is given as 39 : [ where the parameter , dependent on the Lamé parameters and of the two materials, is defined as: We use the standard  2 and energy norms of the relative error to study convergence, given by and respectively, where u is the analytical solution given by Equation (10) and u h is the numerical solution.

F I G U R E 10 Discretization (A) and
deformation (B) of the Eshelby inclusion problem at the coarsest level. The discretization shows that the inclusion is modeled exactly as a circle even with a very coarse mesh; however, the deformed inclusion loses its perfectly circular shape, due to the linear interpolation of the solution field  Figure 10(A) shows the coarsest discretization used to study convergence and the deformation at this mesh size. Noteworthy, although the geometry is represented exactly, the field approximation is linear and therefore the deformed inclusion does not remain perfectly circular. Considering a higher polynomial interpolation for the field would yield considerable improvements. Figure 11 summarizes the results of the convergence study. The relative errors in energy and  2 norms are plotted as a function of mesh size h and the total number of DOFs. We compare the standard FEM solutions with matching meshes, DE-FEM with linear approximation of the interfaces, and DE-FEM with NURBS discontinuities. Only the inner boundary is modeled with discontinuity enrichments; the outer boundary, instead, simply consists in a piecewise linear approximation of the boundary (as done with matching meshes). The curves show that the NURBS-based DE-FEM reaches an optimal convergence rate. The improvement in accuracy compared with the standard FEM and the linearly approximated DE-FEM is only marginal and decreases with mesh refinement. This result is not surprising: the improved accuracy in elements with NURBS edges becomes less and less dominant for increasingly refined meshes, as the linearly approximated boundary tends towards the shape of a circle. Figures 12 and 13 investigate in more detail the accuracy for the different methods. Figure 12 shows, at the second coarsest mesh size, the stress field computed with standard FEM with a uniform mesh, piecewise-linear DE-FEM, and DE-FEM with NURBS discontinuities. The results show that both standard DE-FEM and NURBS-based DE-FEM are able to capture the high stresses near the interface. Although similar, there are some differences in the values of the computed stresses. These differences have direct consequences on the errors computed for the different solutions. Figure 13 shows the field for both error norms at the second coarsest mesh size for the three approaches. In all methods the error is localized near the material interface due to the circle's discretization error. However, as the circle is better represented in DE-FEM, and to an ever higher extent in NURBS-based DE-FEM, the error is considerably lower in these enriched methods.

Stress intensity factors
With this example we aim at studying the capabilities of DE-FEM with NURBS discontinuities in the treatment of curved crack problems. The problem, represented schematically in Figure 14, consists of an infinite plate under biaxial tension containing a circular arc crack, for which an analytical expression for the SIFs is known. [40][41][42][43] Due to symmetry, we analyze half of the domain of interest with L = 5. Symmetry BCs are prescribed on the bottom and the analytical displacement is prescribed on the remaining sides. 44 A Young's modulus E = 1000 and Poisson's ratio = 0 are considered.

F I G U R E 14
Problem definition for the study of stress intensity factors. An infinite plate contains a circular arc crack with radius r that spans an angle 2 . The plate is subjected to biaxial tension 0 . The domain of interest Ω for the numerical study is denoted by a darker color Under the assumption of an infinite plate subjected to uniform tension at infinity, the SIFs are given by: 1 + sin 2 ( ∕2) ) .
In our case, uniform biaxial traction is applied, therefore 11 = 22 = and 12 = 0. Our domain of interest Ω is discretized using a structured mesh of triangular elements defined on a 199 × 99 grid. SIFs are then computed by using the domain integral method, 17,45 which takes the contribution from elements inside an integration domain in the neighborhood of the crack tip-we take a circular area with radius 3h. SIFs are reported for nine different choices of , ranging from 10 • to 90 • . The radius r is kept at a constant value of 1, with the arc cord changing to satisfy the assigned angle. SIF values are also normalized as K i = K i ∕ √ a, i = I, II, to make them independent of the specific choices of dimensions and loading.
The results are reported in Table 2 and visualized in Figure 15, where computed SIFs for both standard DE-FEM and NURBS-based DE-FEM are compared with the analytical values obtained from Equations (14) and (15). For illustration purposes, the deformed configuration and the von Mises stress distribution are shown in Figure 16 on a coarser discretization than that used for extracting SIFs.
It can be seen that numerical results match the analytical values closely, particularly for mode I, which presents a relative error in the order of 1%. SIFs for mode II, however, are less accurate, with an error that is as high as 8.79% for an angle of 10 • . Not including singular enrichments into the formulation can be identified as a limiting factor to the solution's accuracy and to the possibility to achieve optimal convergence for singular problems in fracture mechanics. By comparing the results between NURBS-based DE-FEM and standard DE-FEM, we observe that the former does not introduce significant improvements: this suggests that, in this case, the geometric discretization error does not affect significantly the stress field near crack tips and therefore their stress intensity.
It has been observed that one factor potentially affecting the accuracy of the computed field gradient is the aspect ratio of the integration elements; 21,22 the creation of enriched nodes too close to existing mesh nodes may result, in some cases, in nonphysical stress concentrations; these are more pronounced near material interfaces 46,47 rather than along Dirichlet boundaries 22 or traction-free cracks. 21 The occurrence of tangent discontinuities or close-to-coincident enriched TA B L E 2 Normalized K I and K II , obtained analytically, with NURBS-based DE-FEM, and with standard DE-FEM

F I G U R E 16 Deformation (A) and von
Mises stress (B) of the curved crack spanning a 50 • angle. The figures, obtained with a relatively coarse discretization, illustrate also the interaction between the mesh and the NURBS discontinuity and mesh nodes is more frequent for very refined meshes, which was observed also by Safdari et al., 31 and for shapes that resemble the pattern of the mesh. In these cases, a slight increase in the error may occur.
From this example, it can be concluded that DE-FEM with NURBS discontinuities can lead to accurate results when modeling curved cracks and computing the associated SIFs, although the method can be sensitive to poor aspect ratios in integration elements and to the low order of the field approximation. With the flexibility of NURBS-based shapes, not only simple curves like circumference arcs can be treated: applications can be extended to arbitrary more complex shapes.

F I G U R E 18
Patch test on a complex-shaped immersed domain: a duck-shaped domain is immersed in a rectangular mesh, where the domain outside the duck is void. A uniform expansion is applied to the duck on its boundary Γ i , which is expected to result in a uniform state of stress. In (A), the geometry of the domain is piecewise-linearly approximated, and a constant state of stress is indeed obtained. In (B), a NURBS representation of the domain is employed. It is found that this patch test fails due to the superparametric formulation

Immersed patch test
Finally, to illustrate the shortcomings of the formulation for linear background elements, an immersed patch test on a complex shape is performed. Taking advantage of the flexibility of NURBS in treating arbitrary geometries, a duck-shaped boundary Γ i has been designed, based on a second-order NURBS curve. The boundary has been superposed on a rectangular mesh and two distinct domains were defined: a homogeneous solid domain (E 1 = 2, 1 = 0 ) inside the boundary and a void domain (E 2 = 0, 2 = 0 ) outside. Figure 17 illustrates this immersed problem schematically. On the boundary Γ i , a radial displacement u(r, ) = [−0.5r, 0] ⊺ is prescribed; following the immersed Dirichlet BCs described in van den Boom et al. 22 that mimics a uniform compression of the duck that leads to a uniform state of stress inside it. All original mesh nodes in the void region are fixed. Figure 18 illustrates the result of this problem solved with the piecewise enriched formulation in Figure 18(A), and with the NURBS-based formulation in Figure 18(B). Noteworthy, this example is devoid of strong discontinuities, so the NURBS-based formulation is compared with that of the linear IGFEM. It is clear that the NURBS-based method cannot reproduce the constant state of stress, as artificial stiffness is introduced in the NURBS integration elements, as was also outlined in Section 2.2. The NURBS-based formulation does not pass the immersed patch test.

SUMMARY AND CONCLUSIONS
In this article, we studied the use of NURBS-based discontinuities in combination with DE-FEM to represent exactly the geometry of both weak and strong discontinuities. Building upon standard DE-FEM, where the geometry is approximated piecewise linearly, the adoption of NURBS requires only slight modifications. The differences lie mostly in the nature of the geometric operations needed to find intersections and in the integration of elements intersected by discontinuities, where the NURBS-based mapping is used. Furthermore, some specialized bookkeeping is required.
NURBS-based DE-FEM was tested for a variety of problems in linear elasticity. For weak discontinuities without singularities-where in the absence of strong discontinuities the formulation falls back to NURBS-based IGFEM-we showed that the method successfully recovers optimal convergence rates. The numerical examples show the versatility of the method for handling discontinuities of varying complexity. The simulations performed on fracture problems with arc cracks showed that DE-FEM recovers accurate results in the computation of SIFs, with an error in the order of 1% for K I . Yet, these results were not visibly different from those obtained for standard DE-FEM, confirming that errors in the geometric representation of the cracks had little influence on the singularity at crack tips. Furthermore, we have shown that NURBS-based enriched formulation does not recover rigid body modes in an immersed setting and does not pass an immersed patch test. In fact, adopting a NURBS-based mapping results in a superparametric formulation, which cannot properly describe rigid body displacements and constant states of stress, 3

as demonstrated in the proposed examples.
While this work focused on the analysis of problems with stationary NURBS, more interesting applications may be found in boundary evolution problems. The proposed method could be extended to treat such problems by performing additional operations to modify NURBS discontinuities between iterations. In problems involving moving boundaries, control points can be used to manipulate a NURBS discontinuity and impose the desired displacement on the boundary. 2 The coordinates of control points can be also conveniently chosen as design variables in shape optimization. 33,48 On the other hand, crack growth problems require extending the discontinuity while preserving the shape of the preexisting crack. This can be achieved, for example, by stitching a NURBS crack increment to the original fracture or by knot insertion near the tip, and subsequent displacement of control points to propagate the crack. 49,50 In boundary evolution problems in general, the local level of detail and freedom may be increased in subsequent steps by refining the knot vector by knot insertion and, concurrently, increasing the number of control points.
On first sight, NURBS-based DE-FEM preserves all the key advantages of the standard DE-FEM formulation: the enrichments are local by construction, with enriched DOFs assigned to nodes collocated along discontinuities. In particular for strong discontinuities, the interpretation of the enriched DOFs as the extent of the crack opening is preserved. The original mesh nodes also retain their physical interpretation, facilitating the strong enforcement of Dirichlet BCs. Finally, we showed that extending the method to treat immersed problems was straightforward.
In standard DE-FEM, mesh refinements are required to improve the quality of both the solution field and the geometry. In NURBS-based DE-FEM, geometries are accurate even for the coarsest meshes and thus no refinements are needed for the sole purpose of improving shape representation. Nonetheless, refinements are generally needed to achieve an accurate representation of stresses. Often, this requires meshes that are fine enough to attenuate the advantages of employing an exact geometry compared with a very refined linear approximation. This holds in particular when using linear shape functions: as verified in our implementation, this imposes a considerable limitation to the maximum accuracy that can be achieved. The full potential of DE-FEM with NURBS discontinuities will, therefore, only be exploited when using the technique on higher-order interpolations.
The study of NURBS-based DE-FEM on low-order elements brought the attention also to several disadvantages which are not expected to be eliminated by higher-order field approximations. The NURBS-based mapping of enrichments comes at the price of increased implementation complexity and computational cost for handling geometric operations on splines. For instance, suitable data structures are needed to store element-wise sections of NURBS discontinuities, which are then used for integration purposes and in the postprocessing stage. The complexity of these tasks would increase even further on elements with multiple NURBS edges, created, for example, by intersecting NURBS discontinuities, which would not only require additional processing and storage, but also a dedicated mapping of enrichments.
The introduction of splines in the formulation raises a number of issues also when it comes to geometric operations and numerical quadrature. The presence of curved boundaries opens the possibility to create undesired self-intersecting integration elements that may in turn result in negative Jacobians if not processed correctly. Another factor affecting the sign of the Jacobians lies in that splines are oriented curves by construction: this needs to be taken into account for the purposes of numerical integration, to ensure that the element boundary is oriented coherently. Moreover, the numerical integration of NURBS-based domains cannot be performed exactly by Gauss quadrature: this means that a higher number of points is needed to achieve the desired precision, even in case of low-order splines. Finally, the formulation for a NURBS-based mapping of low-order elements is superparametric, which introduces the aforementioned issues when computing rigid body displacements and constant states of stress.
All in all, the convenience of introducing spline-based modeling will depend on the governing physics, the geometry of the discontinuities, and the order of the approximation used. However, this study brought the attention to several drawbacks associated with this more complex and costly formulation, which do not appear to be justified by the potential minimization of geometric error and mesh refinement. This imposes significant limitations to the convenience and attractiveness of using NURBS geometries for low-order interpolations.

APPENDIX A. COMPUTER IMPLEMENTATION
A comprehensive discussion on the computer implementation aspects of DE-FEM was presented in Aragón and Simone. 17 Here, we outline only the changes that needed to be made to include NURBS in the method; these curves intersect FEs of the background mesh, and as a result, the latter are split in subdomains-triangles in fact-that are used only for integration of local arrays.
The operations performed on NURBS geometries when intersecting discontinuities with the underlying mesh and, then, constructing children elements are summarized in Algorithm 1. A first check verifies whether a mesh element overlaps with the discontinuity bounding box. If this is the case, the edges of the element are considered one by one, to determine whether they intersect the discontinuity. For this purpose, a sign function associated with the discontinuity is defined, such that its sign changes across the discontinuity. Finding opposite signs at the nodes that delimit an edge indicates that the edge has an intersection point; a Newton-Raphson solver finds the coordinates of the intersection, where an enriched node is then created. If a single intersection point is found on the element boundary, the algorithm looks for a discontinuity tip (e.g., a crack tip) inside the element, to create a second enriched node. To conclude, the algorithm extracts the parametrization of the NURBS segment contained in the element and updates the data for the integration elements accordingly. The construction of local arrays in integration elements follows a similar approach to that used in standard elements. As shown in Algorithm 2, the local arrays are first initialized and then updated with the contribution of the integration points, which is obtained from shape and enrichment functions, their derivatives, and the Jacobian. An ordered tree data structure  is passed to the function that contains information on the splitting of mesh elements into integration elements by discontinuities. The differences with respect to integration elements in the originally proposed DE-FEM are mainly in the Jacobian of the transformation, that uses NURBS functions, and in the way enrichment functions are calculated. The other parameters to the function follow the standard FEM approach.

Algorithm 1. Interaction between mesh and NURBS discontinuity
Pseudocode for the call that obtains shape and enrichment functions is given in Algorithm 3. The function summarized iterates over enriched nodes to obtain the contribution from the corresponding enrichment functions. Here, a NURBS-based mapping is used to map the integration point from reference to global coordinates and to compute the enrichment functions and their derivatives. The contribution of enrichments and shape functions, the latter obtained from the parent element, are stacked into matrices and returned along with the Jacobian determinant.

APPENDIX B. KNOT INSERTION
Knot insertion is a form of refinement where one or more knots are added to a spline without changing its shape. 36 The insertion of coincident knots is of particular interest for NURBS-based enriched FEM because it allows decreasing the continuity of the parametrization across the corresponding point: null intervals between knots reduce the blending range of basic functions across their coordinates. In the general case, a NURBS curve with order p is C p − k -continuous at a point mapped from a knot with multiplicity k. Reducing the continuity of the spline parametrization across the enriched nodes, by inserting knots with sufficient multiplicity to achieve C −1 -continuity, allows representing the NURBS segment within an element without retrieving external data.  In a NURBS curve, each point results from the contribution of several basis functions, blended within the range of a few knot intervals by means of control points, according to Equation (7). By construction, a NURBS basis function R p (s) of where w i is the weight factor associated with the ith B-spline basis function. A NURBS curve and the associated B-spline share the same knot vector, while to each control point b i = (x i , y i ) of the NURBS will correspond a higher-dimensional control pointb i = (x i , y i , w i ) of the B-spline. Starting from a B-spline curve with order p, a knot vector = [s 0 , … , s m−1 ] of length m and a set of n control points B = {b 0 , … ,b n−1 }, a new knot s ′ will be inserted, such that s ′ ∈ [s k , s k + 1 ] for an arbitrary integer k ∈ [0, m − 2]. Since each basis function spans a limited number of knot intervals, this insertion affects only a portion of the B-spline: the one that lies within the convex hull defined by {b k−p , … ,b k }. Outside this range, the curve and the control points will remain untouched. The knot insertion algorithm to determine the modified control points reads 36 with coefficients a i given by The newly computed control points will replace the previous points {b k−p+1 , … ,b k−1 } in the set of control points of the B-spline, which then becomes {b 0 , … ,b k−p ,b ′ k−p+1 , … ,b ′ k ,b k , … ,b n−1 }. The new control points and the associated weights for the corresponding NURBS are straightforwardly obtained from the control points of the B-spline.
The example reported in Figure B1 illustrates the transformation of the control points and parametrization of a NURBS arc after multiple knot insertion. New knots are inserted at s = 0.7, with multiplicity 3, to achieve C −1 -continuity at that point. In this way, the NURBS segment between knots s = 0 and s = 0.7 can be fully described by the knot vector