Algebraic Point Projection for Immersed Boundary Analysis on Low Degree NURBS Curves and Surfaces

Abstract

Point projection is an important geometric need when boundaries described by parametric curves and surfaces are immersed in domains. In problems where an immersed parametric boundary evolves with time as in solidification or fracture analysis, the projection from a point in the domain to the boundary is necessary to determine the interaction of the moving boundary with the underlying domain approximation. Furthermore, during analysis, since the driving force behind interface evolution depends on locally computed curvatures and normals, it is ideal if the parametric entity is not approximated as piecewise-linear. To address this challenge, we present in this paper an algebraic procedure to project a point on to Non-uniform rational B-spline (NURBS) curves and surfaces. The developed technique utilizes the resultant theory to construct implicit forms of parametric Bézier patches, level sets of which are termed algebraic level sets (ALS). Boolean compositions of the algebraic level sets are carried out using the theory of R-functions. The algebraic level sets and their gradients at a given point on the domain are then used to project the point onto the immersed boundary. Beginning with a first-order algorithm, sequentially refined procedures culminating in a second-order projection algorithm are described for NURBS curves and surfaces. Examples are presented to illustrate the efficiency and robustness of the developed method. More importantly, the method is shown to be robust and able to generate valid solutions even for curves and surfaces with high local curvature or ${\mathcal{G}}_0$ continuity—problems where the Newton–Raphson method fails due to discontinuity in the projected points or because the numerical iterations fail to converge to a solution, respectively.

1. Introduction

Given a test point and a parametric entity (curve or surface), the generalized point projection problem is to find the closest point (footpoint) on the entity as well as the corresponding parameter value. Since the footpoint is the closest point on the curve or surface, the line connecting the test point to the footpoint is normal to the curve or the surface [1]:

$$\begin{array}{c}\hfill g\left(u\right)={\mathbf{C}}^{\prime }\left(u\right)·(\mathbf{C}\left(u\right)-\mathbf{P})=0\end{array}$$

(1)

Given a parametric curve or surface entity $\mathbf{C}\left(u\right)\in {\mathbb{R}}^n$ (u is treated as a vector when the entity is a surface), the Euclidean distance function $d_E\left(\mathbf{x}\right)$ is defined as the shortest distance from physical test point $\mathbf{x}$ to $\mathbf{C}\left(u\right)$ given by:

$$\begin{array}{c}\hfill d_E\left(\mathbf{x}\right)=\underset{u}{inf}\parallel \mathbf{x}-\mathbf{C}\left(u\right)\parallel \end{array}$$

(2)

where $\mathbf{C}\left(u\right)$ is a physical point on the curve or surface of interest. The distance function $d_E\left(\mathbf{x}\right)$ is continuous for all $\mathbf{x}\in {\mathbb{R}}^n$ and differentiable almost everywhere. The “footpoint” of projection in parametric space $u_f$ is defined as:

$$\begin{array}{c}\hfill u_f\left(\mathbf{x}\right)=\underset{u\in [a,b]}{arg\; min}\parallel \mathbf{x}-\mathbf{C}\left(u\right)\parallel \end{array}$$

(3)

where $[a,b]$ is the parameter range. In general, $u_f\left(\mathbf{x}\right)$ may be non-unique, discontinuous, or non-existent, as illustrated in Figure 1. The footpoint of a test point near the curve or surface segment with high local curvature can be non-unique, leading to discontinuity of point projection process as illustrated in Figure 1a. The non-existence is illustrated in Figure 1b, and occurs around points where $\mathbf{C}\left(u\right)$ has only ${\mathcal{G}}^0$ continuity.

This problem is of importance in geometric modeling. For instance, while fitting a curve or surface to sampled data, one may need to compute corresponding parameter values and errors at data points since the error is the distance between the data point and the fitting curve or surface [2].

Point projection also plays an important role in computer aided engineering (CAE), especially when boundaries are immersed into the domain and evolved. Such immersed boundary analysis [3] uses a non-conforming mesh to significantly reduce computational cost required for mesh generation as the boundaries evolve. Often, the immersed boundaries are represented as parametric splines, and, more recently, in isogeometric analysis, the underlying domain is also approximated by parametric splines. Isogeometric analysis (IGA) [4,5] is aimed at building behavioral approximations isoparametrically on a spline geometry, most commonly on NURBS curves and surfaces. The goal is to eliminate the need to mesh the geometry for analysis and to ensure the exactness of geometry to the CAD model during analysis. Tambat and Subbarayan [6] developed an Enriched Isogeometric Analysis (EIGA) for immersed boundary problems in which both the domain as well as the enrichments are described by NURBS entities, which are then blended to describe the enriched approximation.

In any immersed boundary problem solution, capturing the interaction of the field approximation defined on the immersed (explicitly defined) boundary with the approximations on the enriched domain requires one to determine the nearest point on the boundary from any given point in the underlying domain. This projection from the spatial point to the boundary is necessary to compute the influence of the domain approximation on those approximations defined on the boundary (see Figure 2). For example, in solutions to mechanical contact problems [7,8], point projection is needed to define the normal gap and tangential slip between two bodies. In fluid–structure interaction (FSI) problems, point projection is required to transfer kinematic and traction data between non-matching fluid–structure interface [9]. In the enriched isogeometric analysis mentioned above, point projection is used to enrich the base approximations with those on lower-dimensional geometrical features such as crack surfaces and phase boundaries, enabling simulations of fracture propagation [6,10] and solidification [11]. A fast and robust point projection method is critical to efficiently solving these problems.

The rest of the paper is organized as follows. In Section 2, a review of the literature pertaining to the point projection problem is carried out. In Section 3, the algebraic estimation of distance from a low-degree NURBS curve or surface is reviewed. In Section 4.1, a detailed algorithm for two-dimensional algebraic point projection is presented followed in Section 4.2 by its extension to three-dimensional NURBS surfaces. Several examples are provided in Section 5 to validate the developed algorithm. The paper is concluded with remarks in Section 6.

2. Literature Review

The use of Newton–Raphson (NR) iterations for solving Equation (1) is well established at this time. These iterative methods mainly consist of two steps:

• Seek an initial point or segment.
• Iterate by Newton–Raphson scheme until convergence.

The robustness and the efficiency of Newton–Raphson scheme depends significantly on the initial guess. Therefore, to assure convergence of the second step, careful selection of initial guess is needed. In addition, if the NURBS entity has only ${\mathcal{G}}_0$ continuity at some local point, the derivative based Newton–Raphson scheme would fail to converge to such a point.

To assure robustness of the iterations, a significant focus of the existing literature is on eliminating portions of the curve or surface where the solution cannot lie. Piegl and Tiller [2] developed a non-iterative, heuristic algorithm where a NURBS surface was decomposed into quadrilaterals and test points were projected onto the closest quadrilateral. Ma and Hewitt [12] described a search for the initial guess of the footpoint by recursively sub-dividing the Bézier segment associated with a valid control polygon. However, using the control polygon to eliminate segments of the curve may exclude the correct solution [13]. Selimovic [14] improved Ma and Hewitt’s method using selective sub-division of the NURBS curve (surface) based on distance to the end (corner) point of the entity. Chen et al. [15] pointed out the need for all control points to lie on different sides of a hyperplane in Selimovic’s algorithm and proposed a circular clipping method with a sufficiency condition for a curve to lie outside an elimination circle. Other iterative methods in the physical space have also been proposed for point projection including the geometric first-order iterative [16,17,18] and geometric second-order iterative [19,20,21] methods.

In this paper, a robust and efficient point projection technique for low degree two-dimensional (2D) NURBS curves and three-dimensional (3D) NURBS surfaces is developed. The proposed technique preserves and operates directly on the parametric description of the NURBS curve or surface. Therefore, the technique gives a projected point directly on the curve or surface when query points lie on the parametric curve or surface. In addition, the technique is robust for curves and surfaces with high local curvature or ${\mathcal{G}}_0$ continuity compared to techniques that rely on derivatives. A detailed comparison of existing methods in the literature relative to the proposed method is listed in

Table 1. Comparison of methods used for point projection in literature

Description	Reference	Algebraic	Initial Guess Dependent?	Efficiency	Accuracy	Smoothness
Subdivision method	[22]	Yes	No	High	Medium	No
Subdivision method + Newton–Raphson method	[12,14,15]	No	Yes	High	High	No
Geometric iteration method	[16,17,18,19,20,21]	No	Yes	High	High	No
Proposed method		Yes	No	High	Medium	Yes

.

3. Background on Algebraic Level Sets

In addition to point projection, blending behaviors on immersed boundaries with domain also requires distance estimates from the boundary to a point in the domain. This is because the behavioral influence of the immersed boundaries decays monotonically with distance. While traditionally distances are estimated using Newton–Raphson iterations, recently, Upreti et al. [23,24] developed algebraic procedures for efficient computation of distance estimates from curves and surfaces. The method developed by Upreti et al. relies on converting the parametric NURBS geometry to its implicit form using the Dixon resultant, and constructing level sets on the implicit form of the geometry. Upreti et al. termed these level sets as algebraic level sets. The construction of the algebraic level sets requires one to decompose the NURBS entity into constituent Bézier patches and later to blend the level sets constructed on Bézier patches using R-functions. The algebraic level sets provide monotonic measures of distance that are accurate to exact distance near the boundary. While the algebraic level sets are approximate measures of distance, they are sufficient to capture the interaction of domain approximations with those on the immersed boundary during analysis. The algebraic level sets have the following properties:

• Accurate measure of distance locally near the curve/surface
• Monotonic function of Euclidean (exact) distance
• Sufficiently smooth for engineering applications
• Efficiently obtained without numerical iterations for points close to the curve/surface

For the sake of completeness, we briefly review the computation of algebraic level sets and illustrate the procedure through simple examples.

3.1. Implicitization of a Parametric Curve

Given a rational parametric curve $\mathbf{C}\left(X\right(u),Y(u),W(u\left)\right)$ of degree p with $x=\frac{X\left(u\right)}{W\left(u\right)},y=\frac{Y\left(u\right)}{W\left(u\right)}$, one can construct two auxiliary polynomials:

$$\begin{array}{cc}\hfill g_1(x,u)& =W\left(u\right)x-X\left(u\right)=0\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill g_2(y,u)& =W\left(u\right)y-Y\left(u\right)=0\hfill \end{array}$$

(4b)

The above polynomial equations can be rearranged in descending power of u as follows:

$$\begin{array}{cc}\hfill g_1\left(u\right)& =a_p{u}^p+a_{p-1}{u}^{p-1}+\cdots +a_{1}u+a_0\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill g_2\left(u\right)& =b_p{u}^p+b_{p-1}{u}^{p-1}+\cdots +b_{1}u+b_0\hfill \end{array}$$

(5b)

From the above, the following resultant system may be obtained through algebraic manipulations [25]:

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}\left(a_p{b}_{p-1}\right)& \cdots & \left(a_p{b}_0\right)\\ ⋮& \ddots & ⋮\\ \left(a_p{b}_0\right)& \cdots & \left(a_1{b}_0\right)\end{array}\right]\left(\begin{array}{c}{u}^{p-1}\\ u^{p-2}\\ ⋮\\ 1\end{array}\right)={\left[{\mathbf{M}}^B\right]}_{p\times p}\left(\begin{array}{c}{u}^{p-1}\\ u^{p-2}\\ ⋮\\ 1\end{array}\right)=0\end{array}$$

(6)

where $\left(a_i{b}_j\right)=a_i{b}_j-a_j{b}_i$, ${\mathbf{M}}^B$ is Bezout matrix and is a function of x and y having the following important property:

$$\begin{array}{c}\hfill {\mathbf{M}}^B(x,y)={\mathbf{M}}_x^{B}x+{\mathbf{M}}_y^{B}y+{\mathbf{M}}_w^B\end{array}$$

(7)

where ${\mathbf{M}}_x^B,{\mathbf{M}}_y^B$, and ${\mathbf{M}}_w^B$ depend on control point coordinates and weights. Therefore, these matrices can be pre-computed for a given rational parametric curve and re-used given any new physical point $\mathbf{x}$. The determinant, $\det \left({\mathbf{M}}^B\left(\mathbf{x}\right)\right)$, is defined as the Bezout resultant. Since all the allowable parameter values u for curve $\mathbf{C}\left(X\right(u),Y(u),W(u\left)\right)$ are roots of Equation (6), $\det \left({\mathbf{M}}^B\left(\mathbf{x}\right)\right)=0$ gives the equation of the implicitized curve. Thus, the algebraic level sets corresponding to a rational parametric curve (e.g., Bézier curve) are given by:

$$\begin{array}{c}\hfill \mathsf{\Gamma }\left(\mathbf{x}\right)=\det \left({\mathbf{M}}^B\left(\mathbf{x}\right)\right)\end{array}$$

(8)

An example of algebraic level sets is shown in Figure 3.

3.2. Boolean Operations by R-Functions

As observed in Figure 3, the direct implicitization extends the parametric curve beyond its end points, and yields an invalid distance measure in the extended region. Therefore, it is desirable to trim the curve $\mathbf{C}\left(X\right(u),Y(u),W(u\left)\right)$ within its parameter range $u\in [a,b]$. In related prior work, Biswas and Shapiro [26] constructed an approximate distance from a line segment as:

$$\begin{array}{c}\hfill g=\sqrt{{\mathsf{\Gamma }}^2+\frac{{\left(\right|\varphi |-\varphi )}^2}{4}}\end{array}$$

(9)

with $\mathsf{\Gamma }$ being the distance in the normal direction from a piecewise linear approximation of the curve, and $\varphi$ being distance to the region formed by a circle circumscribing the line that is positive inside and negative outside. This form yields a smooth distance function across the boundary $\varphi =0$. Upreti et al. [23] extended the above idea by carrying out Boolean operations on fields obtained on (individual curve segments of) an arbitrarily shaped parametric curve and an enclosing convex region using R-functions (Figure 4). The R-functions [27,28] enable a smooth and purely algebraic Boolean operation, and result in a continuous distance measure. Two specific R-functions used in this study are:

• R-conjunction, equivalent to Boolean intersection:

  $$\begin{array}{c}\hfill g_1\wedge g_2=g_1+g_2-\sqrt{g_1^2+g_2^2}\end{array}$$

  (10)
• R-disjunction, equivalent to Boolean union:

  $$\begin{array}{c}\hfill g_1\vee g_2=g_1+g_2+\sqrt{g_1^2+g_2^2}\end{array}$$

  (11)

A well known property of Bézier and NURBS geometry is the convex hull property, which assures that the curve/surface is contained within its convex hull constructed using the control points. Upreti et al. [23] used the convex hull property of Bézier and NURBS curves to provide a natural convex region bounded by control points for curve trimming. Assume that the ith hyper-plane of the convex hull is expressed as:

$$\begin{array}{c}\hfill h_i\left(\mathbf{x}\right)={\mathbf{n}}_i·\mathbf{x}+b_i=0\end{array}$$

(12)

where $\mathbf{x}\in {\mathbb{R}}^n$ is a spatial point, ${\mathbf{n}}_i\in {\mathbb{R}}^n$ is inner normal, and $b_i\in \mathbb{R}$ is offset. Thus, the exact distance from any point $\mathbf{x}$ to the hyper-plane $h_i\left(\mathbf{x}\right)=0$ is $h_i\left(\mathbf{x}\right)$. The function $\varphi$ can be obtained by applying R-conjunction operation of Equation (10) to all $h_i\left(\mathbf{x}\right),i=1,2,\dots ,n$. An example of a cubic Bézier curve is shown in Figure 5.

3.3. Normalization and Composition of Algebraic Level Sets

The aforementioned procedure generates a monotonic and continuous distance measure for a basic parametric curve such as a Bézier curve. Piecewise polynomial curves such as NURBS curves, on the other hand, require decomposition to Bézier segments and composition of algebraic level sets of the obtained segments. Further, normalization for individual level sets is desired to yield a monotonically varying composed field. Considering a physical footpoint ${\mathbf{x}}_f$, one can approximate $\mathsf{\Gamma }\left({\mathbf{x}}_f\right)$ to a first order using Taylor expansion:

$$\begin{array}{c}\hfill \mathsf{\Gamma }\left({\mathbf{x}}_f\right)=\mathsf{\Gamma }\left(\mathbf{x}\right)-\frac{\partial \mathsf{\Gamma }\left(\mathbf{x}\right)}{\partial \mathbf{n}}d\end{array}$$

(13)

Since the resultant has exact zero set on a parametric curve, i.e., $\mathsf{\Gamma }\left({\mathbf{x}}_f\right)=0$, one can derive a distance in the normal direction as follows:

$$\begin{array}{cc}\hfill d& =\frac{\mathsf{\Gamma }\left(\mathbf{x}\right)}{\frac{\partial \mathsf{\Gamma }\left(\mathbf{x}\right)}{\partial \mathbf{n}}}=\frac{\mathsf{\Gamma }}{\parallel \nabla \mathsf{\Gamma }\parallel }\hfill \end{array}$$

(14)

After obtaining normalized algebraic level sets for each decomposed Bézier segment, one can compose them using R-conjunction operation (Equation (10)) and thereby generate the desired algebraic level sets. As demonstrated in [23], the R-conjunction operation preserves the normalization of individual Bézier segments. However, an implicitized curve obtained from a Bézier curve of degree p may have as many as $\frac{1}{2}(p-1)(p-2)$ self-intersections or double points [29]. Any double points inside the convex hull will affect the algebraic level sets construction, and therefore need to be moved out by sub-divisions of the Bézier curve. The algorithm to carryout this process is discussed in Reference [23]. Thus, for practical reasons of avoiding more than one double point while enabling sufficient generality in modeling complex geometries, the methodology is restricted to low degree NURBS curves $(p\le 3)$. Figure 6 shows an example of algebraic level sets of an open curve containing two points with ${\mathcal{G}}^0$ continuity.

3.4. Extension to NURBS Surface

The algebraic level sets construction can be extended to three-dimensional NURBS surfaces in a straightforward manner by implicitizing the rational parametric surface with the Dixon resultant [25]. Given a rational parametric surface $\mathbf{S}\left(X\right(u,v),Y(u,v),Z(u,v),W(u,v\left)\right)$ of degree $p\times q$ with $x=\frac{X(u,v)}{W(u,v)},y=\frac{Y(u,v)}{W(u,v)}$ and $z=\frac{Z(u,v)}{W(u,v)}$, three auxiliary polynomials can be formed as follows:

$$\begin{array}{cc}\hfill g_1(x,u,v)& =W(u,v)x-X(u,v)=0\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill g_2(y,u,v)& =W(u,v)y-Y(u,v)=0\hfill \end{array}$$

(15b)

$$\begin{array}{cc}\hfill g_3(z,u,v)& =W(u,v)z-X(u,v)=0\hfill \end{array}$$

(15c)

As before, using algebraic elimination theory, one can derive the corresponding resultant system for surface $\mathbf{S}$:

$$\begin{array}{c}\hfill {\left[{\mathbf{M}}^D\right]}_{2pq\times 2pq}{\left(\begin{array}{ccccccccc}1& u& \cdots & u^{p-1}& \cdots & v^{2q-1}& u{v}^{2q-1}& \cdots & u^{p-1}{v}^{2q-1}\end{array}\right)}^T=0\end{array}$$

(16)

where the vector is indexed lexicographically. ${\mathbf{M}}^D$ is the Dixon matrix, which also possesses a property analogous to Equation (7) of linearity with respect to x, y, and z:

$$\begin{array}{c}\hfill {\mathbf{M}}^D\left(\mathbf{x}\right)={\mathbf{M}}_x^{D}x+{\mathbf{M}}_y^{D}y+{\mathbf{M}}_z^{D}z+{\mathbf{M}}_w^D\end{array}$$

(17)

where, as before, ${\mathbf{M}}_x^D,{\mathbf{M}}_y^D,{\mathbf{M}}_z^D$, and ${\mathbf{M}}_w^D$ depend on control point coordinates and weights. The determinant of the Dixon matrix is the Dixon resultant:

$$\begin{array}{c}\hfill \mathsf{\Gamma }\left(\mathbf{x}\right)=\det \left({\mathbf{M}}^D\left(\mathbf{x}\right)\right)\end{array}$$

(18)

An example of the algebraic level sets from a free surface is illustrated in Figure 7. The pseudo-code in Algorithm 1 shows the generic steps in algebraic level sets computation for NURBS curves and surfaces. Both NURBS curves and surfaces are denoted by $\mathbf{C}\left(\mathbf{u}\right)$ here for notational convenience, with the implicit understanding that $\mathbf{u}=\left(u\right)$ for curves and $\mathbf{u}=(u,v)$ for surfaces.

Algorithm 1 Algebraic level sets algorithm.

Input: NURBS curve or surface $\mathrm{C}\left(\mathrm{u}\right)$ and given test point $\mathrm{x}$ Output: Algebraic distance measure d from $\mathrm{x}$ to the NURBS entity $\mathrm{C}\left(\mathrm{u}\right)$ 1: functionAlgebraic_Distance($\mathrm{C}$, $\mathrm{x}$) 2:     $\mathcal{B}\left(\mathrm{C}\right)$← Split NURBS entity $\mathrm{C}$ into a Bézier set with segments $B_i,i=1,2,\cdots ,n$ 3:     for $i\leftarrow 1,n$ do          ▹ Loops are independent and parallelizable 4:         $h_i$← Create convex hull for $B_i\in \mathcal{B}\left(\mathrm{C}\right)$ 5:         ${\mathsf{\Gamma }}_i$← Compute the Bezout or Dixon resultant using Equation (8) or Equation (18), respectively 6:         ${\mathsf{\Gamma }}_i$←$\frac{{\mathsf{\Gamma }}_i}{\parallel \nabla {\mathsf{\Gamma }}_i\parallel }$       ▹ Normalization of the resultant using Equation (14) 7:         $d_i$← Carryout Boolean union of distance fields of $h_i$ obtained using Equation (9) with ${\mathsf{\Gamma }}_i$ 8:     end for 9:     d← Carryout Boolean intersection of individual level sets $d_i,i=1,2,\cdots ,n$ using Equation (10) 10: end function

3.5. Time Complexity of the Algebraic Level Sets Algorithm

The algorithmic complexity of Algorithm 1 is given in [23] and reproduced here. Consider a NURBS curve of degree p. Further, assume that the NURBS curve decomposes into n Bézier curves. This decomposition process can be carried out in $O\left(n\right)$ time. For each step in constructing the algebraic level sets of the Bézier curve, the time complexity is a function of its degree p and is given in

Table 2. Time complexity of each step in Algorithm 1 for computing the algebraic level sets for a Bézier segment. Time complexities are listed for Bézier curves of degree p and Bézier surfaces of degree $p\times q$.

Step	Time Complexity
	Curve	Surface
Convex hull construction	$O(plog(p\left)\right)$	$O(pqlog(pq\left)\right)$
Distance field of convex hull	$O\left(p\right)$	$O\left(pq\right)$
Computing Dixon resultant	$O\left(p^3\right)$	$O\left(p^3{q}^3\right)$
Normalization of resultant	$O\left(p^3\right)$	$O\left(p^3{q}^3\right)$
Trimming operation	$O\left(1\right)$	$O\left(1\right)$

. Since each Bézier curve comprises of $p+1$ control points, the construction of its convex hull requires $O(plogp)$ time. Since the convex hull contains at most $p+1$ edges, computing the normalized distance field for the hull requires $O\left(p\right)$ time. Finding the Bezout/Dixon resultant involves evaluating the determinant of the Bezout matrix of size $p\times p$ and costs $O\left(p^3\right)$ time. Normalization of this resultant requires the gradient of the resultant, which in turn requires solving a linear system (see Equation (21)) and costs $O\left(p^3\right)$ time. The trimming operation only evaluates Equation (9) and can be carried out in $O\left(1\right)$ time. Thus, constructing algebraic level sets for a Bézier curve segment requires $O\left(p^3\right)$ time. These level sets constructed on Bézier curve segments can then be combined to construct level sets for the NURBS curve in $O\left(n{p}^3\right)$ time. The computational time complexity for a NURBS surface of degree $p\times q$ may be obtained analogously. For a Bézier segment of such a NURBS surface, the convex hull contains at most $(p+1)(q+1)$ points and edges, and the Dixon matrix used for the resultant is of size $2pq\times 2pq$. From the step-by-step time complexities listed in Table 2, construction of algebraic level sets for a Bézier surface segment requires $O\left(p^3{q}^3\right)$ time, while that of the NURBS surface requires $O\left(n{p}^3{q}^3\right)$ time.

4. Methodology of Algebraic Point Projection

In this section, we describe the developed procedure for algebraic point projection. The development of the algorithmic procedure is initially motivated using one-parameter NURBS curves and later extended to two-parameter NURBS surfaces.

4.1. Algebraic Point Projection for a NURBS Curve

As illustrated using Figure 1, iterative numerical solution for point projection may lead to a discontinuity in the projected point or may miss the correct solution. Hence, we develop in this section a purely algebraic point projection algorithm with the following properties:

• Exact at any point on the curve or surface, i.e., exact point inversion
• Controllable accuracy when projected from points near the curve or surface
• Efficient, non-iterative, and non-recursive solution procedure
• Footpoints are continuous even near curve segments with high curvature
• Valid solutions even when projected onto curves with only ${\mathcal{G}}^0$ continuity

The present method consists of two steps: estimation of the footpoint in physical space and parametric inversion using the resultant matrix.

4.1.1. First Order Algebraic Point Projection in Physical Space

From Equation (14), the gradient of the normalized approximate distance function to a Bézier segment is derived as:

$$\begin{array}{cc}\hfill \nabla d& =\left(\mathbf{I}-\frac{\mathsf{\Gamma }}{{\parallel \nabla \mathsf{\Gamma }\parallel }^2}\mathbf{H}\right)\frac{\nabla \mathsf{\Gamma }}{\parallel \nabla \mathsf{\Gamma }\parallel }\hfill \end{array}$$

(19)

where $\mathbf{I}$ is the identity matrix and $\mathbf{H}$ is the Hessian of function $\mathsf{\Gamma }\left(\mathbf{x}\right)$. Using the above gradient, the physical footpoint ${\mathbf{x}}_f$ can now be approximately located as:

$$\begin{array}{c}\hfill {\mathbf{x}}_f=\mathbf{x}-d\frac{\nabla d}{\parallel \nabla d\parallel }\end{array}$$

(20)

To calculate d and $\nabla d$ using Equations (14) and (19), it would appear at first glance as though $\mathsf{\Gamma },\nabla \mathsf{\Gamma }$, and $\mathbf{H}$ need to be evaluated for every test point. However, the following derivation as well as procedural detail show that d and $\nabla d$ can be computed efficiently without explicitly calculating $\mathsf{\Gamma },\nabla \mathsf{\Gamma }$ or $\mathbf{H}$. One can express $\nabla \mathsf{\Gamma }$ and $\mathbf{H}$ in terms of Bezout matrix ${\mathbf{M}}^B$ and its components ${\mathbf{M}}_x^B$ and ${\mathbf{M}}_y^B$ (the superscript B is dropped in the following for ease of reading):

$$\begin{array}{cc}\hfill \nabla \mathsf{\Gamma }& =\left|\mathbf{M}\right|\left[\begin{array}{c}\mathrm{tr}\left({\mathbf{M}}^{-1}\frac{\partial }{\partial x}\mathbf{M}\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}\frac{\partial }{\partial y}\mathbf{M}\right)\end{array}\right]=\left|\mathbf{M}\right|\left[\begin{array}{c}\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\end{array}\right]=\mathsf{\Gamma }\tilde{\mathbf{g}}\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \mathbf{H}& =\left|\mathbf{M}\right|\left[\begin{array}{cc}{\mathrm{tr}}^2\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)& \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\\ -\mathrm{tr}\left({\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)}^2\right)& -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)& {\mathrm{tr}}^2\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\\ -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\right)& -\mathrm{tr}\left({\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)}^2\right)\end{array}\right]=\mathsf{\Gamma }\tilde{\mathbf{H}}\hfill \end{array}$$

(22)

where $\tilde{\mathbf{g}}$ and $\tilde{\mathbf{H}}$ are the vector/matrix multiplying $\left|\mathbf{M}\right|$ in the above equations, respectively. Substituting Equations (21) and (22) back into Equations (14) and (19), one obtains:

$$\begin{array}{cc}\hfill d& =\frac{1}{\parallel \tilde{\mathbf{g}}\parallel }\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \nabla d& =\frac{\tilde{\mathbf{g}}}{\parallel \tilde{\mathbf{g}}\parallel }-\frac{\tilde{\mathbf{H}}\tilde{\mathbf{g}}}{{\parallel \tilde{\mathbf{g}}\parallel }^3}\hfill \end{array}$$

(24)

The efficiency of algebraic point projection in two-dimensional physical space is summarized as follows:

• Component matrices ${\mathbf{M}}_x$, ${\mathbf{M}}_y$, and ${\mathbf{M}}_w$ are constant for a given rational parametric curve. Therefore, they can be pre-computed and repeatedly used at a point $\mathbf{x}$.
• Only matrix $\mathbf{M}$ needs to be factorized, and the procedure extensively reuses the products ${\mathbf{M}}^{-1}{\mathbf{M}}_x$ and ${\mathbf{M}}^{-1}{\mathbf{M}}_y$ when computing $\tilde{\mathbf{H}}$ and $\tilde{\mathbf{g}}$.
• For a Bezout matrix $\mathbf{M}$ of size $p\times p$, where p is the degree of the rational parametric curve, the typical computational cost is low since p is usually small in engineering applications.

4.1.2. Second Order Algebraic Point Projection in Physical Space

For test points near the curve, the first-order algebraic point projection method described in Section 4.1.1 performs well; however, as shown through the example in Figure 8, when points are farther away, the projection is not accurate. The example is of a quadratic Bézier curve with test points on a horizontal line. The failure of first-order algebraic point projection is because the distance function derived in Equation (14) is essentially a first order approximation using Taylor expansion. To improve accuracy and make algebraic point projection possible for points farther from the curve, we present next a second-order algebraic point projection algorithm in the following. While the second-order algorithm is also approximate, in practice, it is sufficient for enriched immersed boundary analysis since only quadrature points near the immersed boundary needs to be projected on to the surface.

The second-order approximation of resultant $\mathsf{\Gamma }\left({\mathbf{x}}_f\right)$ at footpoint ${\mathbf{x}}_f$ can be written as:

$$\begin{array}{cc}\hfill \mathsf{\Gamma }\left({\mathbf{x}}_f\right)& =\mathsf{\Gamma }\left(\mathbf{x}\right)-\frac{\partial \mathsf{\Gamma }\left(\mathbf{x}\right)}{\partial \mathbf{n}}{d}_{quad}+\frac{{\partial }^2\mathsf{\Gamma }\left(\mathbf{x}\right)}{\partial {\mathbf{n}}^2}{d}_{quad}^2\hfill \end{array}$$

(25)

Since the resultant at footpoint is zero, i.e., $\mathsf{\Gamma }\left({\mathbf{x}}_f\right)=0$, Equation (25) can be rearranged as:

$$\begin{array}{cc}\hfill d_{quad}& =\frac{1}{\mathbf{n}·\mathbf{H}·\mathbf{n}}(\nabla \mathsf{\Gamma }·\mathbf{n}\pm \sqrt{{(\nabla \mathsf{\Gamma }·\mathbf{n})}^2-2\mathsf{\Gamma }\mathbf{n}·\mathbf{H}·\mathbf{n}})\hfill \end{array}$$

(26)

where $d_{quad}$ is the distance estimate between the test point and footpoint based on the second-order approximation. $\nabla \mathsf{\Gamma }$ and $\mathbf{H}$ are the gradient and Hessian of algebraic level sets, which are calculated as described in Section 4.1.1. For test points near the curve, a good approximation of normal vector $\mathbf{n}$ could be $\mathbf{n}=\frac{\nabla \mathsf{\Gamma }}{\parallel \nabla \mathsf{\Gamma }\parallel }$ or $\mathbf{n}=\frac{\nabla d}{\parallel \nabla d\parallel }$.

Algebraic point projection in two-dimensional physical space is validated in Figure 9, where the second-order point projection result is compared against first-order point projection as well as Newton–Raphson iterations for various test distances. Both algebraic point projection methods yield the reference solution obtained using the Newton–Raphson method in the limit when the query point is on the curve/surface, and when test distance is small ($\overline{d}=0.1$). The second-order point projection, however, converges to the solution when distances are larger ($\overline{d}=0.3$), where the first-order point projection fails.

4.1.3. Improvement to First Order Algebraic Point Projection

As illustrated in Figure 8 and Figure 9, when the test point is far from the curve, first-order algebraic point projection is inaccurate. In immersed boundary analysis, point projection is only required at quadrature points near the boundary since the physical influence of the boundary on the underlying domain is local. For quadrature points that are far from the boundary, a recursive first-order algebraic point projection must be adopted for accuracy. The main idea is to treat the footpoint estimate of the previous step of first-order algebraic point projection as the starting point for the next step, and to recursively proceed until convergence. If a point is far from the boundary, we combine the distance to the curve d in Equation (14) and distance to convex hull of Bézier segment $\varphi$ to get the composed distance $\tilde{d}$ and its gradient $\nabla \tilde{d}$ as follows:

$$\begin{array}{c}\hfill \tilde{d}=\sqrt{d^2+\frac{{\left(\right|\varphi |-\varphi )}^2}{4}}\end{array}$$

(27)

$$\begin{array}{c}\hfill \nabla \tilde{d}=\left\{\begin{array}{cc}\nabla d,\hfill & \varphi \ge 0\hfill \\ \frac{2\nabla d+2\nabla \varphi }{\sqrt{d^2+{\varphi }^2}},\hfill & \varphi <0\hfill \end{array}\right.\end{array}$$

(28)

The composed distance together with its gradient is then used for estimating the footpoint in Equation (20). The result after recursive execution is illustrated in Figure 10. Comparing against original first-order algebraic point projection in Figure 9b, one can observe that the improved first order procedure is better.

4.1.4. Improvement to Second Order Algebraic Point Projection

While second-order algebraic point projection is more accurate, as should be expected, relative to first-order algebraic point projection, it could be further improved, as discussed below. Since the $\mathbf{n}$ computed locally at a point in the domain is relatively inaccurate for approximating the direction of the vector $\mathbf{x}-{\mathbf{x}}_f$ at greater distances, Equation (26) will fail to provide the correct solution as distance from the curve increases. As illustrated in Figure 11a, at test points along the line where $\mathbf{n}$ does not intersect with the Bézier curve, no valid solution is obtained. This is because the normal vector $\mathbf{n}$ generated using algebraic level sets usually does not coincide with normal vector ${\mathbf{n}}_f$ at footpoint, especially when test point is far away. One solution to this issue is to construct a new vector ${\mathbf{n}}_{corr}$, which is guaranteed to have positive discriminant in Equation (26). A trial ${\mathbf{n}}_{corr}$ that points to the centroid of control polygon achieves the desired result. This choice of the normal direction at points far away from the curve is simple and effective as shown in Figure 11b. While this correction is not general, it is likely to provide sufficient accuracy in practice for quadrature points where the behavior of the immersed boundary is blended with that of the underlying domain.

4.1.5. Inversion to Parametric Space

Given a footpoint ${\mathbf{x}}_f$ in physical space, finding a corresponding parameter $u_f$ such that $\mathbf{C}\left(u_f\right)={\mathbf{x}}_f$ is the point inversion problem. The direct approach to carrying out point inversion is by solving a system of polynomial equations, which may not be easy to do since there is no analytical solution for high-order polynomials of $deg>4$, and since ${\mathbf{x}}_f$ may not lie exactly on curve $\mathbf{C}\left(u\right)$ [22].

This drawback can be overcome by using the Bezout matrix [25], as shown in the following. Evaluate ${\mathbf{M}}^B\left({\mathbf{x}}_f\right)$ with ${\mathbf{x}}_f=(x_f,y_f)$:

$$\begin{array}{c}\hfill {\mathbf{M}}^B\left({\mathbf{x}}_f\right)={\mathbf{M}}_x^B{x}_f+{\mathbf{M}}_y^B{y}_f+{\mathbf{M}}_w^B={\left[m_{ij}\right]}_{p\times p}\end{array}$$

(29)

Then, Equation (6) can be rewritten as the following over-constrained linear system:

$$\begin{array}{c}\hfill A\tilde{\mathbf{u}}\equiv \left[\begin{array}{cccc}{m}_{11}& m_{12}& \cdots & m_{1(p-1)}\\ m_{21}& m_{22}& \cdots & m_{2(p-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ m_{p1}& m_{p2}& \cdots & m_{p(p-1)}\end{array}\right]\left(\begin{array}{c}{u}^{p-1}\\ u^{p-2}\\ ⋮\\ u\end{array}\right)=-\left(\begin{array}{c}{m}_{1p}\\ m_{2p}\\ ⋮\\ m_{pp}\end{array}\right)\end{array}$$

(30)

Matrix A is full ranked if Equation (5) has only one common root, i.e., if ${\mathbf{x}}_f$ is not a double point [30]. Thus, $u_f$ can be obtained by solving a linear least square problem resulting from Equation (30), which requires bounded computational cost unlike numerical iterations using the Newton–Raphson method. In addition, if a physical test point $\mathbf{x}$ is initially on the curve, then ${\mathbf{x}}_f=\mathbf{x}$, and the point inversion can be directly applied, without needing to find the footpoint. Alternative to solving the over-constrained system, one can discard any one row of $\mathbf{A}$ in Equation (30), which will make it full rank and the point inversion solution can be obtained by solving linear system of equations.

To compute $u_f$ on a NURBS curve, which is piece-wise polynomial, projection onto the appropriate Bézier segment is required. For this purpose, one can identify the closest Bézier segment using individual algebraic level sets, and apply algebraic point projection on the closest Bézier segment. Denoting the computed parameter on the closest Bézier curve as $u_f^B$, the corresponding parameter on the original NURBS curve $u_f^N$ can be obtained by purely linear scaling and offset. Unlike iterative or recursive schemes, the algebraic method guarantees the existence of a definite footpoint without needing to manipulate user-controlled parameters such as stop criterion or recursion limit in an effort to coax a solution. If a test point is close to the connection node of two adjacent Bézier segments, a result of $u_f^B<0$ or $u_f^B>1$ may be obtained. In this case, higher projection precision can be achieved when a second point projection to the adjacent Bézier segment is attempted (see Figure 12).

4.2. Extension to NURBS Surfaces

Analogous to the algebraic level sets in three-dimensions, the algebraic point projection can be naturally extended to three-dimensional Bézier and NURBS surfaces by replacing Bezout matrix ${\mathbf{M}}^B$ with Dixon matrix ${\mathbf{M}}^D$. Since the Dixon matrix also has the linearity property, as given in Equation (17), the basic procedure for point projection remains the same, as outlined in Section 4.1.

4.2.1. Projection in Physical Space

Utilizing the linearity property (Equation (17)), one can rewrite $\nabla \mathsf{\Gamma }$ and $\mathbf{H}$ in Equation (19) as follows (as before, superscript D is dropped for ease of reading):

$$\begin{array}{cc}\hfill \nabla \mathsf{\Gamma }& =\left|\mathbf{M}\right|\left[\begin{array}{c}\mathrm{tr}\left({\mathbf{M}}^{-1}\frac{\partial }{\partial x}\mathbf{M}\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}\frac{\partial }{\partial y}\mathbf{M}\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}\frac{\partial }{\partial z}\mathbf{M}\right)\end{array}\right]=\left|\mathbf{M}\right|\left[\begin{array}{c}\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\end{array}\right]=\mathsf{\Gamma }\tilde{\mathbf{g}}\hfill \\ \hfill \mathbf{H}& =\left|\mathbf{M}\right|\left[\begin{array}{ccc}{\mathrm{tr}}^2\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)& \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)& \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\\ -\mathrm{tr}\left({\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)}^2\right)& -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\right)& -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)& {\mathrm{tr}}^2\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)& \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\\ -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\right)& -\mathrm{tr}\left({\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)}^2\right)& -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\right)\\ \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)& \mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\mathrm{tr}\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)& {\mathrm{tr}}^2\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\\ -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_x\right)\right)& -\mathrm{tr}\left(\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)\left({\mathbf{M}}^{-1}{\mathbf{M}}_y\right)\right)& -\mathrm{tr}\left({\left({\mathbf{M}}^{-1}{\mathbf{M}}_z\right)}^2\right)\end{array}\right]\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & =\mathsf{\Gamma }\tilde{\mathbf{H}}\hfill \end{array}$$

(32)

Next, as before, Equations (20), (23) and (24) can be exploited to obtain the physical footpoint ${\mathbf{x}}_f$ in three-dimensional space. Earlier statements on efficiency of the algebraic point projection for rational parametric curves also apply to rational parametric surfaces except that the size of the Dixon matrix ${\mathbf{M}}^D$ is $2pq\times 2pq$, where p and q are the degrees of the rational parametric surface in each dimension. Algebraic point projection in three-dimensional physical space is demonstrated in Figure 13, where test points are projected onto a target Bézier surface using the proposed second-order algebraic point projection method and the Newton–Raphson iterations. Again, one can observe that the proposed method leads to accurate solutions for test points closer to the surface.

4.2.2. Inversion to Parametric Space

The point inversion for rational parametric surfaces can be carried out using Dixon matrix as well. Substituting the physical coordinates of the footpoint ${\mathbf{x}}_f=(x_f,y_f,z_f)$ in ${\mathbf{M}}^D$, we get:

$$\begin{array}{c}\hfill {\mathbf{M}}^D\left({\mathbf{x}}_f\right)={\mathbf{M}}_x^D{x}_f+{\mathbf{M}}_y^D{y}_f+{\mathbf{M}}_z^D{z}_f+{\mathbf{M}}_w^D={\left[m_{ij}\right]}_{2pq\times 2pq}\end{array}$$

(33)

Thus, as before, the homogeneous system in Equation (16) can be converted into an over-constrained non-homogeneous system as follows:

$$\begin{array}{c}\hfill \left[\begin{array}{ccccc}{m}_{12}& \cdots & m_{1(1+i+jp)}& \cdots & m_{1\left(2pq\right)}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ m_{\left(2pq\right)2}& \cdots & m_{\left(2pq\right)(1+i+jp)}& \cdots & m_{\left(2pq\right)\left(2pq\right)}\end{array}\right]\left(\begin{array}{c}{u}^{p-1}{v}^{2q-1}\\ ⋮\\ u^i{v}^j\\ ⋮\\ u\end{array}\right)=-\left(\begin{array}{c}{m}_{11}\\ ⋮\\ m_{\left(2pq\right)1}\end{array}\right)\end{array}$$

(34)

Generally, the parameters $(u_f,v_f)$ of the footpoint ${\mathbf{x}}_f$ can be computed by solving a least square problem or discarding any one row and then solving a linear system of equations as mentioned earlier. As before, the computation of parameters $(u_f^N,v_f^N)$ on a NURBS surface requires sub-division of the surface into a set of Bézier segments. One can apply algebraic point projection to the Bézier segment with smallest algebraic level set measure, and acquire $(u_f^N,v_f^N)$ from the Bézier parameters $(u_f^B,v_f^B)$ by simple linear scaling and offset. As illustrated in Figure 14, a second projection may be necessary when $u_f^B$ or $v_f^B$ is outside the range $[0,1]$.

4.3. Time Complexity of the Algebraic Point Projection Algorithm

The generic pseudo-code of algebraic point projection for both NURBS curves and surfaces is listed in Algorithm 2. As can be seen in Algorithms 1 and 2, algebraic level sets (ALS) and algebraic point projection (APP) are closely connected. Algebraic level sets provide the closest Bézier segment for the first point projection, but also restrict the target curve or surface to be low degree ($p,q\le 3$) so as to avoid double points.

The time complexity for Algorithm 2 is now presented for both curves and surfaces. Consider a NURBS curve of degree p or a NURBS surface of degree $p\times q$, decomposing into n Bézier segments. For each step in finding the point projection, the time complexity is a function of the degree and is given in

Table 3. Time complexity of each step in Algorithm 2 for algebraic point projection. Time complexities are listed for NURBS curves of degree p and NURBS surfaces of degree $p\times q$.

Step	Time Complexity
	Curve	Surface
Computing algebraic level set	$O\left(n{p}^3\right)$	$O\left(n{p}^3{q}^3\right)$
Determining the closest Bézier segment	$O\left(n\right)$	$O\left(n\right)$
Projection in physical space	$O\left(p^3\right)$	$O\left(p^3{q}^3\right)$
Point inversion to parametric space	$O\left(p^3\right)$	$O\left(p^3{q}^3\right)$
Scaling and offset	$O\left(1\right)$	$O\left(1\right)$

. As shown in Section 3.5, computing the algebraic level set at a point costs $O\left(n{p}^3\right)$ time for a curve and $O\left(n{p}^3{q}^3\right)$ time for a surface. Choosing the closest Bézier segment requires finding the minimum of the level sets of n Bézier segments and costs $O\left(n\right)$ time. Projection in the physical space to the closest Bézier segment involves finding the gradient and Hessian of the Bezout/Dixon matrix, which requires the solution to a linear system (see Equations (21) and (22)). This costs $O\left(p^3\right)$ time for curves and $O\left(p^3{q}^3\right)$ time for surfaces. Point inversion to parametric coordinates requires solving a resultant system (Equations (30) and (34)), which also requires $O\left(p^3\right)$ time for curves and $O\left(p^3{q}^3\right)$ time for surfaces. Finally, the parametric coordinates of the Bézier segment are scaled and offset to obtain the parametric coordinates of the NURBS curve/surface. This can be done in $O\left(1\right)$ time. Thus, the total time complexity of algebraic point projection is the same as that of computing the algebraic level sets, i.e., $O\left(n{p}^3\right)$ for NURBS curves and $O\left(n{p}^3{q}^3\right)$ for NURBS surfaces.

Algorithm 2 Algebraic point projection algorithm.

Input: NURBS curve or surface $\mathrm{C}\left(\mathrm{u}\right)$ and given point $\mathrm{x}$ Output: Parameter ${\mathrm{u}}_f^N$ of footpoint on NURBS entity $\mathrm{C}$. 1: functionAlgebraic_Point_Projection($\mathrm{C}$, $\mathrm{x}$) 2:     $\mathcal{B}\left(\mathrm{C}\right)$← Split NURBS entity $\mathrm{C}$ into a Bézier set with segments $B_i,i=1,2,\cdots ,n$ 3:     for $i\leftarrow 1,n$do        ▹ Loops are independent, parallelization is possible 4:         $d_i$←Algebraic_Distance($B_i$, $\mathrm{x}$) 5:     end for 6:     ${\mathrm{x}}_f$←$\mathrm{x}-d_{B_j}\left(\mathrm{x}\right){\mathrm{n}}_{B_j\left(\mathrm{x}\right)}$▹$j={arg\; min}_{i\in [1,n]}{d}_i$, and $d_{B_j}$ is the distance to jth Bezier segemnt using 1st or 2nd order approximation, ${\mathrm{n}}_{B_j\left(\mathrm{x}\right)}$ is normal vector of $d_{B_j}$ 7:     ${\mathrm{u}}_f^B$← Solve Equation (30) with ${\mathrm{M}}^B\left({\mathrm{x}}_f\right)$ or Equation (34) with ${\mathrm{M}}^D\left({\mathrm{x}}_f\right)$ 8:     if ${\mathrm{u}}_f^B$ is out of span then 9:         ${\mathrm{u}}_f^B$← Carryout the second projection based on Figure 12 or Figure 14 10:         if Second projection is infeasible or ${\mathrm{u}}_f^B$ is still out of span then 11:            ${\mathrm{u}}_f^B$← Compute corresponding parameter value on $B_j$ boundary 12:         end if 13:     end if 14:     ${\mathrm{u}}_f^N$← Scale and offset ${\mathrm{u}}_f^B$ based on knot span of $\mathrm{C}$ 15: end function

5. Results and Discussion

Four numerical examples are presented to demonstrate the algebraic point projection methodology of Section 4.1.1 and Section 4.1.2. Two of the examples are projections to curves and the remaining two to surfaces. The curve examples show the performance of algebraic point projection in simple tests to achieve the accuracy of the Newton–Raphson method. The examples also reveal the superiority of second-order point projection over first-order point projection. The surface examples further demonstrate the cheaper computational cost of second-order algebraic point projection when compared against the Newton–Raphson method.

The procedure for algebraic point projection implemented in the examples is summarized as follows and illustrated through a flowchart in Figure 15:

• Decompose NURBS curve or surface and get the corresponding implicit form for each Bézier entity.
• Construct algebraic level sets for the decomposed Bézier entities.
• Compose level sets on Bézier patches using R-functions to construct algebraic level sets on the NURBS entity.
• Find nearest Bézier entity for a given quadrature point in the domain.
• Project test point on the nearest Bézier entity to get foot point in the physical space.
• Apply point inversion and decide if recursive projection is needed.
• Obtain global parametric solution by scaling the parameter value of the Bézier entity to that of the NURBS entity.

5.1. Curve Tests

The first curve example is illustrated in Figure 16, where physical points in the underlying domain are projected onto a given NURBS curve using both Newton–Raphson iterations as well as the algebraic point projection. Contour levels indicate the value of parameter $u_f^N$ of the predicted footpoint.

As can be observed in Figure 16, both first- and second-order algebraic point projection provide exact on-curve solutions and good approximate solutions to the parameter values from points near the curve. In addition, one may observe that there are two regions, at the bottom-left and upper-right of Figure 16a, respectively, where due to high curvature of nearby curve segments, a jump in parameter value occurs when Newton–Raphson iterations are used. Such discontinuities disappear when algebraic point projection is used. The relationship between distance from the curve and relative error is obtained for both first- and second-order algebraic point projection methods using the Newton–Raphson method as a reference (Figure 17). The test points are picked along an arbitrary vertical trace line in Figure 16. The parameter range of the NURBS curve is $[0,1]$ in the example. As the test point moves closer to the target curve, the projection error decreases quadratically (first-order point projection) and cubically (second-order point projection).

The second curve example shown in Figure 18 demonstrates the robustness of the algebraic point projection when the footpoint is either discontinuous or non-existent. One can observe in Figure 18b that, although Points A and B are continuous in terms of the parametric value on the tracing curve, the Newton–Raphson method results in a large jump in parametric solution $A_f^N$ and $B_f^N$, whereas the algebraic method yields a continuous parameter value. In general, the algebraic projected solution is smoother and matches the Newton–Raphson solution well, and second-order point projection is more accurate than first-order point projection.

5.2. Surface Tests

The first and second surface examples demonstrate the robustness of the second-order algebraic point projection for NURBS surfaces involving discontinuous and non-existent footpoints, respectively. In the first example (Figure 19), the discontinuous projection occurs again when the Newton–Raphson method is applied on a surface segment with high mean curvature. In the second example (Figure 20), the Newton–Raphson method failed in regions where the mathematical footpoints do not exist. Not only does the algebraic point projection overcome those problems, but it also produces an accurate and efficient solution. As listed in

Table 4. The results of point projection for NURBS surfaces. The tolerance in Newton–Raphson iterations was chosen as $ϵ={10}^{-8}$. Note that the time of finding an initial point is excluded in the time per iteration.

Example Surface	Newton–Raphson Iterations			2nd order APP
	Time per Point ($\mathsf{\mu }$s)	Average Number of Iterations	Time per Iteration ($\mathsf{\mu }$s)	Time per Point ($\mathsf{\mu }$s)	Time per Point for Algorithm 1 ($\mathsf{\mu }$s)
#1 (Figure 19)	211.91	5.00	36.03	55.28	14.84
#2 (Figure 20)	769.73	10.85	50.47	107.88	16.92

, the computational cost per point of the proposed method is only $26\%$ and $14\%$ of that of the Newton–Raphson method in the first and second surface examples, respectively.

6. Conclusions

In this paper, novel first- and second order-algebraic point projection methods for low degree two-dimensional NURBS curves and three-dimensional NURBS surfaces are proposed. The procedure utilizes the recently developed algebraic level sets. As a first step, the differential property of the resultant matrix is used to obtain the footpoint in the physical space. Next, the parameter value of the footpoint is computed by solving the over-constrained resultant system. Algebraic point projection eliminates inefficient iterative computations and the need for a good initial guess by providing an exact on-curve solution and good approximate solution for points near the curve. Through numerical examples, the algebraic method, especially the second-order point projection is demonstrated to be faster and more smooth than Newton–Raphson iterations.

The proposed algebraic point projection technique possesses two limitations. As illustrated in Section 4, the proposed method is inaccurate or will fail when query points are far away from the curve/surface due to the inaccuracy of the algebraic distance measure. Algebraic point projection will also fail at points where the gradient of the algebraic level-sets is not defined. Such a point occurs, for example, at the center of a circle or sphere. At these locations, a correction to the gradient or a perturbed point is needed. Although algebraic point projection is an approximation, particularly for test points far away from the curve/surface, it has utility for isogeometric analysis, where small inaccuracy in point inversion solution may be acceptable as a trade-off against smooth, robust, and efficient projection. Such a projection is critical given the large number of quadrature points at which point projection is necessary during analysis.