Isogeometric analysis

Isogeometric analysis is a computational approach that offers the possibility of integrating finite element analysis (FEA) into conventional NURBS-based CAD design tools. Currently, it is necessary to convert data between CAD and FEA packages to analyse new designs during development, a difficult task since the two computational geometric approaches are different. Isogeometric analysis employs complex NURBS geometry (the basis of most CAD packages) in the FEA application directly. This allows models to be designed, tested and adjusted in one go, using a common data set.^[1]

The pioneers of this technique are Tom Hughes and his group at The University of Texas at Austin. A reference free software implementation of some isogeometric analysis methods is GeoPDEs.^[2][3] Likewise, other implementations can be found online. For instance, PetIGA^[4] is an open framework for high performance isogeometric analysis heavily based on PETSc. In addition, MIGFEM is another IGA code which is implemented in Matlab and supports Partition of Unity enrichment IGA for 2D and 3D fracture. Furthermore, G+Smo^[5] is an open C++ library for isogeometric analysis. In particular, FEAP^[6] is a finite element analysis program which includes an Isogeometric analysis library FEAP IsoGeometric (Version FEAP84 & Version FEAP85).

Isogeometric analysis presents two main advantages with respect to the finite element method:^[1][7]

• There is no geometric approximation error, due to the fact the domain is represented exactly^[1]
• Wave propagation problems, arising for instance in cardiac electrophysiology, acoustics and elastodynamics, are better described, thanks to the reduction of numerical dispersion and dissipation errors.^[7]

In the framework of IGA, the notions of both control mesh and physical mesh are defined.^[1]

A control mesh is made by the so-called control points and it is obtained by a piecewise linear interpolation of them. Control points play also the role of degrees of freedom (DOFs).^[1]

The physical mesh lays directly on the geometry and it consists of patches and knot spans. According to the number of patches that are used in a specific physical mesh, a single-patch or a multi-patch approach is effectively employed. A patch is mapped from a reference rectangle in two dimensions and from a reference cuboid in three dimensions: it can be seen as the entire computational domain or a smaller portion of it. Each patch can be decomposed into knot spans, which are points, lines and surfaces in 1D, 2D and 3D, respectively. Knots are inserted inside knot spans and define the elements. Basis functions are ${\displaystyle C^{p-m}}$ across the knots, with ${\displaystyle p}$ degree of the polynomial and ${\displaystyle m}$ multiplicity of a specific knot, and ${\displaystyle C^{\infty }}$ between a certain knot and the next or preceding one.^[1]

A knot vector, normally indicated as ${\displaystyle \Xi =\{\xi _{1},\xi _{2},...,\xi _{n+p+1}\}}$, is a set of non-descending points. ${\displaystyle \xi _{i}\in \mathbb {R} }$ is the ${\displaystyle i^{th}}$ knot, ${\displaystyle n}$ is the number of functions, ${\displaystyle p}$ refers to the basis functions order. A knot divides the knot span into elements. A knot vector is uniform or non-uniform according to the fact that its knots, once their multiplicity is not taken into account, are equidistant or not. If the first and the last knots appear ${\displaystyle p+1}$ times, the knot vector is said to be open.^[1][7]

Once a definition of knot vector is provided, several types of basis functions can be introduced in this context, such as B-splines, NURBS and T-splines.^[1]

B-splines can be derived recursively from a piecewise constant function with ${\displaystyle p=0}$:^[1]

${\displaystyle N_{i,0}(\xi )={\mathcal {I}}_{[\xi _{i},\xi _{i+1})}(s)\quad 1\leq i\leq n}$

Using De Boor's algorithm, it is possible to generate B-splines of arbitrary order ${\displaystyle p}$:^[1]

${\displaystyle N_{i,p}(\xi )={\frac {\xi -\xi _{i}}{\xi _{i+p}-\xi _{i}}}N_{i,p-1}(\xi )+{\frac {\xi _{i+p+1}-\xi }{\xi _{i+p+1}-\xi _{i+1}}}N_{i+1,p-1}(\xi )\quad 1\leq i\leq n}$

valid for both uniform and non-uniform knot vectors. For the previous formula to work properly, let the division of two zeros to be equal to zero, i.e. ${\displaystyle {\dfrac {0}{0}}:=0}$.

B-splines that are generated in this way own both the partition of unity and positivity properties, i.e.:^[1]

${\displaystyle \sum _{i=1}^{n}N_{i,p}(\xi )=1\quad \forall \xi }$

${\displaystyle N_{i,p}(\xi )\geq 0\quad \forall \xi }$

So as to calculate derivatives or order ${\displaystyle k}$ of the ${\displaystyle i^{th}}$ B-splines of degree ${\displaystyle p}$, another recursive formula can be employed:^[1]

${\displaystyle {\frac {d^{k}}{d^{k}\xi }}N_{i,p}(\xi )={\frac {p!}{(p-k)!}}\sum _{j=0}^{k}\alpha _{k,j}N_{i+j,p-k}(\xi )}$

where:

${\displaystyle \alpha _{0,0}=1}$

${\displaystyle \alpha _{k,0}={\frac {\alpha _{k-1,0}}{\xi _{i+p-k+1}-\xi _{i}}},}$

${\displaystyle \alpha _{k,j}={\frac {\alpha _{k-1,j}-\alpha _{k-1,j-1}}{\xi _{i+p+j-k+1}-\xi _{i+j}}}\quad j=1,...,k-1,}$

${\displaystyle \alpha _{k,k}={\frac {-\alpha _{k-1,j-1}}{\xi _{i+p+1}-\xi _{i+k}}}.}$

whenever the denominator of an ${\displaystyle \alpha _{j,j}}$ coefficient is zero, the entire coefficient is forced to be zero as well.

A B-spline curve can be written in the following way:^[7]

${\displaystyle {\textbf {C}}(\xi )=\sum _{i=1}^{n}N_{i,p}(\xi ){\textbf {B}}_{i}}$

where ${\displaystyle n}$ is the number of basis functions ${\displaystyle N_{i,p}}$, and ${\displaystyle {\textbf {B}}_{i}\in \mathbb {R} ^{d}}$ is the ${\displaystyle i^{th}}$ control point, with ${\displaystyle d}$ dimension of the space in which the curve is immersed.

An extension to the two dimensional case can be easily obtained from B-splines curves.^[7] In particular B-spline surfaces are introduced as:^[7]

${\displaystyle {\textbf {S}}(\xi ,\eta )=\sum _{i=1}^{n}\sum _{j=1}^{m}N_{i,p}(\xi )M_{j,q}(\eta ){\textbf {B}}_{i,j}}$

where ${\displaystyle n}$ and ${\displaystyle m}$ are the numbers of basis functions ${\displaystyle N_{i,p}}$ and ${\displaystyle M_{j,q}}$ defined on two different knot vectors ${\displaystyle \Xi =\{\xi _{1},\xi _{2},...,\xi _{n+p+1}\}}$, ${\displaystyle {\mathcal {H}}=\{\eta _{1},\eta _{2},...,\eta _{m+q+1}\}}$, ${\displaystyle {\textbf {B}}_{i,j}}$ represents now a matrix of control points (also called control net).

Finally, B-splines solids, that need three sets of B-splines basis functions and a tensor of control points, can be defined as:^[7]

${\displaystyle {\textbf {S}}(\xi ,\eta ,\zeta )=\sum _{i=1}^{n}\sum _{j=1}^{m}\sum _{k=1}^{l}N_{i,p}(\xi )M_{j,q}(\eta )L_{k,r}(\zeta ){\textbf {B}}_{i,j,k}}$

In IGA basis functions are also employed to develop the computational domain and not only for representing the numerical solution. For this reason they should have all the properties that permit to represent the geometry in an exact way. B-splines, due to their intrinsic structure, are not able to generate properly circular shapes for instance.^[1] In order to circumvent this issue, non-uniform rational B-splines, also known as NURBS, are introduced in the following way:^[1]

${\displaystyle R_{i}^{p}(s)={\frac {N_{i,p}(s)\omega _{i}}{W(s)}}}$

where ${\displaystyle N_{i,p}}$ is a one dimensional B-spline, ${\displaystyle W(s)=\sum _{i=1}^{n}N_{i,p}(s)\omega _{i}}$ is referred to as weighting function, and finally ${\displaystyle \omega _{i}}$ is the ${\displaystyle i^{th}}$ weight.

Following the idea developed in the subsection about B-splines, NURBS curve are generated as follows:^[1]

${\displaystyle {\textbf {C}}(s)=\sum _{i=1}^{n}R_{i}^{p}(s){\textbf {B}}_{i}}$

with ${\displaystyle {\textbf {B}}_{i}}$ ${\displaystyle i^{th}}$ vector of control points.

The extension of NURBS basis functions to manifolds of higher dimensions (for instance 2 and 3) is given by:^[1]

${\displaystyle R_{i,j}^{p,q}(s,t)={\frac {N_{i,p}(s)M_{j,q}(t)\omega _{i,j}}{\sum _{a=1}^{n}\sum _{b=1}^{m}N_{a,p}(s)M_{b,q}(t)\omega _{a,b}}}}$

${\displaystyle R_{i,j,k}^{p,q,r}(s,t,w)={\frac {N_{i,p}(s)M_{j,q}(t)L_{k,r}(w)\omega _{i,j,k}}{\sum _{a=1}^{n}\sum _{b=1}^{m}\sum _{c=1}^{l}N_{a,p}(s)M_{b,q}(t)L_{c,r}(w)\omega _{a,b,c}}}}$

There are three techniques in IGA that permit to enlarge the space of basis functions without touching the geometry and its parametrization.^[1]

The first one is known as knot insertion (or h-refinement in the FEA framework), where ${\displaystyle {\overline {\Xi }}=\{{\overline {\xi _{1}}}=\xi _{1},{\overline {\xi _{2}}},...,{\overline {\xi _{n+m+p+1}}}=\xi _{n+p+1}\}}$ is obtained from ${\displaystyle \Xi =\{\xi _{1},\xi _{2},...,\xi _{n+p+1}\}}$ with the addition of more knots, which implies an increment of both the number of basis functions and control points.^[1]

The second one is called degree elevation (or p-refinement in the FEA context), which permits to increase the polynomial order of the basis functions.^[1]

Finally the third method, known as k-refinement (without a counterpart in FEA), derives from the preceding two techniques, i.e. combines the order elevation with the insertion of a unique knot in ${\displaystyle \Xi }$.^[1]

• GeoPDEs: a free software tool for Isogeometric Analysis based on Octave
• MIG(X)FEM: a free Matlab code for IGA (FEM and extended FEM)
• PetIGA: A framework for high-performance Isogeometric Analysis Archived 2014-07-14 at the Wayback Machine based on PETSc
• G+Smo (Geometry plus Simulation modules): a C++ library for isogeometric analysis, aiming at the seamless integration of Computer-aided Design (CAD) and Finite Element Analysis (FEA), maintained by an open-source community of contributors.
• FEAP: a general purpose finite element analysis program which is designed for research and educational use, developed at University of California, Berkeley
• Bembel: An open-source isogeometric boundary element library for Laplace, Helmholtz, and Maxwell problems written in C++
• T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs: "Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement", Computer Methods in Applied Mechanics and Engineering, Elsevier, 2005, 194 (39-41), pp.4135-4195.