Isogeometric Analysis: Knot Insertion for NURBS

Introduction

In the classical Finite Element Method (FEM) the domain of the problem is approximated by defining a mesh and basis functions over its elements. More elements for the mesh means more basis functions and therefore larger space where approximated solutions can be found. The procedure of increasing the number of elements in the mesh is known as h-refinement. In Isogeometric Analysis (IGA), no mesh is required as the basis functions used to define the domain are the same basis functions used to define the approximate solution. Nevertheless, in IGA, the solution is still a combination of basis functions, so increasing the number of basis functions is still needed to implement some kind of h-refinement. A good way to do this is through knot insertion into the NURBS. This post also presents demo written in Octave and TypeScript.

You can find the full source code for the algorithms and all the examples in this GitHub repo.

Problem of Knot Insertion

Given the general form of a NURBS curve as described here:

$$\boldsymbol{C}^{w}\left(\xi\right)=\sum_{i=0}^{n}N^{p}_{i}\left(\xi\right)\boldsymbol{P}_i^w$$

on the knot vector:

$$\Xi=\left[\xi_0,…,\xi_m\right]$$

knot insertion is the problem of adding the knot $\bar{\xi}\in\left[\xi_k,\xi_{k+1}\right)$ in the knot vector $\Xi$, so that the following still stands:

$$\bar{\boldsymbol{C}}^w\left(\xi\right)=\sum_{i=0}^{n+1}N_{i}^{p}\left(\xi\right)\bar{\boldsymbol{P}}^{w}_{i}=\boldsymbol{C}^w\left(\xi\right),\forall\xi\in\left[\xi_0,\xi_{m+1}\right]$$

on the knot vector:

$$\Xi=\left[\xi_0,…,\xi_{m+1}\right]$$

where $\bar{\boldsymbol{P}}_i^w$ are the control points of the new curve in the projective space.

Knot insertion is supposed to keep the curve unaltered by adding and modifying control points and knots.

Solutions

The problem can then be solved by finding the $n+1$ values $\bar{\boldsymbol{P}}^w_i$, which can be done solving the linear system in $n+2$ equations:

$$\sum_{i=0}^{n+1}N_{i}^{p}\left(\xi\right)\bar{\boldsymbol{P}}^{w}_{i}=\sum_{i=0}^{n}N^{p}_{i}\left(\xi\right)\boldsymbol{P}^w_{i},\forall\xi\in\left\{\xi_0,…,\xi_{n+1}\right\}$$

There is also a simpler way to insert the knots. It can be shown that:

$$\displaylines{\bar{\boldsymbol{P}}_{i}^{w}=\alpha_{i}\boldsymbol{P}_{i}^{w}+\left(1-\alpha_{i}\right)\boldsymbol{P}_{i-1}^{w},\\
\alpha_{i}=\left\{ \begin{array}{ll} 1, & i\leq k-p\\ \dfrac{\bar{\xi}-\xi_{i}}{\xi_{i+p}-\xi_{i}}, & k-p+1\leq i\leq k\\ 0, & i\geq k+1 \end{array}\right.}$$

This is the technique used in all the implementations in the repo.

Octave Implementation

The Octave implementation includes the computeKnotInsertion script that takes a NURBS as input, together with the knot value to add and returns a new knot vector and control points. This is an example of insertion into a circle:

The computeSurfKnotInsertion script computes knot insertion for NURBS surfaces instead. With the script, knots can be added to both dimensions. An example is shown below:

TypeScript Implementation

The TypeScript implementation includes a method named insertKnot in the NurbsCurve class. That method allows to insert multiple knots into the vector, returning a new NurbsCurve, identical to the first, but with a different knot vector and different control points.

With a few lines of code:

let circle = new NurbsCirle()
drawNurbsCurve(circle.controlPoints, circle.knotVector, circle.weights, 2, false, true, plot1, null, true, "Knot Insertion 1")
circle.insertKnot(0.6, 6, 0, 1)
drawNurbsCurve(circle.controlPoints, circle.knotVector, circle.weights, 2, false, true, plot2, null, true, "Knot Insertion 2")
circle.insertKnot(0.3, 4, 0, 1)
drawNurbsCurve(circle.controlPoints, circle.knotVector, circle.weights, 2, false, true, plot3, null, true, "Knot Insertion 3")
circle.insertKnot(0.2, 2, 0, 1)
drawNurbsCurve(circle.controlPoints, circle.knotVector, circle.weights, 2, false, true, plot4, null, true, "Knot Insertion 4")

you can draw this result:

The NurbsSurf class includes the methods insertKnotXi and insertKnotEta to insert knots in the two dimensions of the parametric space and returns the new object NurbsSurf. This is the result of the algorithm applied to a plate:

and this is the result on a toroid (in this case knots were not inserted in the entire vector uniformly, but only in a section):

As expected, after each insertion, the NURBS is unaltered. This can also be seen by running the unit tests provided in the repo. For example, one unit test includes this code:

function testNurbs(n1: NurbsSurf, n2: NurbsSurf) {
    for (let xi = 0; xi <= 1; xi += 0.01)
        for (let eta = 0; eta <= 1; eta += 0.01)
            assert(approxEqualPoints(n1.evaluate(xi, eta), n2.evaluate(xi, eta)))
}

[...]

measure("surf_knot_insertion_xi", () => {
    let n1 = new NurbsPlateHole()
    let n2 = new NurbsPlateHole()
    for (let i = 0.1; i <= 1; i += 0.9) {
        if (n2.Xi[i] != i) {
            let k = BsplineCurve.findSpan(n2.Xi, i, n2.p, n2.controlPoints.length - 1)
            testNurbs(n1, n2.insertKnotsXi(i, k, 0, 1))
        }
    }
})

Bye! 😉