Mathematical utilities

B-Splines and NURBS

We start this section by recalling some basic properies about B-splines curves and surfaces. We also recall some fundamental algorithms (knot insertion and degree elevation).

For a basic introduction to the subject, we refer to the books [LP95] and [Far02].

A B-Splines family, (N_i)_{ 1 \leqslant i \leqslant n} of order k, can be generated using a non-decreasing sequence of knots T=(t_i)_{1\leqslant i \leqslant n + k}.

B-Splines series

The j-th B-Spline of order k is defined by the recurrence relation:

N_j^k = w_j^k N_j^{k-1} + ( 1 - w_{j+1}^k ) N_{j+1}^{k-1}

where,

w_j^k (x) = \frac{x-t_j}{t_{j+k-1}-t_{j}} \hspace{2cm} N_j^1(x) = \chi_{ \left[ t_j, t_{j+1} \right[ }(x)

for k \geq 1 and 1 \leq j \leq n.

We note some important properties of a B-splines basis:

  • B-splines are piecewise polynomial of degree p=k-1,
  • Compact support; the support of N_j^k is contained in \left[ t_j, t_{j+k} \right] ,
  • If x \in~ ] t_j,t_{j+1} [, then only the B-splines \{ N_{j-k+1}^k,\cdots,N_{j}^k \} are non vanishing at x,
  • Positivity: \forall j \in \{1,\cdots,n \}~~N_j(x) >0, ~~\forall x \in ] t_j, t_{j+k} [,
  • Partition of unity \sum_{i=1}^n N_i^{k}(x) = 1, \forall x \in \mathbb{R},
  • Local linear independence,
  • If a knot t_i has a multiplicity m_i then the B-spline is \mathcal{C}^{(p-m_i)} at t_i.

Knots vector families

There are two kind of knots vectors, called clamped and unclamped. Both families contains uniform and non-uniform sequences.

The following are examples of such knots vectors

  1. Clamped knots (open knots vector)
  • uniform

T_1 &= \{0, 0, 0, 1, 2, 3, 4, 5, 5, 5 \}
\\
T_2 &= \{-0.2, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 0.8 \}

_images/bsplines_t1_p2.png _images/bsplines_t2_p2.png
  • non-uniform

T_3 &= \{0, 0, 0, 1, 3, 4, 5, 5, 5 \}
\\
T_4 &= \{-0.2, -0.2, 0.4, 0.6, 0.8, 0.8 \}

_images/bsplines_t3_p2.png _images/bsplines_t4_p2.png
  1. Unclamped knots
  • uniform

T_5 &= \{0, 1, 2, 3, 4, 5, 6, 7 \}
\\
T_6 &= \{-0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 \}

_images/bsplines_t5_p2.png _images/bsplines_t6_p2.png
  • non-uniform

T_7 &= \{0, 0, 3, 4, 7, 8, 9 \}
\\
T_8 &= \{-0.2, 0.2, 0.4, 0.6, 1.0, 2.0, 2.5 \}

_images/bsplines_t7_p2.png _images/bsplines_t8_p2.png

B-Spline curve

The B-spline curve in \mathbb{R}^d associated to knots vector T=(t_i)_{1\leqslant i \leqslant n + k} and the control polygon (\mathbf{P}_i)_{ 1 \leqslant i \leqslant n} is defined by :

\mathcal{C}(t) = \sum_{i=1}^n N_i^k(t) \textbf{P}_i

In (Fig. ref{figBSplineCurve}), we give an example of a quadratic B-Spline curve, and its corresponding knot vector and control points.

_images/courbe_bsplines.png _images/basis_fct_p2_N5.png

We have the following properties for a B-spline curve:

  • If n=k, then \mathcal{C} is just a B’ezier-curve,
  • \mathcal{C} is a piecewise polynomial curve,
  • The curve interpolates its extremas if the associated multiplicity of the first and the last knot are maximum (i.e. equal to k), i.e. open knot vector,
  • Invariance with respect to affine transformations,
  • Strong convex-hull property:

if t_i \leq t \leq t_{i+1}, then \mathcal{C}(t) is inside the convex-hull associated to the control points \mathbf{P}_{i-p},\cdots,\mathbf{P}_{i},

  • Local modification : moving the i^{th} control point \mathbf{P}_{i} affects \mathcal{C}(t), only in the interval [t_i,t_{i+k}],
  • The control polygon approaches the behavior of the curve.

Note

In order to model a singular curve, we can use multiple control points : \mathbf{P}_{i}=\mathbf{P}_{i+1}.

Multivariate tensor product splines

Let us consider d knot vectors \mathcal{T} = \{T^1,T^2,\cdots,T^d\}. For simplicity, we consider that these knot vectors are open, which means that k knots on each side are duplicated so that the spline is interpolating on the boundary, and of bounds 0 and 1. In the sequel we will use the notation I=[0,1]. Each knot vector T^i, will generate a basis for a Schoenberg space, \mathcal{S}_{k_{i}}(T^i,I). The tensor product of all these spaces is also a Schoenberg space, namely \mathcal{S}_{\mathbf{k}}(\mathcal{T}), where \mathbf{k}=\{k_1,\cdots,k_d\}. The cube \mathcal{P}=I^d=[0,1]^d, will be referred to as a patch.

The basis for \mathcal{S}_{\mathbf{k}}(\mathcal{T}) is defined by a tensor product :

N_{\mathbf{i}}^{\mathbf{k}} := N_{i_1}^{k_1} \otimes N_{i_2}^{k_2} \otimes \cdots \otimes N_{i_d}^{k_d}

where, \mathbf{i}=\{i_1,\cdots , i_d \}.

A typical cell from \mathcal{P} is a cube of the form : Q_{\mathbf{i}}=[\xi_{i_1}, \xi_{i_1+1}] \otimes \cdots \otimes [\xi_{i_d}, \xi_{i_d+1}].

Deriving a B-spline curve

The derivative of a B-spline curve is obtained as:

\mathcal{C}^{\prime}(t) = \sum_{i=1}^{n} {N_{i}^{k}}^{\prime}(t) \mathbf{P}_i = \sum_{i=1}^{n} \left(\frac{p}{t_{i+p}-t_{i}}N_{i}^{k-1}(t) \mathbf{P}_i - \frac{p}{t_{i+1+p}-t_{i+1}}N_{i+1}^{k-1}(t) \mathbf{P}_i \right)
= \sum_{i=1}^{n-1} {N_{i}^{k-1}}^{\ast}(t) \mathbf{Q}_i

where \mathbf{Q}_i = p \frac{\mathbf{P}_{i+1} - \mathbf{P}_i}{t_{i+1+p}-t_{i+1}}, and \{{N_{i}^{k-1}}^{\ast},~~1 \leq i \leq n-1\} are generated using the knot vector T^{\ast}, which is obtained from T by reducing by one the multiplicity of the first and the last knot (in the case of open knot vector), i.e. by removing the first and the last knot.

More generally, by introducing the B-splines family \{ {N_{i}^{k-j}}^{\ast}, 1 \leq i \leq n-j \} generated by the knots vector T^{j^{\ast}} obtained from T by removing the first and the last knot j times, we have the following result:

proposition

The j^{th} derivative of the curve \mathcal{C} is given by

\mathcal{C}^{(j)}(t) = \sum_{i=1}^{n-j} {N_{i}^{k-j}}^{\ast}(t) \mathbf{P}_i^{(j)}`

where, for j>0

\mathbf{P}_i^{(j)} = \frac{p-j+1}{t_{i+p+1}-t_{i+j}} \left( \mathbf{P}_{i+1}^{(j-1)} - \mathbf{P}_i^{(j-1)} \right)
\\
\mbox{and} ~ ~ ~ \mathbf{P}_i^{(0)} = \mathbf{P}_i.

By denoting \mathcal{C}^{\prime} and \mathcal{C}^{\prime\prime} the first and second derivative of the B-spline curve \mathcal{C}, it is easy to show that:

We have,

  • \mathcal{C}^{\prime}(0) = \frac{p}{t_{p+2}} \left(\mathbf{P}_{2} - \mathbf{P}_1\right),
  • \mathcal{C}^{\prime}(1) = \frac{p}{1-t_{n}} \left(\mathbf{P}_{n} - \mathbf{P}_{n-1}\right),
  • \mathcal{C}^{\prime\prime}(0) = \frac{p(p-1)}{t_{p+2}} \left( \frac{1}{t_{p+2}}\mathbf{P}_{1} - \{ \frac{1}{t_{p+2}} + \frac{1}{t_{p+3}} \} \mathbf{P}_2 + \frac{1}{t_{p+3}}\mathbf{P}_{3} \right),
  • \mathcal{C}^{\prime\prime}(1) = \frac{p(p-1)}{1-t_{n}} \left( \frac{1}{1-t_{n}}\mathbf{P}_{n} - \{ \frac{1}{1-t_{n}} + \frac{1}{1-t_{n-1}} \} \mathbf{P}_{n-1} + \frac{1}{1-t_{n-1}}\mathbf{P}_{n-2} \right).

Example

Let us consider the quadratic B-spline curve associated to the knots vector T=\{000~\frac{2}{5}~\frac{3}{5}~111 \} and the control points \{ P_i, 1 \leq i \leq 5 \}:

\mathcal{C}(t) = \sum_{i=1}^{5} {N_{i}^{3}}^{\prime}(t) \mathbf{P}_i

we have,

\mathcal{C}^{\prime}(t) = \sum_{i=1}^{4} {N_{i}^{2}}^{\ast}(t) \mathbf{Q}_i

where

\mathbf{Q}_1 = 5 \{\mathbf{P}_{2} - \mathbf{P}_1\}, ~~~~\mathbf{Q}_2 = \frac{10}{3} \{ \mathbf{P}_{3} - \mathbf{P}_2\},
\\
\mathbf{Q}_3 = \frac{10}{3} \{ \mathbf{P}_{4} - \mathbf{P}_3\},~~~~\mathbf{Q}_4 = 5 \{\mathbf{P}_{5} - \mathbf{P}_4\}.

The B-splines \{ {N_{i}^{2}}^{\ast},~~1 \leq i \leq 4\} are associated to the knot vector T^{\ast}=\{00~\frac{2}{5}~\frac{3}{5}~11 \}.

Fundamental geometric operations

By inserting new knots into the knot vector, we add new control points without changing the shape of the B-Spline curve. This can be done using the DeBoor algorithm [dB01]. We can also elevate the degree of the B-Spline family and keep unchanged the curve [HHM05]. In (Fig. ref{refinement_curve_B_Spline}), we apply these algorithms on a quadratic B-Spline curve and we show the position of the new control points.

Knot insertion

After modification, we denote by \widetilde{n}, \widetilde{k}, \widetilde{T} the new parameters. (\textbf{Q}_i) are the new control points.

One can insert a new knot t, where t_j \leqslant t < t_{j+1}. For this purpose we use the DeBoor algorithm [dB01]:

\widetilde{n} = n+1
\\
\widetilde{k} = k
\\
\widetilde{T} = \{ t_1,.., t_j, t, t_{j+1},.., t_{n+k}\}
\\
\alpha_i = \left\{\begin{array}{cc}1 & 1 \leqslant i \leqslant j-k+1 \\\frac{t-t_i}{t_{i+k-1}-t_i} & j-k+2 \leqslant i \leqslant j \\0 & j+1 \leqslant i \end{array}\right.
\\
\textbf{Q}_i = \alpha_i \textbf{P}_i + (1-\alpha_i) \textbf{P}_{i-1}

Many other algorithms exist, like blossoming for fast insertion algorithm. For more details about this topic, we refer to [NT93].

Order elevation

We can elevate the order of the basis, without changing the curve. Several algorithms exist for this purpose. We used the one by Huang et al. [PP91], [HHM05].

A quadratic B-spline curve and its control points. The knot vector is T = \{ 000, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, 1 1 1 \}.

_images/curve.png

The curve after a h-refinement by inserting the knots \{ 0.15, 0.35\} while the degree is kept equal to 2.

_images/curve_p0_n9.png

The curve after a p-refinement, the degree was raised by 1 (using cubic B-splines).

_images/curve_p2_n0.png

The curve after duplicating the multiplicity of the internal knots \{ \frac{1}{4}, \frac{1}{2}, \frac{3}{4} \}, this leads to a B’ezier description. We can then, split the curve into 4 pieces (sub-domains), each one will corresponds to a quadratic B’ezier curve.

_images/curve_p0_n3_bezier.png

Translation

Rotation

Todo

not yet available

Scaling

Todo

not yet available

References

[dB01](1, 2) C. de Boor. A Practical Guide to Splines. Applied Mathematical Sciences. Springer New York, 2001. ISBN 9780387953663. URL: https://books.google.de/books?id=m0QDJvBI_ecC.
[Far02]G. Farin. Curves and surfaces for CAGD: a practical guide. Morgan Kaufmann Pub. Inc., San Francisco, CA, USA, 2002. ISBN 1-55860-737-4.
[HHM05](1, 2) Qi-Xing Huang, Shi-Min Hu, and Ralph R. Martin. Fast degree elevation and knot insertion for b-spline curves. Computer Aided Geometric Design, 22(2):183 – 197, 2005. URL: http://www.sciencedirect.com/science/article/B6TYN-4DXBTHR-2/2/d5b3eec2f4c230c8051623c1c000beae, doi:DOI: 10.1016/j.cagd.2004.11.001.
[LP95]W. Tiller L. Piegl. The NURBS Book. Springer-Verlag, Berlin, Heidelberg, 1995. second ed.
[NT93]Goldman R. N. and Lyche T. Knot Insertion and Deletion Algorithms for B-Spline Curves and Surfaces. SIAM, Philadelphia, USA, 1993. ISBN 9780898713060.
[PP91]Hartmut Prautzsch and Bruce Piper. A fast algorithm to raise the degree of spline curves. Comput. Aided Geom. Des., 8:253–265, October 1991. URL: http://portal.acm.org/citation.cfm?id=124930.124932, doi:10.1016/0167-8396(91)90015-4.

DeRham sequences

Section author: A. Ratnani

here without boundary conditions

\mathbb{R} \hookrightarrow \Hgrad  \xrightarrow{\quad \Grad \quad}  \Hcurl  \xrightarrow{\quad \Curl \quad}   \Hdiv  \xrightarrow{\quad \Div \quad}  \Ltwo  \xrightarrow{\quad} 0

Pullbacks

In the case where the physical domain \Omega := \mathcal{F}(\hat{\Omega}) is the image of a logical domain \hat{\Omega} by a smooth mapping \mathcal{F} (at least \mathcal{C}^1), we have the following parallel diagrams

\begin{array}{ccccccc}
\Hgrad & \xrightarrow{\quad \Grad \quad} & \Hcurl & \xrightarrow{\quad \Curl \quad} &  \Hdiv & \xrightarrow{\quad \Div \quad} & \Ltwo \\
\igrad \Bigg\uparrow   &     & \icurl \Bigg\uparrow  &   & \idiv \Bigg\uparrow &  & \iltwo \Bigg\uparrow       \\
\HgradLogical & \xrightarrow{\quad \Grad \quad} & \HcurlLogical & \xrightarrow{\quad \Curl \quad} &  \HdivLogical & \xrightarrow{\quad \Div \quad} & \LtwoLogical \\
%
\end{array}

Where the mappings \igrad, \icurl, \idiv and \iltwo are called pullbacks and are given by

\phi (x) :=& \igrad \hat{\phi} (\hat{x}) = \hat{\phi}(\mathcal{F}^{-1}(x))
\\
\Psi (x) :=& \icurl \hat{\Psi} (\hat{x}) = \left( D \mathcal{F} \right)^{-T} \hat{\Psi}(\mathcal{F}^{-1}(x))
\\
\Phi (x) :=& \idiv \hat{\Phi} (\hat{x})  = \frac{1}{J} D \mathcal{F} \hat{\Phi}(\mathcal{F}^{-1}(x))
\\
\rho (x) :=& \iltwo \hat{\rho} (\hat{x}) = \hat{\rho}(\mathcal{F}^{-1}(x))

where D \mathcal{F} is the jacobian matrix of the mapping \mathcal{F}.

Note

The pullbacks \igrad, \icurl, \idiv and \iltwo are isomorphisms between the corresponding spaces.

Discrete Spaces

Let us suppose that we have a sequence of finite subspaces for each of the spaces involved in the DeRham sequence. The discrete DeRham sequence stands for the following commutative diagram between continuous and discrete spaces

\begin{array}{ccccccc}
\Hgrad & \xrightarrow{\quad \Grad \quad} & \Hcurl & \xrightarrow{\quad \Curl \quad} &  \Hdiv & \xrightarrow{\quad \Div \quad} & \Ltwo \\
\Pigrad \Bigg\downarrow   &     & \Picurl \Bigg\downarrow  &   & \Pidiv \Bigg\downarrow &  & \Piltwo \Bigg\downarrow       \\
\Vgrad & \xrightarrow{\quad \Grad \quad} & \Vcurl & \xrightarrow{\quad \Curl \quad}  & \Vdiv & \xrightarrow{\quad \Div \quad} & \Vltwo    \\
%
\end{array}

When using a Finite Elements methods, we often deal with a reference element, and thus we need also to apply the pullbacks on the discrete spaces. In fact, we have again the following parallel diagram

\begin{array}{ccccccc}
\Vgrad & \xrightarrow{\quad \Grad \quad} & \Vcurl & \xrightarrow{\quad \Curl \quad} &  \Vdiv & \xrightarrow{\quad \Div \quad} & \Vltwo \\
\igrad \Bigg\uparrow   &     & \icurl \Bigg\uparrow  &   & \idiv \Bigg\uparrow &  & \iltwo \Bigg\uparrow       \\
\VgradLogical & \xrightarrow{\quad \Grad \quad} & \VcurlLogical & \xrightarrow{\quad \Curl \quad} &  \VdivLogical & \xrightarrow{\quad \Div \quad} & \VltwoLogical \\
%
\end{array}

Since, the pullbacks are isomorphisms in the previous diagram, we can define a one-to-one correspondance

\phi :=& \igrad \hat{\phi}, \quad \phi \in \Vgrad, \hat{\phi} \in \VgradLogical
\\
\Psi :=& \icurl \hat{\Psi}, \quad \Psi \in \Vcurl, \hat{\Psi} \in \VcurlLogical
\\
\Phi :=& \idiv \hat{\Phi}, \quad \Phi \in \Vdiv, \hat{\Phi} \in \VdivLogical
\\
\rho :=& \iltwo \hat{\rho}, \quad \rho \in \Vltwo, \hat{\rho} \in \VltwoLogical

We have then, the following results

\Grad \phi =& \icurl \Grad \hat{\phi} , \quad \phi \in \Vgrad
\\
\Curl \Psi =& \idiv \Curl \hat{\Psi} , \quad \Psi \in \Vcurl
\\
\Div \Phi =& \iltwo \Div \hat{\Phi} , \quad \Phi \in \Vdiv

Projectors

In some cases, one may need to define projectors on smooth functions

\begin{array}{ccccccc}
%\Hgradzero & \xrightarrow{\quad \Grad \quad} & \Hcurlzero & \xrightarrow{\quad \Curl \quad}  & \Hdivzero & \xrightarrow{\quad \Div \quad} & \Ltwozero    \\
\Cinfinity & \xrightarrow{\quad \Grad \quad} & \Cinfinity & \xrightarrow{\quad \Curl \quad}  & \Cinfinity & \xrightarrow{\quad \Div \quad} & \Cinfinity    \\
\Pigrad \Bigg\downarrow   &     & \Picurl \Bigg\downarrow  &   & \Pidiv \Bigg\downarrow &  & \Piltwo \Bigg\downarrow       \\
\Vgrad & \xrightarrow{\quad \Grad \quad} & \Vcurl & \xrightarrow{\quad \Curl \quad}  & \Vdiv & \xrightarrow{\quad \Div \quad} & \Vltwo    \\
%\Vgradspline &  \xrightarrow{\quad \Grad \quad}&   \Vcurlspline &  \xrightarrow{\quad \Curl \quad}  &  \Vdivspline&   \xrightarrow{\quad \Div \quad} &  \Vltwospline \\
%\Vgrad &\xrightarrow{\Grad}  & \Vcurl &\xrightarrow{\Curl}   & \Vdiv & \xrightarrow{\Div} & \Vltwo   \\
%
\end{array}

Discrete DeRham sequence for B-Splines

Buffa et al [BSV09] show the construction of a discrete DeRham sequence using B-Splines.

  1. DeRham sequence is reduced to

\mathbb{R} \hookrightarrow
\underbrace{\mathcal{S}^{p}}_{\VgradLogical}  \xrightarrow{\quad \Grad \quad}
\underbrace{\mathcal{S}^{p-1}}_{\VltwoLogical}  \xrightarrow{\quad} 0

  1. The recursion formula for derivative writes

{N_i^p}'(t)=D_i^{p}(t)-D_{i+1}^{p}(t)
\quad \mbox{where} \quad
D_{i}^{p}(t) = \frac{p}{t_{i+p+1}-t_i}N_i^{p-1}(t)

  1. we have \mathcal{S}^{p-1} = \mathbf{span}\{ N_i^{p-1}, 1 \leq i \leq n-1 \} = \mathbf{span}\{ D_i^p, 1 \leq i \leq n-1 \} which is a change of basis as a diagonal matrix
  2. Now if $u in S^p$, with and expansion $u = sum_i u_i N_i^p$, we have

\begin{align*}
  u^{\prime} = \sum_i u_i \left( N_i^p \right)^{\prime} = \sum_i (-u_{i-1} + u_i) D_i^p
%  \label{}
  \end{align*}

  1. If we introduce the B-Splines coefficients vector \mathbf{u} := \left( u_i \right)_{1 \leq i \leq n} $ (and $\mathbf{u}^{\star} for the derivative), we have

\mathbf{u}^{\star} = D \mathbf{u}

where D is the incidence matrix (of entries -1 and +1)

Let I be the identity matrix, we have

in the 2D case:

\mathcal{G} =
\begin{pmatrix}
  D \otimes I
  \\
  I \otimes D
\end{pmatrix}

\mathcal{C} =
\begin{pmatrix}
  I \otimes D
  \\
- D \otimes I
\end{pmatrix}
\quad \mbox{[scalar curl],} \quad
\mathcal{C} =
\begin{pmatrix}
- I \otimes D
 &
  D \otimes I
\end{pmatrix}
\quad \mbox{[vectorial curl]}

\mathcal{D} =
\begin{pmatrix}
  D \otimes I
 &
  I \otimes D
\end{pmatrix}

in the 3D case:

\mathcal{G} =
\begin{pmatrix}
  D \otimes I \otimes I
  \\
  I \otimes D \otimes I
  \\
  I \otimes I \otimes D
\end{pmatrix}

\mathcal{C} =
\begin{pmatrix}
  0    &    - I \otimes I \otimes D     &     I \otimes D \otimes I
  \\
  I \otimes I \otimes D   &    0   &   - D \otimes I \otimes I
  \\
  - I \otimes D \otimes I  & D \otimes I \otimes I & 0
\end{pmatrix}

\mathcal{D} =
\begin{pmatrix}
  D \otimes I \otimes I
  &
  I \otimes D \otimes I
  &
  I \otimes I \otimes D
\end{pmatrix}

References

[BSV09]A. Buffa, G. Sangalli, and R. Vazquez. Isogeometric analysis in electromagnetics: b-splines approximation. Comput. Methods Appl. Mech. Engrg, 199:1143–1152, 2009.

GLT

Where do the GLTs come from?

The main aim of this paragraph is to present a crucial example that highlights the importance of the GLT algebra when dealing with linear systems coming from the discretization of PDEs. Let us start with some preliminaries. In detail, we will recall the notion of symbol of a matrix-sequence and the basic idea behind the GLT theory.

Spectral preliminaries

The following one is a rather informal definition of symbol of a matrix-sequence.

example:

When d_n=n, d=1, D=[0,\pi], \{A_n\}_n\sim_\lambda f means

References

[GMP+14]Carlo Garoni, Carla Manni, Francesca Pelosi, Stefano Serra-Capizzano, and Hendrik Speleers. On the spectrum of stiffness matrices arising from isogeometric analysis. Numerische Mathematik, 127(4):751–799, 2014. URL: http://dx.doi.org/10.1007/s00211-013-0600-2, doi:10.1007/s00211-013-0600-2.