Discretization of Non-uniform Rational B-Spline (NURBS) Models for Meshless Isogeometric Analysis

We present an algorithm for fast generation of quasi-uniform and variable-spacing nodes on domains whose boundaries are represented as computer-aided design (CAD) models, more specifically non-uniform rational B-splines (NURBS). This new algorithm enables the solution of partial differential equations within the volumes enclosed by these CAD models using (collocation-based) meshless numerical discretizations. Our hierarchical algorithm first generates quasi-uniform node sets directly on the NURBS surfaces representing the domain boundary, then uses the NURBS representation in conjunction with the surface nodes to generate nodes within the volume enclosed by the NURBS surface. We provide evidence for the quality of these node sets by analyzing them in terms of local regularity and separation distances. Finally, we demonstrate that these node sets are well-suited (both in terms of accuracy and numerical stability) for meshless radial basis function generated finite differences discretizations of the Poisson, Navier-Cauchy, and heat equations. Our algorithm constitutes an important step in bridging the field of node generation for meshless discretizations with isogeometric analysis.


Introduction
The key element of any numerical method for solving partial differential equations (PDE) is discretization of the domain.In traditional numerical methods such as the finite element method (FEM), this discretization is typically performed by partitioning the domain into a mesh, i.e., a finite number of elements that entirely cover it.Despite substantial developments in the field of mesh generation, the process of meshing often remains the most time consuming part of the whole solution procedure while the mesh quality limits the accuracy and stability of the numerical solution [18].In contrast, meshless methods for PDEs work directly on point clouds; in this context, point clouds typically referred to as "nodes".In particular, meshless methods based on radial basis function generated finite difference (RBF-FD) formulas allow for high-order accurate numerical solutions of PDEs on complicated time-varying domains [29,14] and even manifolds [26].The generation of suitable nodes is an area of ongoing research, with much work in recent years [32,7,9,35,28].In this work, we focus primarily on the generation of nodes suitable for RBF-FD discretizations, although our node generation approach is fully independent of the numerical method used.
Node sets may be generated in several different ways.For instance, one could simply generate a mesh using an existing tool and discard the connectivity information [17].However, such an approach is obviously computationally expensive, not easily generalized to higher dimensions, and in some scenarios even fails to generate node distributions of sufficient quality [28].Another possible approach is to use randomlygenerated nodes [17,22]; this approach has been used (with some modifications) in areas such as compressive sensing [1] and function approximation in high dimensions [24].Other approaches include iterative optimization [12,19,15], sphere-packing [16], QR factorization [34], and repulsion [9,35].It is generally accepted that quasi-uniformly-spaced node sets improve the stability of meshless methods [36].In this context, methods based on Poisson disk sampling are particularly appealing as they produce quasi-uniformly spaced nodes, scale to arbitrary dimension, are computationally efficient, and can be fully automated [28,32,7].
A related consideration is the quality of the domain discretization.In the context of meshes, for instance, it is common to characterize mesh quality using element aspect ratios or determinants of Jacobians [11,38].Analogously, the node generation literature commonly characterizes node quality in terms of two measures: the minimal spacing between any pair of nodes (the separation distance), and the maximal empty space without nodes (the fill distance).Once again, in this context, Poisson disk sampling via advancing front methods constitutes the state of the art [32,7].More specifically, the DIVG algorithm [32] allows for variable-density Poisson disk sampling on complicated domains in arbitrary dimension, while its generalization (sDIVG) allows for sampling of arbitrary-dimensional parametric surfaces.DIVG has since been parallelized [5], distributed as a standalone node generator [30], and is also an important component of the open-source meshless project Medusa [33].
Despite these rapid advances in node generation for meshless methods (and in meshless methods themselves), the generation of node sets on domains whose boundaries are specified by computer aided design (CAD) models is still in its infancy.Consequently, the application of meshless methods in CAD supplied geometries is rare and limited either to smooth geometries [23] or to the use of surface meshes [6,10,13].In contrast, mesh generation and the use of FEM in CAD geometries is a mature and well-understood field [11,3].In our experience, current node generation approaches on CAD geometries violate quasi-uniformity near the boundaries and are insufficiently robust or automated for practical use in engineering applications.In this work, we extend the sDIVG method to the generation of variable-density node sets on parametric CAD surfaces specified by non-uniform rational B-splines (NURBS).We then utilize the variable-density node sets generated by sDIVG in conjunction with the DIVG method to generate node sets in the volume enclosed by the NURBS surface.Our new framework is automated, computationally efficient, scalable to higher dimensions, and generates node sets that retain quasi-uniformity all the way up to the boundary.This framework also inherits the quality guarantees of DIVG and sDIVG, and is consequently well-suited for stable RBF-FD discretizations of PDEs on complicated domain geometries.The remainder of the paper is organized as follows.The NURBS-DIVG algorithm is presented in Section 2 along with analysis of specific components of the algorithm.The quality of generated nodes is discussed in Section 3. Its application to the RBF-FD solution of PDEs is shown in Section 4. The paper concludes in Section 5.

The NURBS-DIVG Algorithm
CAD surfaces are typically described as a union of multiple, non-overlapping, parametric patches (curves in 2D, surfaces in 3D), positioned so that the transitions between them are either smooth or satisfying geometric conditions.A popular choice for representing each patch is a NURBS [27], which is the focus of our work.Here, we present a NURBS-DIVG algorithm that has three primary components: 1. First, we extend the sDIVG algorithm [7] (Section 2.2) for sampling parametric surfaces to sampling individual NURBS patches and also the union of multiple NURBS patches (Section 2.3.2).
2. Next, we deploy the DIVG algorithm [32,5] in the interior of the domain using the sDIVG generated samples as seed nodes (Section 2.1).

3.
To ensure that DIVG generates the correct node sets in the interior of the domain whose boundary consists of multiple parametric NURBS patches, we augment sDIVG with a supersampling parameter (Section 2.3.3).
In following subsections, we first briefly present the DIVG algorithm.We then describe the sDIVG algorithm, which generalizes DIVG to parametric surfaces, focusing on sampling a single NURBS surface.We then describe how the sDIVG algorithm is generalized to surfaces consisting of multiple NURBS patches, each of which have their own boundary curves.Finally, we describe the inside check utilized by our algorithm needed to generate nodes within the NURBS patches, and the complications therein.

The DIVG algorithm
We now describe the DIVG algorithm for generating node sets within an arbitrary domain.As mentioned previously, this algorithm forms the foundation of the sDIVG and NURBS-DIVG algorithms.
DIVG is an iterative algorithm that begins with a given set of nodes called "seed nodes"; in our case, these will later be provided by NURBS-DIVG.The seed nodes are placed in an expansion queue.In each iteration i of the DIVG algorithm, a single node p i is dequeued and "expanded".Here, "expansion" means that a set C i of n candidates for new nodes is uniformly generated on the sphere centered at the node p i , with some radius r i and a random rotation.Here, r i stands for target nodal spacing and can be thought of as derived from a density function h, so that r i = h(p i ) [32,31].Of course, the set C i may contain candidates that lie outside the domain boundary or are too close to an existing node.Such candidates are rejected.However, the candidates that are not rejected are simply added to the domain and to the expansion queue; this is illustrated in Figure 1.The iteration continues until the queue is empty.A full description of DIVG can be found in [32].Its parallel variant is described in [5].

The sDIVG algorithm
The sDIVG algorithm is a generalization of DIVG to parametric surfaces.Unlike DIVG (which fills volumes with node sets), sDIVG instead places nodes on a target parametric surface in such a way that the spacing between nodes on the surface follows a supplied density function.While other algorithms typically achieve this through direct Cartesian sampling and elimination [37,28], the sDIVG algorithm samples the parametric domain corresponding to the surface with an appropriately-transformed version of the supplied density function.More concretely, given a domain Ω ⊂ R d , sDIVG iteratively samples its boundary ∂Ω ⊂ R d by sampling a parametrization Λ of its boundary instead.The advantage of this approach over direct Cartesian sampling is obtained from the fact that Λ ⊂ R d−1 (or S d−1 ) is a lower-dimensional representation of ∂Ω, leading to an increase in efficiency.
We now briefly describe the density function transformation utilized by sDIVG.We first define a function r : Λ → ∂Ω, i.e., a map from the parametric domain Λ ⊂ R d−1 to the manifold ∂Ω ⊂ R d ; the Jacobian of this function is denoted by ∇r.As in the DIVG algorithm, let h denote a desired spacing function.Finally, given a node λ i ∈ Λ, the set of n DIVG-generated candidates for λ i can be written as It is important to note that the candidates η i,j all lie in the parametric domain.Our goal is to determine how far from λ i each candidate must lie.From the definition of r and h, we have for all j = 1, . . ., n.The candidates η i,j can be thought of as lying on some manifold around λ i .This allows us to rewrite η i,j as for some α i,j > 0 and unit vector s i,j .Here, α i,j must be determined by some appropriate transformation of h(r(λ i )), i.e., parametric distances are obtained by a transformation of the spacing function h specified on ∂Ω.We may now use (3) to Taylor expand r(η i,j ) as We can now use ( 4) in ( 2) to obtain the following expression for h(r(λ i )) in terms of α i,j : This in turn allows us to express α i,j as We can now use ( 6) within ( 1) to obtain an expression for the candidate set C i purely in terms of the spacing function h and the parametric transformation r.Thus, we have where S i is set of n random unit vectors on a unit ball.All other steps are identical to DIVG algorithm, albeit in the parametric domain Λ.The final set of points on ∂Ω is then obtained by evaluating the function r at these parametric samples; a schematic of this is shown in Figure 1 (right).A full description of the sDIVG algorithm can be found in [7].

NURBS-DIVG
In principle, sDIVG can be used with any map r : Λ → ∂Ω.NURBS-DIVG, however, is a specialization of sDIVG to surfaces comprised of NURBS patches (collections of non-overlapping and abutting NURBS).We first describe the use of NURBS for generating the surface representation r, then discuss the generalization to a surface containing multiple patches.

An overview of NURBS surfaces
To define a NURBS representation, it is useful to first define the corresponding B-spline basis in onedimension.Given a sequence of nondecreasing real numbers T = {t 0 , t 1 , . . ., t k } called the knot vector, the degree-p B-spline basis functions N i,p (u) are defined recursively as [27] N We can now use these basis function to define a NURBS curve in R d .Given a knot vector of the form  NURBS curve is defined as [27] In practice, it is convenient to evaluate s(u) ⊂ R d as a B-spline curve in R d+1 , and then project that curve down to R d .For more details, see [27].We evaluate the B-spline curve in a numerically stable and efficient fashion using the de Boor algorithm [4], which is itself a generalization of the well-known de Casteljau algorithm for Bezier curves [8].It is also important to note that derivatives of the NURBS curves with respect to u are needed for computing surface quantities (tangents and normals, for instance, or curvature).Fortunately, it is well-known that these parametric derivatives are also NURBS curves and can therefore also be evaluated using the de Boor algorithm [27].We adopt this approach in this work.Finally, the map r which represents a surface of co-dimension one in R d can be represented as a NURBS surface; this is merely a tensor-product of two NURBS curves s(u) and s(v), where u and v are two independent parameters in R d−1 .This allows all NURBS surface operations to be computed via the de Boor algorithm applied to each parametric dimension.The sDIVG algorithm can then be used to sample this surface as desired, which in turn allows for node generation in the interior of the domain via DIVG.

Sampling surfaces consisting of multiple NURBS patches
Practical CAD models typically consist of multiple non-overlapping and abutting NURBS surface patches.A NURBS patch meets another NURBS patch at a NURBS curve (the boundary curves of the respective patches).While degenerate situations can easily arise (such as two NURBS patches intersecting at a single point), we concern ourselves in this work with patches that intersect in NURBS curves.Some examples of node sets generated by NURBS-DIVG on CAD surfaces consisting of multiple NURBS patches are shown in Figure 2. We now describe the next piece of the NURBS-DIVG algorithm: extending sDIVG to discretize a CAD model that consists of several NURBS patches ∂Ω i .We proceed as follows: 1. We use sDIVG to populate patch boundaries ∂(∂Ω i ) with a set of nodes.Recall that these patch boundaries are NURBS curves.
Figure 3: Illustration of positioning nodes on a deformed sphere made of five NURBS patches.In the first step, the boundary of the first patch is filled (first), followed by filling of that patch interior (second).Once the first patch is processed, the boundary of the second patch is discretized (third); this process is repeated until all patches are fully populated with nodes (fourth).
2. We then use these generated nodes as seed nodes within another sDIVG run, this time to fill the NURBS surface patches ∂Ω i enclosed by those patch boundaries.
To populate patch boundaries, the boundary NURBS representation obtained from any of the intersecting patches can be used.But, to use nodes from patch boundaries as seed nodes in sDIVG for populating surface patches, the corresponding node from the patch's parametric domain Λ is required.While it is possible to determine the map Ω → Λ through a nonlinear solve, we found it simpler to simply populate the patch boundaries twice, once from each of the NURBS representations obtained from intersecting patches.This produces two sets of seed nodes (one corresponding to each patch), but only the set from one of the representations is used (it does not matter which one, since both node sets are of similar quality).The full process is illustrated in Figure 3.
For a given CAD model consisting of a union of NURBS patches and a desired node spacing h, it is possible that the smallest dimension of the patch becomes comparable to (or even smaller) than h.We now analyze the behavior of sDIVG in this regime.To do so, we construct simple models comprising of Bezier surfaces (NURBS with constant weights); to emulate the existence of multiple patches, we simply subdivide the Bezier surfaces to obtain patches.In the 3D case, each subdivision is performed in a different direction to ensure patches of similar size.The resulting surface, now a union of non-overlapping and abutting NURBS (Bezier) patches, was then discretized with the NURBS-DIVG algorithm using a uniform spacing of h = 10 −4 .We then assessed the quality of the resulting node sets on those patches using the normalized local regularity metric d i defined in Section 3. The models and this metric are shown in Figure 4. Figure 4 shows that NURBS-DIVG works as expected when h is considerably smaller than the patch size.However, once the patch size becomes comparable to the h, the NURBS-DIVG algorithm rejects all nodes except those on the boundaries of the patch, as there is not enough space on the patch itself for additional nodes. 1 In all discussions that follow, we restrict ourselves to the first and most natural regime where h is considerably smaller than the patch size.

NURBS-DIVG in the interior of CAD objects
As our goal is to generate node sets suitable for meshless numerical analysis, it is vital for the NURBS-DIVG algorithm to be able to generate node sets in the interior of volumes whose boundaries are CAD models, in turn defined as a union NURB patches.While it may appear that the original DIVG algorithm is already well-suited to this task, we encountered a problem of nodes "escaping" the domain interior when DIVG was applied naively in the CAD setting.We now explain this problem and the NURBS-DIVG solution more clearly.To discretize the interior of CAD objects, we must accurately determine whether a particular node lies inside or outside the model.The choice of boundary representation can greatly affect the technique used for such an inside/outside test.For instance, if the domain boundary is modeled as an implicit surface of the form f (x) = 0, a node x k is inside if f (x k ) < 0 (up to some tolerance).However, in the case where the domain boundary is modeled as a parametric surface or a collection of parametric patches (as in this work), the analogous approach would instead solve a nonlinear system of equations to find the parameter values corresponding to x k to test if x k is inside.A simpler approach used in recent work has been to simply find the closest point from the boundary discretization p (with the given spacing h) to x k , and use its unit outward normal to decide if x k is inside the domain.More concisely, if n is the unit outward normal vector at p, x k is inside the domain Ω when n This is the approach used by the DIVG algorithm (and many others).However, in our experience, this  does not work well for complex geometries with sharp edges and concavities.For an illustration, see Figure 5 (left); we see nodes marked as "interior" nodes that are visually outside the convex hull of the boundary nodes.
An investigation revealed that a relatively coarse sampling of a patch near its boundary NURBS curves could result in the closest point p and its normal vector n being a bad approximation of the actual closest point and its normal on the domain Ω, thereby resulting in x k being erroneously flagged as inside Ω.This problem is especially common on patch boundaries, where normal vectors do not vary smoothly.NURBS-DIVG uses a simple solution: supersampling.More precisely, we use a secondary set of refined boundary nodes only for the inside check with a reduced spacing ĥ given by ĥ = h/τ, where τ > 1 is a factor that determines the extent of supersampling.Though this potentially requires τ to be tuned, this solution worked well in our tests with a minimal additional implementation complexity and computational overhead (see execution profiles for Poisson's equation in section 4). Figure 5 (right) shows the effect of setting τ = 2 in the same domain; nodes no longer "escape" the boundary.While this approach is particularly useful for boundary represented as a collection of NURBS patches, it is likely to be useful in any setting where the boundary has sharp changes in the derivative of the normal vector (or the node spacing h).In fact, intuitively, it seems that the greater the derivative (or the smaller the value of h), the greater the value of τ required to prevent nodes escaping.To confirm this intuition, we run a simple test both in 2D and 3D.In 2D, we define a parametric curve where s is a parameter controlling the complexity of the curve (s = 1 gives 2 legs, s = 1.5 gives 3 legs, and so forth).These "legs" create sharp changes in the derivative of the normal (notice r(t) is not smooth in t).
In 3D, we simply extrude this curve in the z direction to obtain a surface: Both test domains are depicted in Figure 6.We then plot the minimum value of τ required for a successful inside check as a function of s, the parameter that controls the number of legs, and the node spacing h.The results are shown in Figure 7.As expected, increasing the number of legs via s necessitates a greater degree of supersampling (τ min in the plots) in both 2D and 3D.However, if h is sufficiently small to begin with, smaller values of τ appear to suffice.In the tests presented in later sections, we selected h to be sufficiently small that τ = 2 sufficed.

Node quality
Although node quality in meshfree context is not as well understood as in mesh based methods, we can analyze local regularity by examining distance distributions to nearest neighbors.For each node p i with and, given a constant h, normalize it to h.
We used c = 2 for 2D and c = 3 for 3D analyses.For the analysis following models are selected • 2D duck with 8 patches, • sphere with 6 patches, • deformed sphere with 6 patches, all depicted in Figure 8.The distance distributions to closest neighbors for boundary nodes are presented in Figure 9 and in Figure 10 distributions for all nodes are shown.The quantitative statistics are presented in Table 1.It can be seen that the nodes are quite uniformly distributed as all distributions are condensed around 1. In general, the boundary nodes are less uniformly distributed than interior nodes.This is expected, since the Taylor expansion (of the first order) used in sDIVG is an approximation, which means that candidates are not generated exactly at distance h from the node being expanded.This can be improved by using higher order Taylor expansion terms, but in practice, first order expansion is good enough even in complex geometries, where derivatives diverge [7].Furthermore, in the 2D duck case, the distribution visually seems much better than in the 3D cases.This is a consequence of the candidate generation procedure, which is more optimal when the (parametric) domains are 1D.Another feature characteristic for 1D parametric domains is that outliers with distance to nearest neighbors slightly less than 2h are not uncommon.This happens at nodes, where the advancing fronts of the sDIVG algorithm meet and is rarely a problem in practice.Because of that, the standard deviation for duck case shown in Table 1 is of the same order of magnitude as the 3D cases.If 10 most extreme outliers are not taken into account, the standard deviation reduces by an order of magnitude.See [7] for a deeper analysis and possible solutions.
In the 3D cases, the distribution of nodes for the deformed sphere is slightly worse, which can be attributed to a greater complexity of the model.Nevertheless, the quality of generated nodes is of the same order as the nodes generated by pure DIVG [32] and sDIVG [7].Additionally, there are two quantities often considered as node quality measures, i.e., minimal distance between nodes or separation distance and maximal empty sphere radius within the domain [12,36].The separation distance is defined for set of nodes Ξ = {x 1 , . . ., x N } ⊂ Ω as Quantity r min,Ξ is determined by searching a closest neighbor for all nodes using a spatial search structure, such as a k-d tree.A r max,Ξ is estimated numerically by sampling Ω with higher node density and searching for the closest node among Ξ.Behaviour of normalized maximum empty sphere radius and separation distance for all three cases with respect to target nodal distance h is presented in Figures 11 and 12.In all cases, r max is relatively stable near an acceptable value of 1 and r min approaches the optimal value of 0.5 with decreasing h.This behaviour is consistent with previous results and the analytical bound for r min [32,7].

Solving PDEs on CAD geometry
In this section we focus on solving PDEs using the proposed NURBS-DIVG and RBF-FD, a popular variant of strong form local meshless methods.In the first step the CAD provided domain is populated with nodes using NURBS-DIVG that are in the second step used as discretization nodes for a RBF-FD method using polyharmonic basis functions augmented with monomials to approximate considered the PDE [2,14].
In following numerical experiments we used RBF-FD with augmentation up to order m ∈ {2, 4, 6} on n = 4 m+2 2 closest nodes in 2D and n = 4 m+3 3 in 3D to obtain the mesh-free approximations of the differential operators involved.The sparse system resulting from discretizing the PDE at hand is solved with BiCGSTAB using an ILUT preconditioner.Note that in all cases we used ghost nodes on boundaries with Neumann or Robin boundary conditions.Once the numerical solution û is obtained, we observe the convergence behaviour of the solution through error norms defined as

Poisson's equation
In Figure 14 we can see that for all three geometries the solution converges with the expected order of accuracy according to the order of augmenting monomials.Next, we assess the execution time of solving Poisson's equation with second order monomial augmentation.
In Figure 15 the execution times for all three geometries are broken down to core modules of the solution procedure.We separately measure execution time for generation of nodes in interior of the domain (DIVG including inside check with supersampling) and on the surface of the domain (sDIVG on multiple NURBS patches). 2We measure generation of stencils and computing stencil weights together as a RBF-FD part of the solution procedure.Separately we also measure sparse matrix assembly (which is negligible [5]) and its solution that with increasing number of nodes ultimately dominates the execution time [5].In 2D duck case the computational times of solving the system and filling the interior of the domain are of the same order as the number of nodes is still relatively small.However, in both 3D cases we can clearly see that solving system scales super-linearly and soon dominates the overall computational cost.The RBF-FD part as expected scales almost linearly (neglecting the O(N log N ) resulting from k-d tree in stencil selection).

Linear elasticity -Navier Cauchy equation
In previous section we established confidence in the presented solution procedure by obtaining expected convergence rates in solving Poisson's equation on three different geometries in 2D and 3D.In this section we apply NURBS-DIVG to a more realistic case from linear elasticity, governed by the Navier -Cauchy equation where u stands for the displacement vector, and Young's modulus E = 72.1 • 10 9 Pa and Poisson's ratio ν = 0.33 define material properties.The displacement and the stress tensor (σ) are related via Hooke's law with ε and I standing for strain and identity tensors.We observe a 3D gear object that is subjected to an external torque resulting in the tangential traction t 0 = 1 • 10 3 Pa on axis, while the gear teeth are blocked, i.e. displacement u = 0 m.Top and bottom surfaces are free, i.e. traction free boundary conditions apply The case is schematically presented in Figure 16 together with von Mises stress scatter plot.The stress is highest near the axis, where the force is applied, and gradually fades towards blocked gear teeth.Displacement and von Mises stress are further demonstrated in Figure 17 at z = 0 m cross section, where we can see how the gear is deformed due to the applied force.Presented results were computed using 41210 scattered nodes generated by NURBS-DIVG.

Transient heat transport
Last example is focused on dimensionless transient heat equation where T stands for temperature and q for the heat source.The goal is to solve heat transport within the duck model subjected to the Robin boundary conditions ∂T ∂t + T = 0 (35) and heat source within the domain q = 5e 10||x−x0|| (36) with x 0 = (0, 0, 0.2) and initial temperature set to 0 throughout the domain.Time marching is performed via implicit stepping using time step ∆t = 3 • 10 −4 and 3000 iterations to reach the steady state at t = 0.9 Figure 18 shows the temperature scatter plot computed with RBF-FD on 21956 nodes generated with the proposed NURBS-DIVG at two different times (first at the beginning of the simulation and second at the steady state).In Figure 19, the time evolution of the temperature at five control points P 1 − 5 is shown.Control point P 1 is located at the heat source, P 2 − 4 at most distant points from the source and P 5 asymmetric with respect to the y-axis.As one would expect, at the source the temperature rises immediately after the beginning of the simulation and also reaches the highest value, while the rise is a bit delayed and lower at the distant points that are closer to the boundary where the heat exchange with surroundings takes place.Once the heat heat exchange with surrounding matches the heat generation at the source, the system reaches steady-state.
Figure 18: Heat transport within a 3D duck.

Conclusions
In this paper, we presented a meshless algorithm, NURBS-DIVG, for generating quasi-uniform nodes on domains whose boundaries are defined by CAD models consisting of multiple NURBS patches.The NURBS-DIVG algorithm is able to deal with complex geometries with sharp edges and concavities, supports refinement, and can be generalized to higher dimensions.We also demonstrated that node layouts generated with NURBS-DIVG are of sufficiently high quality for meshless discretizations, first by directly assessing the quality of these node sets, then by

Figure 2 :
Figure 2: Node sets generated by NURBS-DIVG on a CAD Utah Teapot (left) and a CAD model of tiger (right).The Utah Teapot model is made of 32 patches and has 7031 nodes; the tiger has 124 patches and 4753 boundary nodes.

Figure 5 :
Figure 5: Demonstration of the supersampling approach for the inside/outside test in NURBS-DIVG.The figure on the left shows nodes generated by the naive test used in the DIVG algorithm, with nodes escaping the domain boundary.The figure on the right shows the nodes generated using boundary supersampling in NURBS-DIVG; all non-boundary nodes are enclosed within the volume defined by the boundary.

Figure 6 :
Figure 6: Shapes used to test the supersampling approach in NURBS-DIVG.

6 Figure 7 :
Figure 7: Minimal τ required to appropriately fill a model.

Figure 8 :
Figure 8: Geometries used in the node quality analysis.

Figure 10 :
Figure 10: Local regularity distributions for boundary and interior nodes.

Figure 11 :Figure 12 :
Figure 11: Minimal and maximal spacing on the boundary.

First, we solve
a Poisson's equation∇ 2 u = f(21)with known closed-form solutionu a (x, y) = sin π 100 x cos 2π 100 y , in 2D,(22)u a (x, y, z) = sin(πx) cos(2πy) sin(0.5πz), in 3D, (23)using mixed Neumann-Dirichlet boundary conditions u = u a , on Γ d ,(24)∂u ∂n = ∂u a ∂n , on Γ n , (25) where Γ d and Γ n stand for Dirichlet and Neumann boundaries.In all cases a domain is divided into two halves, where we apply a Dirichlet boundary condition to one half and a Neumann boundary condition to another.The numerical solution of the problem is presented in Figure 13.

Figure 13 :
Figure 13: Solution of Poisson's equation on all three considered geometries.

Figure 14 :
Figure 14: Poisson's equation solution error with respect to the number of nodes.

Figure 15 :
Figure 15: Execution times broken down to separate solution procedure modules.

Figure 16 :
Figure 16: Scheme of the linear elasticity example (left) accompanied with RBF-FD solution in terms of von Misses stress (right).

Figure 17 :
Figure 17: von Mises stress (left) and displacement magnitude (right) at z = 0 cross section accompanied with quiver plot of displacement field.
using RBF-FD to solve the Poisson equation with mixed Dirichlet-Neumann boundary conditions on different domains to high-order accuracy.Finally, we demonstrated NURBS-DIVG in conjunction with RBF-FD in tackling two more challenging test cases: first, the stress analysis of a gear subjected to an external force governed by the Navier-Cauchy equation; and second, a time-dependent heat transport problem inside a duck.This work advances the state of the art in fully-autonomous, meshless, isogeometric analysis.All algorithms presented in this work are implemented in C++ and included in our in-house open-source meshfree library Medusa [20, 33], see the Medusa wiki [21] for usage examples.The interface to CAD files was implemented via Open Cascade [25].

Figure 19 :
Figure 19: Time evolution of temperature at five control points.