Isogeometric analysis: an overview and computer implementation aspects
Vinh Phu Nguyena,1,∗, Stéphane P.A. Bordasa,2 , Timon Rabczukb,3
arXiv:1205.2129v2 [cs.NA] 27 Sep 2013
a School
of Engineering, Institute of Mechanics and Advanced Materials, Cardiff University, Queen’s Buildings, The Parade, Cardiff
CF24 3AA
b Institute of Structural Mechanics, Bauhaus-Universität Weimar, Marienstraße 15 99423 Weimar
Abstract
Isogeometric analysis (IGA) represents a recently developed technology in computational mechanics that offers the
possibility of integrating methods for analysis and Computer Aided Design (CAD) into a single, unified process. The
implications to practical engineering design scenarios are profound, since the time taken from design to analysis is
greatly reduced, leading to dramatic gains in efficiency. The tight coupling of CAD and analysis within IGA requires
knowledge from both fields and it is one of the goals of the present paper to outline much of the commonly used
notation. In this manuscript, through a clear and simple Matlab R implementation, we present an introduction to IGA
applied to the Finite Element (FE) method and related computer implementation aspects. Furthermore, implementation of the extended IGA which incorporates enrichment functions through the partition of unity method (PUM) is
also presented, where several examples for both two-dimensional and three-dimensional fracture are illustrated. The
open source Matlab R code which accompanies the present paper can be applied to one, two and three-dimensional
problems for linear elasticity, linear elastic fracture mechanics, structural mechanics (beams/plates/shells including
large displacements and rotations) and Poisson problems with or without enrichment. The Bézier extraction concept
that allows FE analysis to be performed efficiently on T-spline geometries is also incorporated. The article includes a
summary of recent trends and developments within the field of IGA.
Keywords: isogeometric analysis, linear elasticity, Matlab R , NURBS, finite elements, partition of unity, enrichment,
3D cracks, boundary conditions, CAD, large displacements and rotations, shells
∗ Corresponding
author
1 nguyenpv@cardiff.ac.uk
2 stephane.bordas@alum.northwestern.edu
3 timon.rabczuk@uni-weimar.de
Preprint submitted to Elsevier
September 30, 2013
1. Introduction
1.1. Underlying concepts of isogeometric analysis
The predominant technology that is used by CAD to represent complex geometries is the Non-Uniform Rational
B-spline (NURBS). This allows certain geometries to be represented exactly that are only approximated by polynomial
functions, including conic and circular sections. There is a vast array of literature focused on NURBS (e.g. [1], [2])
and as a result of several decades of research, many efficient computer algorithms exist for their fast evaluation and
refinement. The key concept outlined by Hughes et al. [3] was to employ NURBS not only as a geometry discretisation
technology, but also as a discretisation tool for analysis, attributing such methods to the field of ‘Isogeometric Analysis’
(IGA). Since this seminal paper, a monograph dedicated entirely to IGA has been published [4] and applications can
now be found in several fields including structural mechanics, solid mechanics, fluid mechanics and contact mechanics.
We give in this section an overview of some of these recent developments while outlining the benefits and present
shortcomings of IGA. It should be emphasized that the idea of using CAD technologies in finite elements dates back
at least to [5, 6] where B-splines were used as shape functions in FEM. In addition, similar methods which adopt
subdivision surfaces have been used to model shells [7]. We also review some of the recent attempts at simplifying the
CAD-FEA integration by separating boundary and domain discretisations.
1.2. Applications
In contact formulations using conventional geometry discretisations, the presence of faceted surfaces can lead to
jumps and oscillations in traction responses unless very fine meshes are used. The benefits of using NURBS over
such an approach are evident, since smooth contact surface are obtained, leading to more physically accurate contact
stresses. Recent work in this area includes [8, 9, 10, 11, 12].
IGA has also shown advantages over traditional approaches in the context of optimisation problems [13, 14, 15, 16]
where the tight coupling with CAD models offers an extremely attractive approach for industrial applications. Another
attractive class of methods include those that require only a boundary discretisation, creating a truly direct coupling
with CAD. Isogeometric boundary element methods for elastostatic analysis were presented in [17, 18], demonstrating
that mesh generation can be completely circumvented by using CAD discretisations for analysis.
Shell and plate problems are another field where IGA has demonstrated compelling benefits over conventional
approaches [19, 20, 21, 22, 23, 24, 25]. The smoothness of the NURBS basis functions allows for a straightforward
2
construction of plate/shell elements. Particularly for thin shells, rotation-free formulations can be easily constructed
[20, 26]. Note that for multi-patch NURBS surfaces, rotation-free IGA elements require special treatment at patch
boundaries where the basis functions are found to be C 0 continuous. Furthermore, isogeometric plate/shell elements
exhibit much less pronounced shear-locking compared to standard FE plate/shell elements. Elements with smooth
boundaries such as circular and cylindrical elements were successfully constructed using the IGA concept [27, 28].
The smoothness of NURBS basis functions is attractive for analysis of fluids [29, 30, 31] and for fluid-structure
interaction problems [32, 33]. In addition, due to the ease of constructing high order continuous basis functions, IGA
has been used with great success in solving PDEs that incorporate fourth order (or higher) derivatives of the field
variable such as the Hill-Cahnard equation [34], explicit gradient damage models [35] and gradient elasticity [36]. The
high order NURBS basis has also found potential applications in the Kohn-Sham equation for electronic structure
modeling of semiconducting materials [37].
NURBS provide advantageous properties for structural vibration problems [38, 39, 40, 41] where k-refinement
(unique to IGA) has been shown to provide more robust and accurate frequency spectra than typical higher-order FE
p-methods. Particularly, the optical branches of frequency spectra, which have been identified as contributors to Gibbs
phenomena in wave propagation problems (and the cause of rapid degradation of higher modes in the p-version of
FEM), are eliminated. However when lumped mass matrices were used, the accuracy is limited to second order for any
basis order. High order isogeometric lumped mass matrices are not yet available. The mathematical properties of IGA
were studied in detail by Evans et al.[42].
The isogeometric concept has also spread to the field of meshfree methods such as [43, 44] in which spline-based
meshfree methods were presented. Industrial applications of IGA have been presented in [45, 18] along with applications
in experimental mechanics [46] where NURBS-based DIC (Digital Image Correlation) was shown to outperform standard
FE DIC.
1.3. Shortcomings of NURBS and alternative geometry discretisations
NURBS are ubiquitous in CAD but are known to exhibit major shortcomings from a computational geometry
standpoint. Perhaps the greatest difficulty encountered is the inability of NURBS to produce watertight geometries,
often complicating mesh generation. From an analysis perspective, the tensor product structure of NURBS proves to
be inefficient, caused by the global nature of refinement operations. In turn, this leads to inefficient error estimation
3
and adaptivity algorithms. One solution which has gathered momentum from both the computational geometry and
analysis communities is the use of T-splines [47] which overcome the limitations of NURBS while retaining the familiar
structure of NURBS algorithms. T-splines correct the deficiencies of NURBS by creating a single patch, watertight
geometry which can be locally refined and coarsened. Buffa et al. [48] note that linear independence of the T-spline
basis functions is not guaranteed on generic T-meshes leading to the definition of analysis-suitable T-splines [49], a
mildly restricted subset of T-splines which meet the demands of both design and analysis. Utilisation of T-splines in
an IGA framework has been illustrated in [50, 51], and by adopting a Bézier extraction process, Scott et al.[52] showed
that T-splines can be incorporated efficiently into existing FE codes.
Alternatives to T-splines include polycube splines [53], PHT-splines [54] and LR-splines [55]. PHT-splines (polynomial spline over hierarchical T-meshes) have been extended to rational splines and applied in [56, 57] to problems
in elasticity for continua and thin structures. Adaptive refinement with PHT-splines is particularly simple. Although
T-splines allow for local adaptive refinement, the complexity of knot insertion under adaptive refinement is complex,
particularly in 3D. However, we note that research is currently being pursued on hierarchical T-spline refinement
algorithms that address this issue.
Another direction of IGA research includes hierarchical B-splines [58, 59, 60] and unstructured Powell-Sabin splines
[61]. The hierarchical B-splines finite cell method [59] furnishes a seamless CAD-FEA integration for very complex
geometries. We refer also to [62] for IGA combined with finite element based local refinement capabilities. Different
subdivision surface techniques (Catmull-Clark, Loop) have also been utilized for solid and shell modeling [63, 64].
In computer aided geometric design, patching multiple NURBS parameterizations to form complex topologies is
far from trivial if certain continuity requirements are to be maintained. Trimming techniques provide a promising
alternative for representing complex NURBS domains. In [65], a trimmed surface based analysis framework has been
proposed where NURBS-enhanced FEM [66, 67] was applied to define a suitable integration domain within parameter
space. In a recent contribution [68], the authors presented an alternative method to handle trimmed NURBS geometries.
1.4. Discontinuities and fracture
IGA has been applied to cohesive fracture [69], outlining a framework for modeling debonding along material
interfaces using NURBS and propagating cohesive cracks using T-splines. The method relies upon the ability to specify
the continuity of NURBS and T-splines through a process known as knot insertion. As a variation of the eXtended Finite
4
Element Method (XFEM) [70], IGA was applied to Linear Elastic Fracture Mechanics (LEFM) using the partition of
unity method (PUM) to capture two dimensional strong discontinuities and crack tip singularities efficiently [71, 72].
The method is usually referred to as XIGA (eXtended IGA). In [73] an explicit isogeometric enrichment technique
was proposed for modeling material interfaces and cracks exactly. Note that this method is contrary to PUM-based
enrichment methods which define cracks implicitly.
A phase field model for dynamic fracture was presented in [74] using adaptive T-spline refinement to provide an
effective method for simulating fracture in three dimensions. In [75] high order B-splines were adopted to efficiently
model delamination of composite specimens and in [76], an isogeometric framework for two and three dimensional
delamination analysis of composite laminates was presented where the authors showed that using IGA can significantly
reduce the usually time consuming pre-processing step in generating FE meshes (solid elements and cohesive interface
elements) for delamination computations. A continuum description of fracture using explicit gradient damage models
was also studied using NURBS [35].
1.5. Alternatives to IGA
Other techniques which integrate CAD and analysis include the use of subdivision surfaces to model shells [7],
NURBS-enhanced finite elements [66, 67] and NURBS for BEM shape optimisation [77]. Immersed boundary methods
[78] (and references therein), the finite cell method [59, 79] and the structured XFEM [80, 81] are yet other alternatives
which aim to combine analysis and design technologies. In general, it is found that for CAD and analysis technologies to
work seamlessly together, the underlying discretisation must either be directly compatibile or easily converted between
the two.
IGA has offered significant advances towards the goal of a unified design and analysis framework, but much research
is still needed before this goal is realised. There are several indications of the future promise of IGA for industrial
design but ultimately, the litmus test of success for IGA will be whether the approach is widely adopted by industry.
1.6. Computational aspects
Some major computational aspects of IGA which have been studied so far include (i) locking issues, (ii) sensitivity
to mesh distortion, (iii) impact of high continuity of NURBS on direct solvers, (iv) collocation methods, (v) competing demands of analysis and computational geometry discretisations, (vi) construction of trivariate solids from given
bivariate surface representations and (vii) optimal quadrature rules.
5
(i) Although the smoothness of NURBS basis functions reduces to some extent the locking phenomena for constrained
problems such as incompressible media, thin-walled structures, NURBS-based FEs are not locking free [82, 83,
84, 85, 86]. Existing locking removal techniques in standard FEMs were successfully adapted to IGA such as the
Discrete Strain Gap method [84], the F/B-bar method [83, 85], the enhanced assumed strain method [86] and
mixed formulations [87].
(ii) The effect of mesh distortion on the performance of IGA for solid mechanics was discussed in [88] in which it was
found that higher-order NURBS functions are able to somewhat alleviate the impact of the distortions.
(iii) The high order continuity offered by NURBS has a negative impact on the performance of direct solvers as pointed
out in [89]. The authors found that for a fixed number of unknowns and basis degree, a higher degree of continuity
drastically increases the CPU time and RAM needed to solve the problem when using a direct solver.
(iv) In an attempt to compete with low order FEs with one-point quadrature that are extensively used in industrial
applications, isogeometric collocation methods were developed [90, 91].
(v) Due to the fact that meshes in an isogeometric framework are defined by the parametrisation of the object of
interest, the quality of the geometry parametrisation plays an important role in ensuring mesh quality. This issue
has, however, been addressed by only a few researchers [92, 93, 94, 95, 96]. In particular, in [94], the authors
proposed the concept of “analysis-aware geometry modeling”.
(vi) In CAD, solids are defined as boundary surfaces in which the interior is not explicitly modeled. In FEA, a solid
representation is necessary and therefore, the transition from CAD solids to FEA solids demands a step in which
the CAD representations are converted to solid FEA representations. Initial developments have been reported in
[97, 98, 99, 100, 101, 102]. Note that in this regard, the isogeometric boundary element method (IGABEM) can
be considered a truly isogeometric method [17, 18, 103, 104] since BEM analysis requires only the definition of a
boundary discretisation, completely defined by CAD.
(vii) Gaussian quadrature is not optimal for IGA. Research is currently focussed on optimal integration techniques
such as that in [105, 106] in which (nearly) optimal quadrature rules have been presented.
6
1.7. Available implementations
Some implementation aspects of IGA were reported in [4] and more recently, an open source IGA Matlab R code was
described in [107] with a restriction to 2D scalar PDEs. An excellent open source IGA code written in Matlab R is given
in [108]. Incorporating IGA within an object-oriented C++ FE code was discussed in [109]. Implementation details
for enriched formulations within an IGA framework are reported in [110] using commercial FE software. An IGA BEM
code written in Matlab R was presented in [104]. Isogeometric analysis was also incorporated into FEAP [111, 87]. A
high performance IGA code was given in [112] which is based on PETSc, the Portable, Extensible Toolkit for Scientific
Computation. The Python programming language has also been adopted to implement IGA e.g., [113, 114].
1.8. Contributions and outline
In this paper, we present in detail an isogeometric (Bubnov-)Galerkin finite element method applied to two- and
three-dimensional elasto-static solid/structural mechanics problems and traction-free crack problems. The discussion
is confined to NURBS for the sake of simplicity. Although no new fundamental findings are presented the contribution
of the paper are
• an overview of IGA, applications and recent developments are presented;
• implementation details for 1D/2D/3D solid mechanics and rotation-free plate elements are provided;
• implementation for 2D and 3D XIGA for traction-free cracks;
• visualization techniques for (X)IGA are discussed;
• techniques to impose Dirichlet boundary conditions including the penalty method, the Lagrange multipliers
method, the least square projection method are implemented;
• different techniques to model discontinuities in the context of IGA are discussed;
where some implementations are not available in the literature.
The paper is structured as follows: Section 2 outlines B-spline and NURBS technology used to construct curves,
surfaces and solids; the use of NURBS for discretisation within a finite element framework is treated in Section 3;
implementation of an isogeometric finite element method for two-dimensional elasticity is described in Section 4 and an
7
extended isogeometric formulation based on the concept of the PUM is the subject of Section 5. A detailed description
of our IGA Matlab code is given in Section 6 and Application of IGA to structural mechanics problems is presented in
Section 7. Numerical examples including 2D/3D fracture mechanics and large deformation shell problems are given in
Section 8.
1.9. Notation
We use lowercase indices to indicate a local index and uppercase indices to indicate a global index. We denote
dp and ds as the number of parametric directions and spatial directions respectively. Boldfont is used to indicate
matrices and vectors where the number of components is implied.
2. A brief introduction to B-splines/NURBS
As a precursor to Section 3, an understanding of the discretisation technology which underpins IGA is essential, since
the beneficial properties which are found through analysis emanate directly from the underlying CAD basis functions.
We give a brief outline of parametric functions, then state the terminology and basis function definitions which allow
curves, surfaces and solids to be represented by B-splines and NURBS. The algorithms which define B-splines and
NURBS in the present Matlab R code are based on those given in [1] in which explicit implementations are given for
all algorithms.
2.1. Parametric representation of geometry
What is fundamental to all the discretisation technology used in the present paper (and the majority of technology
in the CAGD community), is the representation of geometry through parametric functions. These define a mapping
from a given parameter to the desired geometry. We can imagine that as we move through the parameter domain, the
parametric function ‘sweeps’ out the desired shape. For example, if we consider the case of a circle of radius 1, the
implicit form of this equation is given by
x2 + y 2 = 1
(1)
or alternatively, if we define a mapping f : [0, 2π] → R2 , the parametric form of the same circle is given by
f (t) = (cos t, sin t) .
8
(2)
This is found to be much more conducive for graphical implementation. To see this, consider the task of plotting
the circle through both Eq. (1) and Eq. (2). In the case of the latter, by simply determining t at a discrete set of
points in the interval [0, 2π], the desired result is obtained. In addition, many algorithms which perform geometrical
transformations become much simpler when parametric functions are used. Both B-splines and NURBS are based on
parametric functions.
2.2. B-splines
The demands of Computer Aided Geometrical Design place certain requirements on the types of discretisation that
can be used to represent geometrical objects where, for example, in the case of car design, surfaces of G2 continuity4
are required to avoid unnatural reflections. Additional requirements include local control of geometrical features, the
ability to apply refinement algorithms and numerical stability of high order curves. It is found that representations
of objects using Bézier curves or Lagrange polynomials do not meet many of these requirements and alternatives are
sought. The CAGD community have settled on the use of B-spline-based technology since it is found to provide many
of the desired properties for interactive geometrical design, realised through the properties of the underlying B-spline
basis functions.
B-splines can be considered as a mapping from parametric space Ω̂ ⊂ Rdp to physical space Ω ⊂ Rds 5 . In this sense,
a B-spline can be considered to ‘sweep’ out a curve, surface or volume as we move through the range of parameter
values. In the present work we define coordinates in parameter space as ξ = (ξ, η, ζ) = (ξ 1 , ξ 2 , ξ 3 ) and coordinates in
physical space as x = (x, y, z) = (x1 , x2 , x3 ). These are simplified accordingly in the case of one- and two-dimensional
problems. What remains is to determine the particular form of the mapping x : Ω̂ → Ω.
To construct a B-spline, a knot vector must be specified and is defined as an ordered set of increasing parameter
values Ξ = {ξ1 , ξ2 , . . . , ξn+p+1 }, ξi ≤ ξi+1 where ξi is the i th knot, n is the number of basis functions and p is the
polynomial order. The knot vector divides the parametric space into intervals usually referred to as knot spans with
the number of coincident knots for a particular knot value referred to as a knot with a certain multiplicity k. That is,
a knot has a multiplicity k if it is repeated k times in the knot vector. Most commonly, open knot vectors are used
where the first and last knots have a multiplicity k = p + 1. Appendix A outlines alternative knot vector notations
B outlines the differences between C k and Gk continuity
discussion of these spaces in given in Section 3.1.
4 Appendix
5 Further
9
which remove redundant information, but it can be assumed that the above knot vector notation is used throughout
the present work.
An important property of B-splines formed from open knot vectors is that they are interpolatory at their start
and end points, greatly facilitating the imposition of boundary conditions for analysis. In geometrical terms, this
corresponds to a B-spline that coincides with its start and end control points. A uniform knot vector is associated to
evenly distributed knots. Otherwise it is classified as a non-uniform knot vector.
Knot vectors are not commonly used by CAD designers, and in most CAD software the ability to modify knot
vectors is not provided. It is much more common to tailor the geometry through modification of the polynomial order,
control points and weightings.
2.2.1. B-spline basis functions
Given a knot vector Ξ, the associated set of B-spline basis functions {Ni,p }ni=1 are defined recursively by the
Cox-de-Boor formula, starting with the zeroth order basis function (p = 0)
1 if ξi ≤ ξ < ξi+1 ,
(3)
ξi+p+1 − ξ
ξ − ξi
Ni,p−1 (ξ) +
Ni+1,p−1 (ξ).
ξi+p − ξi
ξi+p+1 − ξi+1
(4)
Ni,0 (ξ) =
and for a polynomial order p ≥ 1
Ni,p (ξ) =
0 otherwise,
in which fractions of the form 0/0 are defined as zero.
Some salient properties of B-spline basis functions include: (1) they constitute a partition of unity, (2) each basis
function is non-negative over the entire parametric domain, (3) they are linearly independent, (4) the support of a
B-spline basis function Ni,p is given by p + 1 knot spans denoted by [ξi , ξi+p+1 ], (5) they exhibit C p−mi continuity
across knot ξi where mi is the multiplicity of knot ξi and (6) B-spline basis functions in general are not interpolatory,
equivalent to saying that the Kronecker-delta property is not guaranteed at control points. This last point requires
careful treatment when imposing non-homogeneous Dirichlet boundary conditions, and is discussed later in Section 4.2.
To demonstrate the effect of alternative knot vectors, Fig. 1 illustrates a set of quadratic (p = 2) B-spline basis
10
functions for a uniform knot vector. Fig. 2 illustrates a corresponding set of basis functions for an open, non-uniform
knot vector. Of particular note is the interpolatory nature of the basis function at each end of the interval created
through an open knot vector, and the reduced continuity at ξ = 4 due to the presence of the location of a repeated
knot where C 0 continuity is attained. Elsewhere, the functions are C 1 continuous (C p−1 ). The ability to control
continuity by means of knot insertion is particularly useful for modeling discontinuities such as cracks or material
interfaces [69, 76].
Figure 1: Quadratic B-spline basis functions defined for the uniform knot vector Ξ = {0, 1, 2, 3, 4, 5, 6}.
During implementation, the derivatives of B-splines are required in both CAD and analysis to compute quantities
such as tangent and normal vectors and field variable derivatives. The first derivative of a B-spline basis function is
computed recursively from lower order basis functions as
d
p
p
Ni,p−1 (ξ) −
Ni+1,p−1 (ξ).
Ni,p (ξ) =
dξ
ξi+p − ξi
ξi+p+1 − ξi+1
(5)
Derivatives of higher order can be found in [1].
2.3. Refinement algorithms
A feature which is essential for both computational geometry and analysis is the ability to successively refine
discretisations to allow intricate geometries to be modelled or to capture rapid variations in the field solution. A
11
Figure 2:
Quadratic B-spline basis functions defined for the open, non-uniform knot vector Ξ =
{0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5}. Note the flexibility in the construction of basis functions with varying degrees of regularity.
significant advantage of adopting B-spline-based functions as a discretisation tool is the ability to apply a variety of
refinement algorithms in a simple manner. We restrict ourselves to a brief discussion of the available algorithms, with
a more detailed treatment given in [4].
B-spline refinement algorithms in the context of isogeometric analysis can be assigned to one of three types:
1. Knot insertion is the process whereby additional knots are inserted into the knot vectors thereby creating
additional knot intervals or elements, in the context of analysis. This is directly analogous to h-refinement seen
in the conventional FEM.
2. Degree elevation is the process of raising the order of the underlying basis, directly analogous to p-refinement
in the conventional FEM.
3. k-refinement is unique to IGA and consists of a combined process of degree elevation and knot insertion. These
processes are not commutative and therefore the order in which these refinements are applied will change the
final basis. k-refinement first applies degree elevation proceeded by knot insertion, offering a reduction in degrees
of freedom over its counterpart [3].
12
2.4. B-spline curves
Given a set of B-spline basis functions {NA,p }nA=1 and a set of control points {PA }nA=1 where PA ∈ Rds , we can
construct a piecewise-polynomial B-spline curve as
C(ξ) =
n
X
PA NA,p (ξ).
(6)
A=1
An example quadratic B-spline curve is shown in Fig. 3 constructed from a uniform open knot vector. Its associated
control polygon, constructed by piecewise linear interpolation of the control points is also shown. A B-spline curve
inherits all of the continuity properties of its underlying basis and, as mentioned previously, the use of open knots
ensures that the first and last points are interpolatory. Note that the tangent to the curve at the first control point is
defined by the first two control points. This property is exploited in rotation-free element formulations for thin-walled
structures (section 7).
4
3.5
3
2.5
2
1.5
1
0.5
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Figure 3: A quadratic (p = 2) B-spline curve with a uniform open knot vector Ξ = {0, 0, 0, 1, 2, 3, 4, 5, 5, 5}. Control
points are denoted by filled circles and the control polygon is denoted by the red line. Note that the curve is at least
C 1 continuous everywhere except at the extremities where, due to the presence of repeated knots, C −1 continuity is
obtained.
13
2.5. B-spline surfaces and volumes
Given two knot vectors for each parametric direction Ξ1 = {ξ1 , ξ2 , . . . , ξn+p+1 } and Ξ2 = {η1 , η2 , . . . , ηm+q+1 }, and
a control net Pi,j ∈ Rds , a tensor-product B-spline surface is defined as
S(ξ, η) =
n X
m
X
Ni,p (ξ)Mj,q (η)Pi,j ,
(7)
i=1 j=1
where Ni,p (ξ) and Mj,q (η) are the univariate B-spline basis functions of order p and q corresponding to knot vectors
Ξ1 and Ξ2 , respectively. Defining a global index as
A = n(j − 1) + i,
(8)
Eq. (7) can be rewritten in a more compact form as
S(ξ) =
n×m
X
PA NAp,q (ξ),
(9)
A=1
in which NAp,q is a bivariate B-spline basis function defined as NAp,q (ξ) = Ni,p (ξ)Mj,q (η). Figs. 4 illustrates an example
bicubic B-spline basis function.
The extension to B-spline volumes is straightforward, where a trivariate basis is formed through a tensor product
of B-spline basis functions as
V(ξ, η, ζ) =
l
m X
n X
X
Ni,p (ξ)Mj,q (η)Lk,r (ζ)Pi,j,k ,
(10)
i=1 j=1 k=1
or, by defining a global index A through
A = (n × m)(k − 1) + n(j − 1) + i,
14
(11)
0.7
0.8
0.9
1
r
0.4
0.35
0.5
0.6
0.3
0.25
0.2
0
1
0.9
m
0.8
0.7
0.6
0.5
0.4
0.3
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.1
0.3
0.1
0.05
0.2
0.4
0.15
0.2
0.1
0
0.2
0.1
0
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.7
0.6
0.5
0.4
0.3
n
0.2
0.1
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 4: A bivariate cubic B-spline basis function for knot vectors Ξ1 = Ξ2 = {0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1}.
15
a simplified form of Eq. (10) can be written as
V(ξ) =
n×m×l
X
PA NAp,q,r (ξ).
(12)
A=1
To compute derivatives of B-spline surfaces and volumes, the chain rule is applied to (5) with a detailed treatment
of this topic provided in [1].
2.6. NURBS
B-splines are convenient for free-form modelling, but they lack the ability to exactly represent some simple shapes
such as circles and ellipsoids. This is why today, the de facto standard technology in CAD is a generalisation of
B-splines referred to as NURBS (Non-Uniform Rational B-Splines). NURBS are formed through rational functions
of B-splines, forming a superset of B-splines. They inherit all the favourable properties of B-splines and are favoured
over their counterpart due to their ability to form exact representations of conic sections such as spheres, ellipsoids,
paraboloids and hyperboloids. In addition, there exist efficient algorithms for their evaluation and refinement.
2.6.1. NURBS basis functions
NURBS basis functions are defined as
Ri,p (ξ) =
Ni,p (ξ)wi
Ni,p (ξ)wi
= Pn
W (ξ)
î=1 Nî,p (ξ)wî
(13)
where {Ni,p }ni=1 is the set of B-spline basis functions of order p and {wi }ni=1 , wi > 0 is the set of NURBS weights.
Selecting appropriate weights permits the description of many different types of curves including polynomials and
circular arcs. For the special case in which all weights are equal, the NURBS basis reduces to the B-spline basis.
NURBS weights for certain simple geometries are given in [1], but in general, weights are user-defined through CAD
packages such as Rhino [115].
The first derivative of a NURBS basis function Ri,p is computed using the quotient rule as
′
(ξ)W (ξ) − Ni,p (ξ)W ′ (ξ)
Ni,p
d
Ri,p (ξ) = wi
,
dξ
W (ξ)2
16
(14)
′
(ξ) ≡
where Ni,p
d
dξ Ni,p (ξ)
and
W ′ (ξ) =
n
X
′
Nî,p
(ξ)wî .
(15)
î=1
2.6.2. NURBS curves, surfaces and volumes
In a similar fashion to B-spline curves, the NURBS curve associated with a set of control points and weights
{PA , wA }nA=1 and basis functions {RA,p }nA=1 is defined as
C(ξ) =
n
X
PA RA,p (ξ).
(16)
A=1
NURBS surfaces are constructed from a linear combination of bivariate NURBS basis functions, control points Pi,j ∈
Rds and weights wi,j > 0 as
S(ξ, η) =
m
n X
X
p,q
Pi,j Ri,j
(ξ, η),
(17)
i=1 j=1
where the bivariate NURBS basis functions are defined as
N (ξ)Mj (η)wi,j
Pim
.
î=1
ĵ=1 Nî (ξ)Mĵ (η)wî,ĵ
p,q
Ri,j
(ξ, η) = Pn
(18)
Alternatively, using the mapping defined by Eq. (8), Eq. (17) can be written more succinctly as
S(ξ) =
n×m
X
p,q
PA R A
(ξ).
(19)
A=1
NURBS volumes are constructed from control points Pi,j,k ∈ Rds , weights wi,j,k > 0 as
V(ξ, η, ζ) =
l
m X
n X
X
p,q,r
Pi,j,k Ri,j,k
(ξ, η, ζ)
(20)
i=1 j=1 k=1
p,q,r
where the trivariate NURBS basis functions Ri,j,k
are given by
Ni (ξ)Mj (η)Pk (ζ)wi,j,k
.
Pm Pl
î=1
ĵ=1
k̂=1 Nî (ξ)Mĵ (η)Pk̂ (ζ)wî,ĵ,k̂
p,q,r
Ri,j,k
(ξ, η, ζ) = Pn
17
(21)
Using the mapping given by Eq. (11), Eq. (20) can also be written as
V(ξ) =
n×m×l
X
p,q,r
PA R A
(ξ)
(22)
A=1
Remark 2.1. As mentioned in Section 1.6, CAD representations are usually composed of surface models or boundaryrepresentations. Trivariate discretisations defined by (22) are not normally explicitly given, and therefore some preprocessing is required before domain based numerical methods such as the FEM can be applied. Except for the case of
simple cases such as extruded-surface models and swept models, this task is far from trivial. This is currently an open
research topic in IGA.
3. NURBS as a basis for analysis: isogeometric finite element formulation
Our attention now focusses on the use of B-splines and NURBS as a discretisation tool for analysis, outlining the
core concepts of isogeometric analysis. In this section the important spaces and mappings are defined, followed by the
isogeometric FEM formulation in which we use NURBS as a basis for analysis. The discussion is made more concrete
through a one-dimensional example.
3.1. Relevant spaces
Familiarity must be gained with the spaces that are commonplace in isogeometric analysis and the relationships
that exist between each. Those that are considered presently in the context of B-splines and NURBS include: index,
parameter, physical and parent space.
3.1.1. Index space
Index space is formed through the specified knot vectors by giving each knot value a distinct coordinate, regardless
of whether the knot is repeated or not. As an example, consider a NURBS patch defined through bivariate NURBS
basis functions with knot vectors Ξ1 = {0, 0, 0, 1, 2, 3, 3, 3}, Ξ2 = {0, 0, 1, 1} in each of the parametric directions ξ, η
respectively. This will form the index space as illustrated in Fig. 5 where the presence of repeated knots leads to several
regions of zero parametric area. Index space is often used during implementation6 , discarding elements of non-zero
6 T-splines
being one notable example, with the local basis function mesh directly analogous to index space
18
non-zero parametric area
Figure 5: Creation of index space from knot vectors Ξ1 = {0, 0, 0, 1, 2, 3, 3, 3}, Ξ2 = {0, 0, 1, 1} with non-zero parametric
area highlighted.
parametric area, but in the present work we choose to only consider those elements that have a non-zero parametric
area, obviating the need for index space.
3.1.2. Parametric space
Parametric space (sometimes referred to as the ‘pre-image’ of the NURBS mapping) is formed by considering
only the non-zero intervals between knot values. For the knot vectors considered previously, the parametric space
is illustrated in Fig. 6 which can subsequently be reduced to a unit square through appropriate normalisation. All
parametric spaces can be reduced to a unit interval (dp = 1), square (dp = 2) or cube (dp = 3) in this manner. We
define the parametric space as Ω̂ ⊂ Rdp with a associated set of parametric coordinates ξ = (ξ, η, ζ) = (ξ 1 , ξ 2 , ξ 3 ) ∈ Ω̂
(dp = 3). If normalisation is performed, Ω̂ = [0, 1]dp .
Fig. 6 also reveals that regions bounded by knot lines with non-zero parametric area lead to a natural definition of
element domains. More formally, a set S ⊂ Ξ of unique knot values can be defined as
S = {ξ1 , ξ2 , · · · , ξns }
ξi 6= ξi+1 for 1 ≤ i ≤ ns − 1
(23)
where ns is the number of unique knot values. This is generalised to Si ⊂ Ξi which represents the unique knot values
19
for each parametric direction i = 1, 2, . . . dp . Elements can now be defined in the general multivariate case as
Ω̂e = [ξi , ξi+1 ] ⊗ [ηj , ηj+1 ] ⊗ [ζk , ζk+1 ]
1 ≤ i ≤ n1s − 1,
(24)
1 ≤ j ≤ n2s − 1,
1 ≤ k ≤ n3s − 1,
ξ i ∈ S 1 , η j ∈ S2 , ζ k ∈ S3
where n1s , n2s and n3s represent the number unique knots in the ξ,η and ζ parametric directions respectively. This leads
to a natural numbering scheme for elements over a patch as
e = k(n2s − 1)(n1s − 1) + j(n1s − 1) + i.
(25)
Eq. (24) and (25) can be simplied accordingly for dp = 1, 2.
3.1.3. Physical space
The B-spline and NURBS mappings of Eqs. (9) and (17) transform coordinates in parameter space to physical space
Ω ⊂ Rds . For three-dimensional domains, we associate a coordinate system x = (x, y, z) = (x1 , x2 , x3 ) for physical
space, which appropriate modifications for one- and two-dimensional problems. Fig. 7 illustrates a NURBS mapping
for the parametric space shown in Fig. 6 for an arbitrary set of control points and weights. The control grid (which
defines the connectivity between control points) is also shown. The non-interpolatory nature of control points in the
interior of the domain is evident, and represents a notable difference over conventional Lagrangian meshes.
3.1.4. Parent space
The previous three spaces are inherent to B-splines and NURBS, but for analysis to be performed we require the
definition of an additional space, commonly referred to as parent space Ω̃ = [−1, 1]dp . This is required for the use of
numerical integration routines which are often defined over the interval [−1, 1]. Parent space coordinates are denoted
˜ η̃, ζ̃) = (ξ˜1 , ξ˜2 , ξ˜3 ) with corresponding simplifications for dp = 1, 2.
as ξ̃ = (ξ,
20
unnormalised
normalised
Figure 6: Parametric space defined by non-zero knot intervals. All parametric spaces defined for a B-spline or NURBS
patch can be normalised to a unit interval, unit square, unit cube in 1D, 2D, 3D respectively. Knot lines provide a
natural definition of element boundaries.
Figure 7: A 2D NURBS surface defined for knot vectors Ξ1 = {0, 0, 0, 1, 2, 3, 3, 3}, Ξ2 = {0, 0, 1, 1}. The control mesh
is shown in red with control points denoted by black circles. Knot lines shown in blue indicate element boundaries.
21
3.2. Isogeometric formulation
Before outlining the details of isogeometric analysis in a FE context, it is instructive to consider the similarities and
differences of IGA over conventional discretisation technology. Lagrangian basis functions are most commonly used
to discretise both the geometry and unknown fields in an isoparametric fashion. In this way, exactly the same basis
functions are used for both. For the majority of cases, the geometry is always approximated incurring a geometrical
error which may lead to erroneous results, particularly for highly oscillatory problems. In addition, once a discretisation
is generated from a given CAD model, the geometry information that has been lost can never be retrieved, forming a
one-way process from discretisation to analysis. This has serious repercussions for efficiency, attenuated by the iterative
nature of design.
Isogeometric analysis also makes use of an isoparametric formulation, but a key difference over its Lagrangian
counterpart is the use of basis functions generated by CAD to discretise both the geometry and unknown fields. Not
only is task of discretisation (meshing) greatly reduced or eliminated entirely, but a direct link is made with CAD,
forming a bi-directional process. In addition, the use of CAD discretisations ensures that the exact geometry is used
at all stages of analysis, incurring no geometrical error. But the key compelling feature of IGA is the unified nature of
design and analysis.
We can summarise these points as:
• Conventional finite element analysis: the basis which is chosen to approximate the unknown field is also
used to approximate the known geometry. This most commonly takes the form of (low order) Lagrangian basis
functions. In most cases the geometry is only approximated. CAD and analysis are disparate.
• Isogeometric analysis: the basis is generated by CAD which captures the geometry exactly. This basis is
also used to approximate the unknown field. Refinement may be required for the unknown fields, but the exact
geometry is maintained at all stages of analysis. CAD and analysis are combined to form a unified process.
3.3. Isogeometric discretisation
The B-spline and NURBS discretisations outlined in Section 2 are written in terms of parametric coordinates, but
to use such discretisations for analysis, we must provide a mapping that allows us to operate at the parent element
level. We will outline the appropriate mappings in Section 3.4, but for now let us assume that B-spline and NURBS
22
basis functions can be written in terms of parent coordinates. This allows us to state the isoparametric discretisation
used to approximate both the geometry and fields in IGA. For a given element e, the geometry is expressed as
xe (ξ̃) =
nen
X
Pea Rae (ξ̃)
(26)
a=1
where a is a local basis function index, nen = (p + 1)dp is the number of non-zero basis functions over element e and
Pea ,Rae are the control point and NURBS basis function associated with index a respectively. We employ the commonly
used notation of an element connectivity mapping [116] which translates a local basis function index to a global index
through
A = IEN(a, e).
(27)
Global and local control points are therefore related through PA ≡ PIEN(a,e) ≡ Pea with similar expressions for Rae . A
field u(x) which governs our relevant PDE can also be discretised in a similar manner to (26) as
e
u (ξ̃) =
nen
X
dea Rae (ξ̃)
(28)
a=1
where dea represents a control (nodal) variable. In contrast to conventional discretisations, these coefficients are not
in general interpolatory at nodes. This is similar to the case of meshless methods built on non-interpolatory shape
functions such as the moving least squares (MLS) [117, 118, 119].
3.4. Mappings (change of variables)
The use of NURBS basis functions for discretisation introduces the concept of parametric space which is absent in
conventional FE implementations. The consequence of this additional space is that an additional mapping must be
performed to operate in parent element coordinates. As shown in Fig. 8, two mappings are considered for IGA with
NURBS: a mapping φ̃e : Ω̃ → Ω̂e and S : Ω̂ → Ω. The mapping xe : Ω̃ → Ωe is given by the composition S ◦ φ̃e .
Taking the case dp = ds = 2, an element defined by Ω̂e = [ξi , ξi+1 ] ⊗ [ηi , ηi+1 ] is mapped from parent space to
23
parent space
parameter space
physical space
Figure 8: Diagrammatic interpretation of mappings from parent space through parametric space to physical space.
24
parametric space through
φ̃e (ξ̃) =
1 [(ξi+1 − ξi )ξ˜ + (ξi+1 + ξi )]
2
1 [(ηj+1 − ηj )η̃ + (ηj+1 + ηj )]
(29)
2
with an associated Jacobian determinant given by
|Jξ̃ | =
1
(ξi+1 − ξi )(ηj+1 − ηj ).
4
(30)
Similarly, the mapping from parametric space to physical space is given by the NURBS expressions of Eqs. (16),
(19) and (22). In the case dp = ds = 2, the Jacobian of transformation for this mapping is represented by the matrix
∂x
∂ξ
Jξ =
∂y
∂ξ
in which the components are calculated as
∂x
∂η
∂y
∂η
(31)
n
en
∂Rae (ξ)
∂x X
Pea
=
.
∂ξ
∂ξ
a=1
(32)
The associated Jacobian determinant is denoted by |Jξ |.
The mapping xe : Ω̃ → Ωe , formed through the composition of the previous mappings, can be written as
e
x (ξ̃) =
=
N
X
PA RA φ̃ (ξ̃)
A=1
nen
X
e
Pa Ra (φ̃ (ξ̃))
a=1
=
e
nen
X
Pea Rae (ξ̃)
!
!
(33)
e
(34)
e
(35)
a=1
(36)
where N is the number of global basis functions and (·)|e denotes that the expression (·) is restricted to element e. The
25
Jacobian determinant for this mapping is given by
|J| = |Jξ ||Jξ̃ |.
(37)
With this final mapping and Jacobian determinant, it is possible to integrate a function f : Ω → R over the physical
domain as
Z
f (x) dΩ =
Ω
=
=
=
nel Z
X
e=1 Ω
nel Z
X
e=1 Ω̂
nel Z
X
e=1
nel Z
X
e=1
f (x) dΩ
(38)
f (x(ξ)) |Jξ | dΩ̂
(39)
e
e
f (x(φ̃e (ξ̃))) |Jξ | |Jξ̃ | dΩ̃
(40)
Ω̃
f (ξ̃) |J|dΩ̃
(41)
Ω̃
with the final integral in a suitable form for application of standard Gauss-Legendre quadrature (hereafter named
Gaussian quadrature). As detailed in [116] for Lagrangian basis functions, a rule of (p+1)×(q +1) Gaussian quadrature
can be applied for two-dimensional elements in which p and q denote the orders of the chosen basis functions in the ξ
and η direction. The same procedure is also used for NURBS basis functions in the present work, although it should
be emphasised that Gaussian quadrature is not optimal for IGA. Research is currently focussed on optimal integration
techniques such as that in [105, 106] in which an optimal quadrature rule, known as the half-point rule, has been
applied.
Spatial derivatives of basis functions are also required for element assembly algorithms, and are calculated as
∂Rae
∂x
=
∂Rae
∂ξ
∂ξ
∂x
∂η
∂x
with the derivatives ∂Rae /∂ξ obtained through Eq. (14).
26
∂ξ
e
∂y
= ∂Ra J−1 ,
ξ
∂ξ
∂η
∂y
(42)
3.5. One-dimensional IGA formulation
To illustrate the use of the isogeometric concept in a finite element context, a one-dimensional IGA formulation is
developed through an example. Efforts are made to maintain the notation adopted in [116] and [4]. In the following,
we define the domain Ω ⊂ R with boundary Γ ≡ ∂Ω. The boundary is partitioned as Γ = ΓD ∪ ΓN , ΓN ∩ ΓD = ∅
where ΓD and ΓN denote the Dirichlet and Neumann boundaries respectively with an overline representing a closed
set. Robin boundary conditions are not considered in the present work.
Letting Ω = (0, 1) and ΓD = Γ = {0, 1} with ΓN = ∅, we seek the solution u : Ω → R such that
d2 u(x)
+ b(x) = 0,
dx2
x ∈ Ω,
(43)
with the Dirichlet boundary conditions specified through g : ΓD → R as
g(0) = 0,
g(1) = 0.
(44)
Choosing b(x) = x, it can be shown that the exact solution to Eq. (43) subject to Eq. (44) is given by
1
1
u(x) = − x3 + x.
6
6
(45)
To construct the Galerkin formulation of the preceding problem, the infinite dimensional spaces U and V which
represent the usual trial and test spaces respectively7 are required. By multiplying Eq. (43) by a test function w ∈ V
and applying integration by parts, the weak form of the problem reads: find u ∈ U , such that
Z
Ω
dw du
dΩ =
dx dx
Z
wb dΩ
Ω
for all w ∈ V .
(46)
This is reduced to a finite-dimensional problem by creating finite-dimensional subspaces U h ⊂ U and V h ⊂ V both
formed through a NURBS basis. Letting g h ∈ U denote the finite-dimensional representation of g, the solution we
7 Specifically,
U = {u : u ∈ H 1 (Ω), u = g on ΓD }, V = {w : w ∈ H 1 (Ω), w = 0 on ΓD } where H 1 (Ω) represents a Sobolev space of order
1
27
wish to seek can be written as
uh = v h + g h
(47)
which is valid for all v h ∈ V h . Expression (47) underpins a key concept of the Galerkin formulation in that the solution
uh and test (weighting) function wh are constructed from the same space of functions. We also see that the boundary
condition data is ‘built in’ to the solution by the inclusion of the term g h .
The approximation to our solution can be written in terms of NURBS basis functions as
h
u =
neq
X
dB RB (x) +
nnp
X
gB RB (x)
(48)
B=neq +1
B=1
where neq is the number of equations (number of unknowns) and nnp is the total number of nodal points (total degrees
n
n
eq
np
of freedom) in the system. {dB }B=1
represents the set of all unknown control variables and {gB }B=n
is the set of
eq +1
known Dirichlet control variables. Likewise, the test function is discretised as
h
w =
neq
X
cA RA (x)
(49)
A=1
Substitution of Eqs. (48) and (49) into Eq. (46) and noting that the values of the set {cA } are arbitrary, the discrete
form is now written as
neq neq Z
X
X
dRA dRB
dΩ
dx
Ω dx
A=1 B=1
!
dB =
neq Z
X
A=1
Ω
RA b dΩ −
neq
X
nnp
X
A=1 B=neq +1
Z
Ω
dRA dRB
dΩ gB
dx dx
(50)
or, by defining the element stiffness matrix and element force vector as
dRA dRB
dΩ
dx
Ω dx
Z
nnp
X
RA b dΩ −
FA =
KAB =
Z
Ω
B=neq +1
(51)
Z
Ω
dRA dRB
dΩ,
dx dx
(52)
(50) can expressed in matrix notation as
Kd = F
28
(53)
where
K = [KAB ]
(54)
d = {dB }
(55)
F = {FA }
A, B = 1, 2, . . . neq .
(56)
Eq. (53) is in a form amenable for computation.
b2
b1
b3
b4
curve
1
n1
n4
0.8
n2
n3
0.6
0.4
0.2
0
0
0.5
1
basis
e1
xi1
xi
e2
xi2
xi3
para
Figure 9: One dimensional isogeometric analysis example: exact geometry, quadratic basis functions and mesh in the
parametric space. IEN(:,1)=[1 2 3] and IEN(:,2)=[2 3 4].
29
3.5.1. Assembly of system of equations
The global stiffness matrix K and global force vector F are assembled by looping over each element (non-zero knot
interval) and inserting the appropriate element entries into their relevant rows and columns in the global system of
equations. To allow integration over each element, a parent coordinate system is adopted making use of the mappings
and Jacobian matrices outlined in Section 3.4. For example, for a particular element e with local basis function indices
a, b = 1, 2, . . . , p + 1, we can compute a local element stiffness matrix keab using parent coordinates as
keab =
Z
Ω̃
dRae dRbe
|J| dΩ̃
dx dx
(57)
where use has been made of expressions (42) and (37). This can be inserted into the global stiffness matrix through
the element mapping given by Eq. (27).
To make this discussion more concrete, the element assembly process for a simple quadratic (p = 2) NURBS
discretisation with Ξ = {0, 0, 0, 0.5, 1, 1, 1} is outlined. We apply the discretisation to the problem defined by Eqs. (43)
and (44). Fig. 9 illustrates the relevant NURBS discretisation including control points, basis functions and element
definitions. Note that there are p + 1 non-zero basis functions over an element e given by {Rae }p+1
a=1 . We can therefore
summarise the element data for this discretisation as
element
knot interval
non-zero basis functions
control variables
control points
1
[ξ3 , ξ4 ]
R1 , R2 , R3
d 1 , d 2 , d3
P1 , P2 , P 3
2
[ξ4 , ξ5 ]
R2 , R3 , R4
d 2 , d 3 , d4
P2 , P3 , P 4
(58)
The implementation required to compute the element stiffness matrix is shown in Box 1. The code follows the same
standard FEM routines presented in [120] where the notation ||·|| is used to signify the L2 norm8 .
Fig. 10 illustrates the solution obtained for a NURBS discretisation for both two and four elements. Convergence
to the exact solution is observed, verifying the present IGA implementation.
In order to assess how high order B-spline elements behave for problems with localized gradients, let us consider
8 That
is, for a vector x = (x1 , x2 , . . . , xn ), ||x|| = (x21 + x22 + . . . + x2n )1/2
30
Box 1 Element stiffness matrix evaluation for element Ω̂e = [ξi , ξi+1 ]
1. P = [P1 ; P2 ; P3 ]
2. sctr = [IEN(1, e) IEN(2, e) IEN(3, e)]
3. Set ke = 0
4. Loop over Gauss points (GPs) {ξ˜j , w̃j }
j = 1, 2, . . . , ngp ∗
(a) Compute parametric coordinate ξ = φ̃e (ξ˜j ) (see Eq. (29))
e
(b) Compute derivatives Ra,ξ
(a = 1, 2, 3) at point ξ
e
e
e
(c) Define vector Rξ = [R1,ξ
R2,ξ
R3,ξ
]
(d) Compute |Jξ | = ||Rξ P||
(e) Compute |Jξ̃ | = 0.5(ξi+1 − ξi )
(f) Compute shape function derivatives Rx = Jξ−1 RT
ξ
(g) ke = ke + Rx RT
x |Jξ | |Jξ̃ | w̃j
5. End loop over GPs
6. Assemble ke into the global matrix as K(sctr, sctr) = K(sctr, sctr) + ke
w̃j denotes a Gauss point weight and ngp is the total number of Gauss points.
0.07
0.06
0.06
0.05
0.05
0.04
0.04
u
0.07
u
∗
0.03
0.03
0.02
0.02
numerical solution
exact solution
0.01
numerical solution
exact solution
0.01
0
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
x
0
0
1
(a) Two element NURBS mesh
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
x
1
(b) Four element NURBS mesh.
Figure 10: Quadratic NURBS solution for one-dimensional Poisson equation.
31
the following problem [119]
u,xx (x) + b(x) = 0 x ∈ [0, 1];
u(0) = 0,
u(1) = 1,
(59)
with
b(x) =
2α2 − 4[α2 (x − 0.5)]2 exp −[α(x − 0.5)]2
0
x ∈ [0.42, 0.58]
(60)
otherwise
The exact solution of this problem is
u(x) = x + exp −[α(x − 0.5)]2 x ∈ [0, 1]
(61)
We use a value of 50 for α and the exact solution exhibits a sharp peek at location x = 0.5. We are going to solve this
problem using elements of order ranging from one (linear elements) to five (quintic elements). To remove error in the
numerical integration of the body force term, Eq. (??), 10 GPs were used for each element. We use the k-refinement (to
be discussed in detail in Section 6.5) in building meshes of different basis orders. The initial mesh consists of one single
linear element with knot vector Ξ = {0, 0, 1, 1}. The parametrization is thus linear and after performing k-refinement,
the parametrization is still linear. Therefore, a point with x = 0.3 corresponds to ξ = 0.3 in the parameter space. The
Matlab file for this problem is iga1DStrongGradient.m in the iga folder of our source code.
Figs. 11a,b shows the results obtained with meshes consisting of 16 and 32 elements of C p−1 continuity where p
denotes the B-spline basis order. It is obvious that smooth basis is not suitable to problems with sharp gradients.
Linear elements which have a C 0 continuity at location x = 0.5 (precisely at every knots) give better results than high
order B-spline elements.
Knot insertion was made to insert the value 0.5 p times to the initial knot vector so that the basis is C 0 continuous
at x = 0.5. For example for p = 2, a vector {0.5, 0.5} was inserted to the knot vector. The results are given in Fig.
11c,d where the peak was captured better. Finally 0.42, 0.5, 0.58 were added to the knots p times so that the basis
are C 0 continuous at locations x = 0.42, 0.5, 0.58. The corresponding results are given in Fig. 12. With only 16 cubic
elements (25 CPs) the exact solution was well captured. This example showed the flexibility of B-splines/NURBS–high
32
1.5
1
u
u
1
1.5
exact solution
p=1
p=2
p=3
p=4
p=5
0.5
0
0
0.5
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
x
0
0
1
(a) 16 elements
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
x
1
(b) 32 elements
1.5
1.6
exact solution
p=1
p=2
p=3
p=4
1.4
1.2
1
exact solution
p=1
p=2
p=3
p=4
u
u
1
exact solution
p=1
p=2
p=3
p=4
p=5
0.8
0.5
0.6
0.4
0
0
0.2
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
x
1
0.4 0.420.440.460.48 0.5 0.520.540.560.58 0.6
x
(c) 32 elements with C 0 continuity at x = 0.5
(d) 32 elements with C 0 continuity at x = 0.5 (close up)
Figure 11: Comparison of the IGA result against the exact solution for the one dimensional PDE given in Eq. (59).
33
1.6
exact solution
p=1
p=2
p=3
1.4
1.2
1.6
exact solution
p=1
p=2
p=3
1.4
1
u
1.2
0.8
1
u
0.6
0.4
0.8
0.6
0.2
0.4
0
0
0.2
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
x
(a) 16 elements with
C0
1
0
0.42
0.44
0.46
0.48
0.5
0.52
0.54
0.56
x
continuity at x = 0.42, 0.5, 0.58
(b) 16 elements with
up)
C0
continuity at x = 0.42, 0.5, 0.58 (close
Figure 12: Comparison of the IGA result against the exact solution for the one dimensional PDE given in Eq. (59).
order functions and any level of continuity can be easily achieved.
4. Elasticity: two-dimensional implementation
The application of the FEM to elasticity is common and represents a familiar language to many researchers.
We therefore outline a two-dimensional implementation of IGA for linear elasticity, highlighting the differences over
conventional discretisations.
4.1. Assembly process for two dimensional elastostatic analysis
Consider a domain Ω, bounded by Γ. The boundary is partitioned into two sets: Γu and Γt with displacements
prescribed on Γu and tractions t prescribed on Γt : Γ = Γt ∪ Γu , Γt ∩ Γu = ∅. The weak form of a linear elastostatics
problem is to find u in the trial space 9 , such that for all test functions δu in the test space
Z
9 contains
10 contains
ε(u) : D : ε(δu)dΩ =
Ω
Z
Γt
C 0 functions
C 0 functions vanishing on Γu
34
t · δudΓ +
Z
Ω
b · δudΩ,
10
,
(62)
where the elasticity matrix is denoted by D and b refers to a body force. Using the Galerkin method where the same
shape functions Ra (ξ̃) are used for both u and δu, we can write
u(x) =
nnp
X
RA (ξ̃)uA ,
δu(x) =
nnp
X
RA (ξ̃)δuA ,
(63)
A=1
A=1
where uA , δuA denote the nodal displacement and its variations, respectively and nnp is the total number of control
points. In 2D, each control point has two unknowns–the x and y displacements, hence one writes uA = {uxA , uyA }.
Proper modification can be made for 3D problems.
Substitution of these approximations into Eq. (62) and using the arbitrariness of the nodal variations gives the
standard discrete set of equations K u = f with
KAB =
Z
Ω
BT
A DBB dΩ,
fA =
Z
RA tdΓ +
Γt
Z
RA bdΩ,
A, B = 1, 2, . . . , nnp .
(64)
Ω
In two dimensions, the strain-displacement matrix BA is given by
RA,x
BA =
0
RA,y
0
RA,y
,
RA,x
(65)
where the shape function derivatives are computed according to Eq. (42) and RA,x ≡ dRA /dx.
We now consider a concrete two dimensional problem shown in Fig. 13. In this case, the knot vectors are Ξ1 =
[0, 0, 1, 1] and Ξ2 = [0, 0, 0, 0.5, 1, 1, 1], respectively. The orders of the basis functions are p = 1 and q = 2. There two
control points (n = 2) along the ξ direction and four control points (m = 4) along the η direction that results in 8
control points (nnp = n × m = 8).
The domain and non-zero basis functions for Element 1 are given by
direction
knot interval
non zero basis
ξ
[ξ2 , ξ3 ]
N 1 , N2
η
[η3 , η4 ]
M1 , M2 , M3
35
(66)
a
b
b24
2
b22
2
b23
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
1
1
b14
0.8
0.8
b13
0.6
0.6
b12
0.4
0.4
0.2
0.2
b11
b21
0
−2
−1.5
−1
−0.5
0
−2
0
−1.5
−1
−0.5
0
c
1
8
N42
eta3
7
12
6
eta2
0.5
N32
5
N22
11
4
N12
3
0
0
1
eta1
0.2
0.4
0.6
0.8
1
2
xi2
xi1
1
N11
N21
0
0
1
Figure 13: Two dimensional isogeometric analysis example: (a) exact geometry, (b) mesh in physical space and (c)
mesh in the parametric space. The two rectangles are used to illustrate the control points belonging to each element.
36
Hence there are six non-zero basis functions on element e = 1 which can be assembled into a vector R as follows
R1 = [N1 M1 , N2 M1 , N1 M2 , N2 M2 , N1 M3 , N2 M3 ].
(67)
These six basis functions are associated with six global basis indices given by (in Matlab R notation)
IEN(:, 1) = [1, 2, 3, 4, 5, 6].
(68)
Similarly, for element 2, the shape function vector is given by
R2 = [N1 M2 , N2 M2 , N1 M3 , N2 M3 , N1 M4 , N2 M4 ],
(69)
IEN(:, 2) = [3, 4, 5, 6, 7, 8].
(70)
with the associated global indices
Control points are stored in a two dimensional matrix of dimensions nnp × 2. The connectivity data is stored in a two
dimensional matrix of dimensions (p + 1) ∗ (q + 1) × nel where nel denotes the number of elements. For the example
under consideration, these two matrices are given by
controlPts = B11
B21
B12
B22
B13
B23
B14
B24
T
,
T
1
2
3
4
5
6
IEN =
.
3 4 5 6 7 8
(71)
The knot intervals along the ξ and η directions are stored in the following matrices
elRangeU = 0 1 ,
0
elRangeV =
0.5
0.5
,
1
(72)
where the number of rows is equal to the number of elements in each direction (n1s − 1, n2s − 1 respectively).
With this vector of basis functions R, the control points associated to this element, we can compute the derivatives
37
of the basis functions. The derivatives of the basis functions with respect to x are stored in the following vector
R,x = R1,x
R2,x
R3,x
R4,x
R5,x
R6,x
T
,
(73)
Similarly, we have R,y for the derivatives of the basis functions with respect to y. Having these basis function derivatives,
we are now ready to define the B matrix for any element e
R[1],x
Be =
0
R[1],y
0
R[2],x
0
R[1],y
0
R[2],y
R[1],x
R[2],y
R[2],x
···
···
,
···
(74)
where the control point displacement vector is stored in the following order u = [ux1 , uy1 , ux2 , uy2 , . . . , uxnnp , uynnp ]T .
The element stiffness matrix is then given by
Ke =
Z
Ωe
BT
e DBe dΩe ,
(75)
which is then assembled to the global stiffness matrix using the element connectivity matrix and the fact that a control
point I corresponds to positions 2 ∗ I − 1 and 2 ∗ I in the global displacement vector u.
For implementation convenience, Box 2 gives a procedure of an isogeometric analysis of 2D elasticity problems. Note
that this aims for a Matlab implementation and the whole elementary stiffness matrix is computed in one go. In order
to compute the external force vector, it is convenient to define a boundary mesh as shown in Fig. 14. The computation
R
of the external force vector Γ RT t̄dΓ then follows the one dimensional assembly procedure given in Section 3.5.
4.2. Boundary condition enforcement
Fig. 15 illustrates two kinds of Dirichlet boundary conditions: on edge AD, ux = 0 and on edge BC, uy = ū. The
former is a homogeneous Dirichlet boundary condition (BC) while the latter is referred to as uniform inhomogeneous
Dirichlet BCs. Homogeneous Dirichlet BCs can be enforced by setting the corresponding control variables as zeros
(in this example, setting uxI = 0, I = 2, 4, 6, 8). For edge BC, the inhomogeneous Dirichlet BCs can also be satisfied
by setting uyI = ū, I = 1, 3, 5, 7. This is due to the partition of unity property of the NURBS basis i.e., uBC
=
y
M1 (η)uy1 + M2 (η)uy3 + M3 (η)uy5 + M4 (η)uy7 = (M1 (η) + M2 (η) + M3 (η) + M4 (η))ū = ū. Note that a boundary of a
38
Box 2 Procedure for isogeometric analysis of 2D elasticity problems.
1. Loop over elements, e = 1, . . . , nel
(a) Determine NURBS coordinates [ξi , ξi+1 ] × [ηj , ηj+1 ] using elRangeU and elRangeV
(b) Get connectivity array, sctr = IEN(:, e)
(c) Define sctrB(1, 1 : 2 : 2 ∗ nn) = 2 ∗ sctr − 1; sctrB(1, 2 : 2 : 2 ∗ nn) = 2 ∗ sctr for assembly∗
(d) Ke = 0
(e) Loop over Gauss points, {ξ̃ j , w̃j }
i.
ii.
iii.
iv.
v.
Compute
Compute
Compute
Compute
Compute
j = 1, 2, . . . , ngp ∗∗
ξ corresponding to ξ̃ j (Eq. 29)
|Jξ̃ | (Eq. 30)
derivatives of shape functions R,ξ and R,η at ξ
Jξ using controlPts(sctr, :), R,ξ and R,η (Eq. 31)
Jacobian inverse J−1
ξ and determinant |Jξ |
vi. Compute derivatives of shape functions R,x = [R,ξ R,η ]J−1
ξ (Eq. 42)
vii. Use R,x to build the strain-displacement matrix B (Eq. 74 )
viii. Compute Ke = Ke + w̃j |Jξ̃ ||Jξ |BT DB
(f) End loop over Gauss points
(g) Assemble Ke : K(sctrB, sctrB) = K(sctrB, sctrB) + Ke
2. End loop over elements
∗
nn denotes the number of control points per element i.e., nn=length(sctr).
w̃j denotes a Gauss point weight and ngp is the total number of Gauss points.
∗∗
39
12
11
12
2
8
7
8
1
3
10
4
4
6
2
9
5
1
Figure 14: Boundary mesh for external force computation. A linear basis is used in the η direction. Assume a traction
is applied on the edge containing nodes 4, 8 and 12. The boundary mesh is composed of two linear isogeometric
elements 1 and 2 .
NURBS surface is a NURBS curve due to the use of open knots. Inhomogeneous Dirichlet BCs applied on corner control
points (black points in Fig. 15) are enforced by simply setting the corner control variables equal to the prescribed
values since the NURBS shape functions at these points satisfy the Kronecker delta property (assuming the use of open
knot vectors). This is called direct imposition of Dirichlet BCs.
For cases other than the ones previously discussed such as a prescribed displacement imposed at interior control
point 3 or a non-uniform Dirichlet BC–Dirichlet BCs that vary from point to point applied on edge BC, special
treatment of Dirichlet BCs have to be employed as is the case for meshless methods. Techniques available include
the Lagrange multiplier method, the penalty method, the augmented Lagrangian method and we refer to [119] for
an overview of these techniques in the context of meshless methods. In [121] a transformation method was proposed
to impose inhomogeneous Dirichlet BCs in IGAFEM. However this method requires modifications to the stiffness
matrix which breaks the usual structure of a FE code. The authors in [71] presented a weak enforcement of general
inhomogeneous Dirichlet BCs using a least squares minimization and the implementation is described in what follows.
The same procedure was used in imposing BCs in meshless methods when coupling a fluid to a solid domain through a
master-slave concept [122]. Imposing Dirichlet boundary conditions with Nitsche’s method was presented in [123] for
spline-based finite elements.
40
1
N22
N12
N42
N32
0.8
0.6
0.4
0.2
0
0
0.5
1
eta3
eta2
eta1
d
a
8
6
4
2
b
c
5
1
7
3
Figure 15: Imposing Dirichlet BCs: black points denote corner control points where the NURBS basis satisfies the
Kronecker delta property.
41
The basic idea of the least squares method is to find the parameters of the boundary control points that minimize
the following quantity
J=
1X
2
||u(xC ) − ū(xC )||
2
C
1X X
=
RA (xC )qA − ū(xC )
2
C
2,
(76)
A
where xC denotes a set of collocation points distributed on the essential boundary Γu , qA are the parameters of the
control points defining Γu , ū(x) represents the prescribed displacements and RA represents the NURBS basis functions
that are non-zero at xC which are univariate NURBS functions RA (ξ) for the boundary of a NURBS surface is a
NURBS curve. For the sake of clarity, let us consider the case where there is only one collocation point and a quadratic
basis (thus there are 3 non-zero RA at xC ). So, we have
J=
1
2
||R1 (xC )q1 + R2 (xC )q2 + R3 (xC )q3 − ū(xC )|| .
2
(77)
The partial derivatives of J with respect to qi are given by
∂J
= [R1 (xC )q1 + R2 (xC )q2 + R3 (xC )q3 − ū(xC )]R1 (xC )
∂q1
∂J
= [R1 (xC )q1 + R2 (xC )q2 + R3 (xC )q3 − ū(xC )]R2 (xC )
∂q2
∂J
= [R1 (xC )q1 + R2 (xC )q2 + R3 (xC )q3 − ū(xC )]R3 (xC ).
∂q3
The condition
∂J
∂q
(78)
= 0 thus gives the following linear system
R1 R1
R R
1 2
R1 R3
R2 R1
R2 R2
R2 R3
R3 R 1
R3 R 2
R3 R 3
x
q1
xC
q x
2
q3x
q1y
ūx (xC )R1 (xC )
=
ūx (xC )R2 (xC )
y
q3
ūx (xC )R3 (xC )
q2y
ūy (xC )R1 (xC )
ūy (xC )R2 (xC )
.
ūy (xC )R3 (xC )
(79)
By collecting all the NURBS basis at point xC in a column vector N(xC ), the control points displacements in the x
42
and y directions in qx and qy , respectively, Eq. (79) can be written in a compact form as
[N(xC )NT (xC )]qx = ūx (xC )N(xC )
(80)
[N(xC )NT (xC )]qy = ūy (xC )N(xC ).
Repeating the same analysis for other collocation points xC on the Dirichlet boundary, one obtains the linear system
Aq = b with two different b (one for the x component and the other for the y component). The dimension of A
is nD × nD where nD denotes the number of control points defining the Dirichlet boundary. Having these boundary
control point displacements, the enforcement of Dirichlet BCs (when solving Ku = f ) are then treated as in standard
FEM.
We note that this procedure involves only control points that define the essential boundary. This is in sharp contrast
to meshless shape functions such as the MLS used in the Element Free Galerkin method [124] where the displacements
at a point on the essential boundary depend not only on the nodes on that boundary but also the neighbouring interior
nodes. It is said that NURBS therefore satisfy the so-called weak Kronecker delta property as local Maxent interpolants
approximations [125].
Listing 1 gives the Matlab R implementation of the least-squares method. Note that in our implementation, the
collocation points are uniformly distributed in the parameter space. We refer to [71] for a discussion on the influence
of the collocation points on the accuracy. We refer to Fig. 16 for an illustration of this method.
Listing 1: Matlab R implementation of the least-squares method
A
= zeros ( noDispNodes , noDispNodes ) ;
bx = zeros ( noDispNodes , 1 ) ;
by = zeros ( noDispNodes , 1 ) ;
noxC
= 4 ; % number o f c o l l o c a t i o n p t s / e l e m e n t
% l o o p o v e r bottom edge
f o r i e =1: noElemsU
sctr
= bottomEdgeMeshIGA ( i e , : ) ; % s t a n d a r d c o n n e c t i v i t y a r r a y ( g l o b a l numbering )
pts
= controlPts ( sctr , : ) ;
% c o n t r o l pts of element i e
sctrA
= bndElement ( i e , : ) ;
% c o n n e c t i v i t y array f o r A matrix
xiE
= elRangeU ( i e , : ) ;
% parametric coords of i e
43
7
y
6
2
x
5
control
1
col
1
2
3
1
4
2
Figure 16: Illustration of the implementation of the least squares method: a quadratic NURBS surface with knot
vectors Ξ1 = {0, 0, 0, 0.5, 1, 1, 1} and Ξ2 = {0, 0, 0, 0.5, 1, 1, 1}. Essential BCs are imposed on the bottom and right
edges. A one dimensional mesh for this boundary is created. The control points defining the essential boundary are
numbered from one to the total number of boundary control points (seven in this example). Note that this numbering
is required only for assembling the matrix A.
xiArr
= linspace ( xiE ( 1 ) , xiE ( 2 ) , noxC ) ; % c o l l o c a t i o n p t s x C
f o r i c =1:noxC
xi
= xiArr ( i c ) ;
[N dNdxi ]
= NURBS1DBasisDers ( x i , p , uKnot , w e i g h t s ) ;
A( sctrA , s c t r A ) = A( sctrA , s c t r A ) + N’ ∗N;
x
= N ∗ pts ;
% exact displacements
[ ux , uy ]
= e x a c t G r i f f i t h ( x , E0 , nu0 , sigmato , xTip , seg , c r a c k l e n g t h ) ;
bx ( s c t r A ) = bx ( s c t r A ) + ux∗N ’ ;
by ( s c t r A ) = by ( s c t r A ) + uy∗N ’ ;
end
end
% loop over other edges i f neccessary
% s o l v e t h e system Aq x=bx and Aq y=by
[ LL UU] = lu (A ) ;
qxTemp
= LL\bx ;
qyTemp
= LL\by ;
44
qx
= UU\qxTemp ; qy
= UU\qyTemp ;
% l a t e r , b e f o r e s o l v i n g Ku=f , qx and qy w i l l be used
t o e n f o r c e BCs a s i n c o n v e n t i o n a l FEM.
In our Matlab R code implementation are provided for the penalty method, the Lagrange multiplier method and
the least squares method. The implementation of the two former methods are considered standard and the reader is
referred to [119] for details.
5. Extended isogeometric finite element method
There are basically two ways in which discontinuities can be modeled in the context of IGA: PUM based enrichment
and knot insertion. For the former, there are works of [110, 71, 72, 73]. For the latter, we refer to [69, 76]. PUM
based methods are general for they can be applied to any kinds of discontinuity such as weak and strong discontinuities
whereas knot insertion used to produce the C −1 continuity is only suitable for cracks. Knot insertion found great
applications for delamination analyses [76] where the crack path is known a priori. Note that in [69], knot insertion
was used in combination with T-splines to model cracks of which trajectory is not known in advance. However the
implementation is tedious.
In the category of PUM based IGA methods which are sometimes called XIGA (eXtended IGA), the method in
[73] is different from existing XIGA formulations [110, 71, 72] in a way that discontinuities (holes/inclusions/cracks)
are exactly isogeometrically represented by NURBS. Note that in most of XIGA techniques, discontinuities are defined
implicitly using level set method. This section briefly presents the XIGA formulation of [71, 72] and implementation
aspects are deferred to Section 6.7.
The extended finite element method (XFEM) (see e.g., [70], [126] for a recent review and open XFEM library
[127]) is a local PUM enrichment method in which internal boundaries such as holes, inclusions, cracks are modeled
independently of the FE discretisation which allows for crack growth modeling without remeshing. Two dimensional
extended isogeometric finite element formulation (XIGA) was presented in [71, 72] in which the displacement field is
enriched for traction-free crack modelling using the following approximation
h
u (x) =
X
I∈S
RI (x)uI +
X
RJ (x)H(x)aJ +
J∈Sc
X
K∈Sf
45
RK (x)
4
X
α=1
Bα b α
K
!
,
(81)
where RI,J,K are the NURBS basis functions. In addition to the standard degrees of freedom (dofs) uI , additional
c
dofs aJ and bα
K are introduced. The set S includes the control points/nodes whose support is cut by the crack and
the set Sf are control points whose support contains the crack tip xtip , see Fig. 17. Note that we use a topological
tip enrichment and in the literature another tip enrichment scheme called geometrical enrichment with a fixed area (to
ensure that the role of enrichment in the approximation space does not vanish as the mesh is refined) is present see
e.g., [128]. The Heaviside function H is given by
H(x) =
+1
if (x − x∗ ) · n ≥ 0
,
(82)
−1 otherwise
where x∗ is the projection of point x on the crack and n denotes the outward normal vector to the crack. And the
branch functions, which span the crack tip displacement field, are given by
[B1 , B2 , B3 , B4 ] (r, θ) =
√
√
θ √
θ √
θ
θ
r sin , r cos , r sin cos θ, r cos cos θ ,
2
2
2
2
(83)
where r and θ are polar coordinates in the local crack front coordinate system (see e.g., [70] for details). It is noted
that B1 is discontinuous along the crack face.
hea
ver
D
C
P
N
A
B
Q
M
con
bra
Figure 17: Illustration of enriched node sets Sf and Sc for a quadratic NURBS mesh. The thick red line denotes the
crack.
Using the standard Galerkin procedure as outlined in Section 4.1 the discrete system of equations Ku = f are
formed by an enlarged B matrix given by
B=
B
std
46
B
enr
,
(84)
where Bstd is the standard strain-displacement matrix B (cf. Eq. (65)) and Benr is the enriched B matrix of which
components are given by
Benr
I
(RI ),x ΨI + RI (ΨI ),x
=
0
(RI ),y ΨI + RI (ΨI ),y
0
(RI ),y ΨI + RI (ΨI ),y
(RI ),x ΨI + RI (ΨI ),x
,
(85)
ΨI may represent either the Heaviside function H or the branch functions Bα depending whether control point I is
Heaviside or near tip enriched. The unknown vector u contains both displacements and enriched dofs. This extended
IGAFEM can be implemented within an available IGAFEM code with little modification following the ideas given in
the meshless review and computer implementation aspects paper [119]. Some implementation aspects will be given in
Section 6.7. In our code enrichment for holes and material interfaces are also given. They are, however, not covered
here because the implementation follows the ideas given here.
6. MIGFEM- A Matlab IGA (X)FEM code
In this section we describe shortly the open source IGA Matlab (X)FEM program which can be downloaded from
https://sourceforge.net/projects/cmcodes/. The code supports one, two and three dimensional linear elasticity
problems. Extended IGA for crack and material interface modelling is also implemented. Geometrically nonlinear solid
and structural mechanics models are available. The features of the code include:
• Global h,p and k-refinement is provided for one, two and three dimensional meshes.
• Extended IGA for 2D/3D stationary traction-free cracks and material interfaces.
• Visualization of displacements and stresses in Paraview.
• Inhomogeneous Dirichlet boundary conditions are treated with the penalty method, the Lagrange multiplier
method and the least squares method.
• Compatible multi-patch isogeometric formulation for two dimensional problems.
• Support for T-splines via the Bézier extraction operators.
47
• Structural elements including beams, plates and thin shells.
• Implicit Newmark scheme and explicit central difference scheme for time discretisation.
• Python scripts to extract Rhino3d NURBS surfaces to be used for IGA.
6.1. Data structure
MIGFEM follows the Matlab FEM code described in [120]. The main data structures include (1) element (store the
element connectivity), (2) controlPts (store control point coordinates), (3) weights (store the weights) (4) K (stiffness
matrix) and (5) f (external force vector). Contrary to FEM in which the element connectivity and nodal coordinates
are inputs which have been created by a meshing program, in IGA, the input consists of CAD data including knots,
control points, order of basis functions. Therefore, one has to construct the element matrix based on the knots and
the basis orders.
Construction of the element matrix is illustrated by a 2D example shown in Fig. 18. Given the knot vectors uKnot
and vKnot together with the orders of the basis p and q, one can compute the number of control points along the ξ and
ζ directions, denoted by n and m. Then we define a two dimensional matrix of dimension m × n called node pattern
given by for the illustrated example shown in Fig. 18 which has 4 control points along ξ and 4 control points along η
directions
1
5
node pattern =
9
13
2
3
6
7
10
11
14
15
4
8
,
12
16
(86)
which is simply an application of Eq. (8) for defining global indexes.
The number of elements along the two directions n1s − 1, n2s − 1 are noElemsU = length(unique(uKnot)) − 1
and noElemsV = length(unique(uKnot)) − 1 (in the Matlab language). The connectivity matrix for the ξ direction,
denoted by connU which is a noElemsU × (p + 1) matrix and the connectivity matrix for the η direction, denoted by
48
node
14
15
16
1
13
N42
4
3
9
10
11
12
3
2
connV
N32
2
3
1
0.5
4
2
N22
1
5
6
7
8
1
1
2
2
3
0
0
0.2
0.4
0.6
2
0.8
1
N12
1
2
1
2
3
1
4
connU
n22
n12
n32
0.8
basis
3
4
n42
xi
0.6
0.4
et
0.2
0
0
0.5
1
Figure 18: Two dimensional isogeometric analysis: mesh generation for a bi-quadratic NURBS surface (2 × 2 elements).
The circles denote the control points.
connV which is a noElemsV × (q + 1) matrix are given by
1
2
3
connU =
,
2 3 4
1
2
3
connV =
,
2 3 4
(87)
for the example under consideration. Matrix connU stores the indices of the non-zero basis functions at each knot span
for the ξ direction. Having this information and the global indexes of the NURBS basis given in Eq. (86), we are able
49
to define the connectivity matrix for the whole mesh, called element which is a noElemsU ∗ noElemsV × (p + 1)(q + 1)
matrix. For the example being considered, this matrix reads
1 2 3
5
6
7
9
10
11
2 3 4 6
7
8 10 11 12
element =
.
5 6 7 9 10 11 13 14 15
6 7 8 10 11 12 14 15 16
(88)
Note that the matrix element is the transpose of the IEN matrix. We decided to use element instead of IEN to be
compatible with the FEM code [120] on which MIGEM is built.
Finally in order to compute the mapping from the parent domain Ω̃ to the parametric space, we define the knot
intervals in two directions as
0
rangeU =
0.5
0.5
,
1
0
rangeV =
0.5
0.5
.
1
(89)
In order to retrieve the parametric coordinates of a specific element, the matrix index, that is a noElemsU ∗
noElemsV ×2 matrix, is used. For a given element e, its parametric coordinates are determined by rangeU(index(e,1),:);
rangeV(index(e,2),:).
The above discussion, together with the illustration given in Fig. 18 is implemented in Matlab in the file generateIGA2DMesh.m, located in folder meshing, and shown in Listing 2. In that folder, one can find similar M files
for generating 1D and 3D meshes.
Listing 2: Mesh generation for two dimensional problems. Input are uKnot, vKnot, and p, q.
uniqueUKnots
= u n i q u e ( uKnot ) ; uniqueVKnots
= u n i q u e ( vKnot ) ;
noElemsU
= length ( uniqueUKnots ) −1;%#o f e l e m e n t s x i d i r .
noElemsV
= length ( uniqueVKnots ) −1;%#o f e l e m e n t s e t a d i r .
noPtsX
= length ( uKnot)−p−1; noPtsY
nodePattern
= zeros ( noPtsY , noPtsX ) ;
= length ( vKnot)−q −1;
count = 1 ;
f o r i =1: noPtsY
f o r j =1: noPtsX
50
n o d e P a t t e r n ( i , j ) = count ; count
= count + 1 ;
end
end
% d e t e r m i n e our e l e m e n t r a n g e s and t h e c o r r e s p o n d i n g knot i n d i c e s a l o n g each d i r e c t i o n
[ elRangeU , elConnU ] = b u i l d C o n n e c t i v i t y ( p , uKnot , noElemsU ) ;
[ elRangeV , elConnV ] = b u i l d C o n n e c t i v i t y ( q , vKnot , noElemsV ) ;
noElems = noElemsU ∗ noElemsV ;
e l e m e n t = zeros ( noElems , ( p +1)∗( q + 1 ) ) ;
e = 1;
f o r v =1: noElemsV
vConn = elConnV ( v , : ) ;
f o r u=1: noElemsU
c = 1 ; uConn = elConnU ( u , : ) ;
f o r i =1: length ( vConn )
f o r j =1: length ( uConn )
e l e m e n t ( e , c ) = n o d e P a t t e r n ( vConn ( i ) , uConn ( j ) ) ; c = c + 1 ;
end
end
e = e + 1;
end
end
i n d e x = zeros ( noElems , 2 ) ;
count = 1 ;
f o r j =1: s i z e ( elRangeV , 1 )
f o r i =1: s i z e ( elRangeU , 1 )
i n d e x ( count , 1 ) = i ; i n d e x ( count , 2 ) = j ; count = count + 1 ;
end
end
6.2. Shape function routines
The shape function routines (evaluate NURBS basis functions and derivatives with respect to parametric coordinates) are implemented using MEX files to improve the performance. They are located in folder C files. Listing 3
gives some commonly used shape function routines to compute the NURBS basis function and their first derivatives in
51
1D, 2D and 3D at a given point. In the last line, second derivatives are also computed.
Listing 3: Shape function routines: uKnot,vKnot,wKnot store Ξ1 , Ξ2 , Ξ3 .
% R, dRdxi : row v e c t o r s i . e . dRdxi d i m e n s i o n i s ( 1 , 3 ) f o r p=2
[R dRdxi ]
= NURBS1DBasisDers ( Xi , p , uKnot , w e i g h t s ) ;
[R dRdxi dRdeta ]
= NURBS2DBasisDers ( [ Xi ; Eta ] , p , q , uKnot , vKnot , w e i g h t s ’ ) ;
[R dRdxi dRdeta dRdzeta ] = NURBS3DBasisDers ( [ Xi ; Eta ; Zeta ] , p , q , r , uKnot , vKnot , wKnot , w e i g h t s ’ ) ;
[R dRdxi dRdeta dR2dxi dR2det dR2dxe ] = NURBS2DBasis2ndDers ( [ Xi ; Eta ] , p , q , uKnot , vKnot , w e i g h t s ’ ) ;
Listing 4: Matlab code for 2D spatial derivatives of shape functions.
1
% i n p u t : e l e m e n t e and p a r a m e t e r c o o r d i n a t e s o f a Gauss p o i n t Xi and Eta
2
sctr
= element ( e , : ) ;
%
element s c a t t e r vector
3
pts
= controlPts ( sctr , : ) ;
%
o f d i m e n s i o n nn x 2 where nn : number o f nodes
4
5
6
% d e r i v a t e o f R=Nxi ∗ Neta w . r . t x i and e t a
[ dRdxi dRdeta ] = NURBS2Dders ( [ Xi ; Eta ] , p , q , uKnot , vKnot , w e i g h t s ’ ) ;
% Jacobian matrix
7
jacob
= pts ’ ∗ [ dRdxi ’ dRdeta ’ ] ;
8
invJacob
= inv ( j a c o b ) ;
9
10
% nn x 2 m a t r i x :
dRdx
f i r s t column d e r i v a t i v e s wrt x . . .
= [ dRdxi ’ dRdeta ’ ] ∗ i n v J a c o b ;
6.3. Assembly process
The assembly of an IGA-FEM code is given in Listing 5 where it can be seen that the procedure is almost identical
to that used in the conventional FEM. The minor differences lie in (1) the need of the elements in the parameter space
(lines 4 to 7) and (2) the second map (from the parent domain to the parametric domain) in the numerical integration
of the stiffness matrix (line 35).
Listing 5: Matlab code for IGA (2D elasticity problems).
1
[W,Q]= q u a d r a t u r e ( 4 , ’ GAUSS ’ , 2 ) ; % 4 x4 p o i n t q u a d r a t u r e
2
f o r e =1: noElems % Loop o v e r e l e m e n t s ( knot s p a n s )
3
idu
= index ( e , 1 ) ;
4
idv
= index ( e , 2 ) ;
52
5
xiE
= elRangeU ( idu , : ) ; % [ x i i , x i i +1]
6
etaE
= elRangeV ( idv , : ) ; % [ e t a j , e t a j +1]
7
sctr
= element ( e , : ) ;
8
nn
= length ( s c t r ) ; nn2 = 2∗ nn ;
9
s c t r B ( 1 , 1 : 2 : 2 ∗ nn ) = 2∗ s c t r −1;
%
10
s c t r B ( 1 , 2 : 2 : 2 ∗ nn ) = 2∗ s c t r
11
B
= zeros ( 3 , 2 ∗ nn ) ;
12
pts
= controlPts ( sctr , : ) ;
13
f o r gp =1: s i z e (W, 1 )
element s c a t t e r vector
;
% l o o p o v e r Gauss p o i n t s
14
pt
= Q( gp , : ) ;
15
wt
= W( gp ) ;
16
% compute c o o r d s i n p a r a m e t e r s p a c e
17
Xi
= p a r e n t 2 P a r a m e t r i c S p a c e ( xiE , pt ( 1 ) ) ;
18
Eta
= p a r e n t 2 P a r a m e t r i c S p a c e ( etaE , pt ( 2 ) ) ;
19
J2
= jacobianPaPaMapping ( xiE , etaE ) ;
% d e r i v a t e o f R=Nxi ∗ Neta w . r . t x i and e t a
20
[ dRdxi dRdeta ] = NURBS2Dders ( [ Xi ; Eta ] , p , q , uKnot , vKnot , w e i g h t s ’ ) ;
21
% Jacobian matrix
22
23
jacob
= pts ’ ∗ [ dRdxi ’ dRdeta ’ ] ;
24
J1
= det ( j a c o b ) ; i n v J a c o b
25
dRdx
= [ dRdxi ’ dRdeta ’ ] ∗ i n v J a c o b ;
26
% B matrix
27
B ( 1 , 1 : 2 : nn2 )
= dRdx ( 1 , : ) ; B ( 2 , 2 : 2 : nn2 )
= dRdx ( 2 , : ) ;
28
B ( 3 , 1 : 2 : nn2 )
= dRdx ( 2 , : ) ; B ( 3 , 2 : 2 : nn2 )
= dRdx ( 1 , : ) ;
29
% compute e l e m e n t a r y s t i f f n e s s m a t r i x and
30
K( s c t r B , s c t r B ) = K( s c t r B , s c t r B ) + B’ ∗C∗B∗ J1 ∗ J2 ∗wt ;
31
32
= inv ( j a c o b ) ;
assemble i t to the g l o b a l matrix
end
end
6.4. Post-processing
We present here a simple technique to visualize the IGA results that reuse available visualization techniques for
finite elements. To simplify the exposé, only 2D cases are considered here. In the first step, a mesh consisting of
four-noded quadrilateral (Q4) elements is generated, see Fig. 19. We call this mesh the visualization mesh (whose
53
connectivity matrix is stored in elementV and nodal coordinates are stored in node) whose nodes are images of the
knots ξi , ηj in the physical space. In the second step, quantities of interest e.g., stresses are computed at the nodes
of the Q4 mesh. This mesh together with the nodal values can then be exported to a visualization program such as
Paraview, see [129], for visualization. It should be emphasized that due to the high order continuity of the NURBS
basis, there is no need to perform nodal averaging as required in standard C 0 finite element analysis to obtain smooth
fields. Listing 6 gives the Matlab code (we have removed some code due to the similarity with the assembly code)
for building the Q4 visualization mesh, computing the stresses at the nodes of this mesh and exporting the result to
Paraview. The source code can be found in the file plotStress1.m located in folder post-processing. The results
can be then visualized directly in Matlab or exported to Paraview, see Listing 7.
Listing 6: Matlab code for computing stresses and displacements at nodes of the visualization mesh.
b u i l d V i s u a l i z a t i o n M e s h ; % b u i l d v i s u a l i z a t i o n Q4 mesh
s t r e s s = zeros ( noElems , 4 , 3 ) ;
= zeros ( noElems , 4 , 2 ) ;
disp
f o r e =1: noElems
idu
= index ( e , 1 ) ; idv
= index ( e , 2 ) ;
xiE
= elRangeU ( idu , : ) ; etaE = elRangeV ( idv , : ) ;
sctr
= element ( e , : ) ;
sctrB
= [ s c t r s c t r+noCtrPts ] ; %v e c t o r s c a t t e r s B m a t r i x
% element s c a t t e r vector
uspan = FindSpan ( noPtsX −1 ,p , xiE ( 1 ) , uKnot ) ;
vspan = FindSpan ( noPtsY −1 ,q , etaE ( 1 ) , vKnot ) ;
% l o o p o v e r 4 nodes ( i n s t e a d o f l o o p o v e r GPs)
gp = 1 ;
f o r i v =1:2
% Q4 e l e m e n t s , t h e nodes a r e numbered c o u n t e r −c l o c k w i s e
i f ( i v ==2) xiE = sort ( xiE , ’ descend ’ ) ; end
f o r i u =1:2
Xi
= xiE ( i u ) ; Eta = etaE ( i v ) ;
[N dRdxi dRdeta ]= NURBS2DBasisDersSpecial ( [ Xi ; Eta ] , p , q , uKnot , vKnot , w e i g h t s ’ , [ uspan ; vspan ] ) ;
% B matrix as usual
strain
= B∗U( s c t r B ) ;
s t r e s s ( e , gp , : )
= C∗ s t r a i n ;
gp = gp +1;
54
(a) NURBS mesh
(b) approximate Q4 mesh
(c) stress visualization on a refined mesh
Figure 19: Exact NURBS mesh (top left) and approximate Q4 mesh (top right) for visualization purpose. The nodes
in the Q4 mesh are the intersections of the ξ and η knot lines. The bottom figure shows a contour plot of a stress field
in Paraview. It should be emphasized that the mesh in (b) does not provide a sufficiently smooth contour plot for it
to be directly usable. The result given in (c) was obtained with a refined NURBS mesh (hence a refined Q4 mesh).
55
end
end
end
Listing 7: Matlab code for post-processing.
% p l o t s i g m a x x c o n t o u r d i r e c t l y i n Matlab
figure
p l o t f i e l d ( node , elementV , ’ Q4 ’ , s t r e s s ( : , : , 1 ) ) ;
% e x p o r t t o VTK f o r m a t t o p l o t i n Mayavi o r Paraview
sigmaXX = zeros ( s i z e ( node , 1 ) , 1 ) ;
sigmaYY = zeros ( s i z e ( node , 1 ) , 1 ) ;
sigmaXY = zeros ( s i z e ( node , 1 ) , 1 ) ;
dispX = zeros ( s i z e ( node , 1 ) , 1 ) ;
dispY = zeros ( s i z e ( node , 1 ) , 1 ) ;
f o r e =1: s i z e ( elementV , 1 )
c o n n e c t = elementV ( e , : ) ;
f o r i n =1:4
nid = connect ( in ) ;
sigmaXX ( n i d ) = s t r e s s ( e , in , 1 ) ;
sigmaYY ( n i d ) = s t r e s s ( e , in , 2 ) ;
sigmaXY ( n i d ) = s t r e s s ( e , in , 3 ) ;
end
end
VTKPostProcess ( node , elementV , 2 , ’ Quad4 ’ , ’ result . vtu ’ , . . .
[ sigmaXX sigmaYY sigmaXY ] , [ dispX dispY ] ) ;
For three-dimensional problems, the same procedure is used where a mesh of tri-linear brick elements is created and
the values of interest are computed at the nodes of this mesh (see file plotStress3d.m). The results are then exported
to Paraview under a structured grid format (*.vts files), see the file mshToVTK.m).
56
6.5. h, p, k-refinement
For the refinement of NURBS, we reuse the NURBS Toolbox described in [108]. We construct a NURBS surface as
shown in Fig.20a. The corresponding Matlab code is given in Listing 8. Using a uniform h-refinement (Listing 9) that
divides a knot span into two, we obtain the mesh given in Fig.20b. If one needs to use k−refinement (p−refinement
followed by h−refinement), then the code in Listing 10 can be used (see Fig.20c). Finally, Fig.20d gives the mesh which
is obtained by the process in which h−refinement is employed first and then p−refinement is performed (see Listing
11). After using the NURBS toolbox, the NURBS object is then converted to MIGFEM data structures using the
function convert2DNurbs located in folder nurbs-util.
Listing 8: Construct a NURBS surface using the NURBS toolbox.
a = 0 . 3 ; b = 0 . 5 ; % inner / outer radius
%% knots , c o n t r o l p o i n t s
uKnot = [ 0 0 0 1 1 1 ] ; %q u a d r a t i c b a s i s
vKnot = [ 0 0 1 1 ] ;
%l i n e a r b a s i s
% homogeneous c o o r d s ( x∗w, y∗w, z ∗w, w)
% 3 p t s i n u−d i r , 2 p t s i n v−d i r
controlPts
= zeros ( 4 , 3 , 2 ) ;
controlPts (1:2 ,1 ,1) = [ a ; 0 ] ;
controlPts (1:2 ,2 ,1) = [ a ; a ; ] ;
controlPts (1:2 ,3 ,1) = [ 0 ; a ] ;
controlPts (1:2 ,1 ,2) = [ b ; 0 ] ;
controlPts (1:2 ,2 ,2) = [ b ; b ] ;
controlPts (1:2 ,3 ,2) = [ 0 ; b ] ;
controlPts ( 4 , : , : )
= 1;
c o n t r o l P t s ( 4 , 2 , 1 ) = 1/ sqrt ( 2 ) ;
c o n t r o l P t s ( 4 , 2 , 2 ) = 1/ sqrt ( 2 ) ;
% homogenous c o o r d i n a t e s ( x∗w, y∗w, z ∗w)
controlPts (1:2 ,2 ,1) = controlPts (1:2 ,2 ,1)∗ fac ;
controlPts (1:2 ,2 ,2) = controlPts (1:2 ,2 ,2)∗ fac ;
57
(a) Initial NURBS mesh
(b) After h-refinement
(c) After k-refinement
(d) After hp-refinement
Figure 20: Illustration of the utilization of the NURBS toolbox in building NURBS object (a). From the initial mesh,
different meshes can be obtained using either h−refinement, p−refinement or combination thereof. As can be seen from
(c) and (d), k−refinement (c) is more efficient than hp−refinement (d). The function plotMesh.m in folder meshing
is used to plot NURBS mesh and control polygon.
58
%% b u i l d NURBS o b j e c t
s o l i d = nrbmak ( c o n t r o l P t s , { uKnot vKnot } ) ;
Listing 9: h-refinement using the NURBS toolbox.
1
refineLevel = 2;
2
f o r i =1: r e f i n e L e v e l
3
uKnotVectorU = u n i q u e ( uKnot ) ;
4
uKnotVectorV = u n i q u e ( vKnot ) ;
5
% new k n o t s a l o n g two d i r e c t i o n s
6
newKnotsX =uKnotVectorU ( 1 : end−1)+0.5∗ d i f f ( uKnotVectorU ) ;
7
newKnotsY =uKnotVectorV ( 1 : end−1)+0.5∗ d i f f ( uKnotVectorV ) ;
8
newKnots
= {newKnotsX newKnotsY } ;
10
solid
= n r b k n t i n s ( s o l i d , newKnots ) ;
11
uKnot
= cell2mat ( s o l i d . knots ( 1 ) ) ;
12
vKnot
= cell2mat ( s o l i d . knots ( 2 ) ) ;
9
13
end
14
convert2DNurbs % c o n v e r t t o t h e MIGFEM f o r m a t
15
plotMesh ( c o n t r o l P t s , w e i g h t s , uKnot , vKnot , p , q , r e s , ’r - ’ , ’ try . eps ’ ) ;
Listing 10: k-refinement using the NURBS toolbox.
1
% code from L i s t i n g 5
2
% o r d e r e l e v a t i o n , p=p+2 , q=q+1
3
4
5
s o l i d = nrbdegelev ( solid , [ 2
1]);
% then , knot i n s e r t i o n u s i n g t h e code from L i s t i n g 6
convert2DNurbs % c o n v e r t t o t h e MIGFEM f o r m a t
Listing 11: h-refinement followed by p-refinement using the NURBS toolbox.
1
% code from L i s t i n g 5
2
% then , knot i n s e r t i o n u s i n g t h e code from L i s t i n g 6
3
% o r d e r e l e v a t i o n , p=p+2 , q=q+1
4
s o l i d = nrbdegelev ( solid , [ 2
1]);
5
convert2DNurbs % c o n v e r t t o t h e MIGFEM f o r m a t
59
Listing 12: A typical input file for MIGFEM.
1
% code from L i s t i n g 5
2
% code from L i s t i n g 6 ( i f o n l y h−r e f i n e m e n t )
3
convert2DNurbs%c o n v e r t t o t h e MIGFEM f o r m a t ( i f not done y e t )
4
generateIGA2DMesh % b u i l d t h e mesh ( e l e m e n t c o n n e c t i v i t y )
6.6. Input file for MIGFEM
A typical input file is given in Listing 12. This file replaces the standard FE mesh file. Note that, for backward
compatibility with older versions of the MIGFEM code, some input files do not use the NURBS toolbox to create the
NURBS object. For those input files, it is however impossible to perform order elevation and hence k−refinement. It is,
therefore, recommended to create the NURBS objects using the NURBS toolbox [108] which, besides the aforementioned
refinement functionalities, also supports many useful operations such as extrusion, rotation etc.
Incorporating NURBS into an existing FE code cannot be considered a trivial task due to the presence of various
spaces (index, parameter) that are not present in conventional FE codes, Bézier extraction [130, 52] provides a FE
data structure that allows for a straightforward implementation of NURBS/T-splines into any FE codes. Appendix D
briefly presents this concept and its implementation in MIGFEM.
6.7. XFEM implementation
There are certainly different ways to implement XIGA. We present a way that reuse most of the tools present in an
existing XFEM code. We use the approximate Q4 mesh used for visualization as discussed in Section 6.4 for selection
of enriched control points. The level set values defining the crack at the vertices of this mesh are then computed.
Based on these level sets, elements cut by the crack and elements containing the crack tip can be determined [131]. For
example, element ABCD in Fig. 17 is cut by the crack. In a FEM context, its four nodes are then enriched using the
Heaviside function. In an isogeometric framework, however, the control points associated to this element are enriched.
Listing 13 details the Matlab R implementation of this process. Note that this listing is taken directly from our XFEM
code [119] with only one small modification and thus the proposed technique is considered simpler than that adopted
in [72]. We emphasize that the crack geometry is defined in the physical space to keep the usual XFEM notation.
Listing 13: Selection of enriched control points/nodes.
60
e n r i c h n o d e = zeros ( noCtrPts , 1 ) ;
c oun t1 = 0 ;
c oun t2 = 0 ;
f o r i e l = 1 : numelem
sctr
= elementV ( i e l , : ) ;
sctrIGA = e l e m e n t ( i e l , : ) ;
phi
= l s ( s c t r , 1 ) ; % normal l e v e l s e t
psi
= l s ( sctr , 2 ) ; % tangent l e v e l s e t
i f ( max( p h i ) ∗min( p h i ) < 0 )
i f max( p s i ) < 0
c ou nt 1 = co un t1 + 1 ; % one s p l i t e l e m e n t
s p l i t e l e m ( c ou nt 1 ) = i e l ;
e n r i c h n o d e ( sctrIGA )
= 1;
e l s e i f max( p s i ) ∗min( p s i ) < 0
c ou nt 2 = co un t2 + 1 ; % one t i p e l e m e n t
t i p e l e m ( c oun t 2 ) = i e l ;
e n r i c h n o d e ( sctrIGA )
= 2;
end
end
end
s p l i t n o d e s = find ( e n r i c h n o d e == 1 ) ;
tip nodes
= find ( e n r i c h n o d e == 2 ) ;
Remark 6.1. Note that for simple geometries as those tackled in this paper, the use of level sets is certainly not
necessary. More generally, describing open surfaces (cracks) with level sets requires two level sets functions that, for
crack growth simulations, must be reinitialised for stability (this decreases accuracy) and reorthogonalised every few
time steps. This is particularly cumbersome and probably explains along with the difficulties in dealing intersecting
and branching cracks the recent trend of research efforts in the area of phase field models of fracture, see e.g., [132, 74],
and the thick level set method [133].
61
√
Crack visualization We have H(x+ ) − H(x− ) = 2 and B1 (x+ ) − B1 (x− ) = 2 r, therefore the displacement jump
at a point x on the crack face is given by
[[u]](x) = 2
X
J∈Sc
√ X
RK (x)b1K .
RJ (x)aJ + 2 r
(90)
K∈Sf
Note that the other branch functions Bα , α = 2, 3, 4 are continuous functions and thus do not contribute to the
displacement jump.
crack
(a)
(b)
Figure 21: Crack visualization in XIGA: (a) build a mesh that is compatible to the crack by introducing double nodes
along the crack (square nodes) and (b) assign the displacement jumps to these new nodes. Note that the square nodes
in the tip element are used only for compatibility purposes. Exact mode I displacements are imposed on the bottom,
right and top edge using the Lagrange multiplier method while Neumann BCs from the exact stress field are enforced
on the left edge. We refer to [119] for a detailed description of this standard problem.
Fig. 21 illustrates the idea for crack visualization with the script post-processing/crackedMeshNURBS.m
providing implementation details. The contour plots of the displacement and stress field of a mode I cracked sample
are given in Fig. 22. Note that the stresses at points on the crack are simply set to zero (traction-free cracks) and the
stresses of the new nodes of the tip element are interpolated from the values of the four nodes of this Q4 element.
Remark 6.2. As is the case for XFEM, integration over elements cut by the cracks usually requires subdivision of
the elements into integration subcells. We refer to [134] for a recent discussion on this issue. In addition to this
popular technique, a simple integration rule is also provided–elements crossed by the crack and tip-enriched elements
62
Figure 22: Contour plots on a cracked mesh: (a) vertical displacement and (b) normal stress in the vertical direction.
are numerically integrated using a regular Gauss-Legendre quadrature with a large number of Gauss points as done in
[135].
7. Structural mechanics
Thanks to the high order continuity provided by NURBS/T-splines, the implementation of rotation-free thin beam/plate/shell elements becomes direct and simple. In this section, we are going to present the implementation of a rotationfree IGA Kirchhoff plate formulation (rotation free shell elements can be found in folder structural-mechanics). The
plate geometry and the deflection are both approximated by NURBS. At control points there is only one unknown-the
deflection or transverse displacement. For simplicity only isotropic elastic plates are considered. We refer to [136] for
a treatment of plate theories.
The element stiffness matrix is defined as
Ke =
Z
Ωe
BT
e DBe dΩ,
63
(91)
where the constitutive matrix D reads
1
Eh3
ν
D=
12(1 − ν 2 )
0
ν
0
1
0
0
0.5(1 − ν)
,
(92)
where E, ν are the Young’s modulus and Poisson’s ratio, respectively; h denotes the plate thickness and the element
displacement-curvature matrix Be that contains the second derivatives of the shape functions is given by
R1,xx
Be =
R1,yy
2R1,xy
R2,xx
R2,yy
2R2,xy
· · · Rn,xx
· · · Rn,yy
,
· · · 2Rn,xy
(93)
where n denotes the number of basis functions of element e and RA,xx ≡ d2 RA /dx2 .
7.1. Boundary conditions
For clamped BCs one needs to fix the rotations. The nodal unknowns are, however, only the transverse displacements
w. To fix the rotation of a boundary, we simple fix two row of control points at the boundary [20] because these control
points define the tangent of the surface at the boundary (cf. Fig. 3), see Fig. 23.
Figure 23: Enforcing BCs for a fully clamped plate: simply fixing the deflections of two rows of control points around
the clamped boundary. Note that the set of CPs next to the boundary CPs are not artificially added to impose the
rotations. They are simply the CPs defining the geometry of the plate.
64
7.2. Symmetry boundary conditions
Fig. 24 illustrates the use of symmetry boundary conditions when only 1/4 of the plate is modelled. Along the
symmetry lines, the rotation should be zero which can be enforced by constraining the deflection (w) of two rows of
control points along these lines together. These contraints can be implemented using a simple penalty technique as
shown in Listing 14.
cou
sym
fix
sym
fix
Figure 24: A fully clamped square plate: 1/4 model is analyzied using appropriate symmetry BCs. Along the symmetry
lines, the rotation is fixed which can be achieved by enforcing the deflection of two rows of control points that define
the tangent of the plate to have the same value.
Listing 14: Enforcing symmetry BCs.
w
= 1 e7 ;
p e n a l t y S t i f f n e s s = w∗ [ 1 −1;−1 1 ] ;
f o r i =1: length ( topNodes )
sctr
= [ topNodes ( i ) nextToTopNodes ( i ) ] ;
K( s c t r , s c t r ) = K( s c t r , s c t r ) + p e n a l t y S t i f f n e s s ;
end
% t h e same f o r t h e l e f t two rows
For completeness, we give the implementation of the rotation free Kirchhoff plate elements in Listing 15. Note that
we have skipped codes that are common with IGA code for 2D continua.
Listing 15: Computation of stiffness matrix for the rotation free Kirchhoff plate.
1
f o r gp =1: s i z e (W, 1 )
65
pt
2
= Q( gp , : ) ; wt
% shape f u n c t i o n s ,
3
= W( gp ) ;
f i r s t and s e c o n d d e r i v a t i v e s w . r . t n a t u r a l c o o r d s
4
[R dRdxi dRdeta dR2dxi dR2det dR2dxe]=NURBS2DBasis2ndDers ( [ Xi ; Eta ] , p , q , uKnot , vKnot , w e i g h t s ’ ) ;
5
jacob
6
j a c o b 2 = [ dR2dxi ; dR2det ; dR2dxe ] ∗ p t s ; % 3 x2 m a t r i x
7
J1
8
dxdxi = j a c o b ( 1 , 1 ) ; dydxi = j a c o b ( 1 , 2 ) ;
9
dxdet = j a c o b ( 2 , 1 ) ; dydet = j a c o b ( 2 , 2 ) ;
j33
= [ dRdxi ; dRdeta ]
∗ p t s ; % 2 x2 m a t r i x
= det ( j a c o b ) ;
= [ dxdxi ˆ2
dydxi ˆ2
2∗ dxdxi ∗ dydxi ;
11
dxdet ˆ2
dydet ˆ2
2∗ dxdet ∗ dydet ;
12
dxdxi ∗ dxdet dydxi ∗ dydet dxdxi ∗ dydet+dxdet ∗ dydxi ] ;
10
% J a c o b i a n i n v e r s e and s p a t i a l 1 s t and 2nd d e r i v a t i v e s
13
14
invJacob
= inv ( j a c o b ) ;
15
dRdx
= i n v J a c o b ∗ [ dRdxi ; dRdeta ] ;
16
dR2dx
= inv ( j 3 3 ) ∗ ( [ dR2dxi ; dR2det ; dR2dxe]− j a c o b 2 ∗dRdx ) ;
17
% B matrix
18
B
= dR2dx ;
19
B( 3 , : )
= B( 3 , : ) ∗ 2 ;
20
% compute e l e m e n t a r y s t i f f n e s s matrix , s c t r=e l e m e n t ( e , : )
21
K( s c t r , s c t r ) = K( s c t r , s c t r ) + B’ ∗ C ∗ B ∗ J1 ∗ J2 ∗ wt ;
22
end
8. Verification examples
In this section, numerical examples in linear elasticity and linear elastic fracture mechanics in 2D and 3D are
presented with the purpose to serve as a set of verification examples for MIGFEM. They include an infinite plate
with a circular hole under constant in-plane tension, the pinched cylinder, an edge cracked plate in tension, a threedimensional mode I fracture problem and a large deformation thin shell problem. Unless otherwise stated, standard
direct imposition of Dirichlet BCs is used. Units are standard International System (SI) units.
66
8.1. Two and three dimensional solid mechanics
8.1.1. Infinite plate with a circular hole
The problem considered is that of an infinite plate with a circular hole in the centre under constant in-plane tension
at infinity as shown in Fig. 25 where, due to symmetry, only a quarter of the plate is modeled. The plate dimension
is taken to be L × L and the circular hole has a radius R. The exact stress field in the plate is given by
σxx (r, θ) = 1 −
R2
r2
3
3 R4
cos 2θ + cos 4θ +
cos 4θ
2
2 r4
(94a)
R2
σyy (r, θ) = − 2
r
3 R4
1
cos 2θ − cos 4θ −
cos 4θ
2
2 r4
(94b)
R2
σxy (r, θ) = − 2
r
3 R4
1
sin 2θ + sin 4θ +
sin 4θ,
2
2 r4
(94c)
where r, θ are the usual polar coordinates centered at the center of the hole.
The material properties are specified as E = 103 , Poisson’s ratio ν = 0.3 and the geometry is such that L = 4,
R = 1. A plane stress condition is assumed. The problem is solved with quadratic NURBS meshes such as those shown
in Fig. 26. The control points and weights for the coarsest mesh can be found in [3] or file plateHoleCkData.m.
Fig. 27, generated in Paraview, illustrates the contour plot of numerical σxx . Note that the stress concentration at
point (R, 3π/2) is well captured and a smooth stress field is obtained throughout.
Remark 8.1. Using the visualization technique described in Section 6.4 for this problem, a note should be made on the
evaluation of the stress field at the top left corner where there are two coincident control points. This causes a singular
Jacobian matrix. Therefore at this corner, the stresses at a point slightly shifted from the original position are used.
8.1.2. Pinched cylinder
In order to demonstrate the performance of the 3D IGA implementation, we consider the pinched cylinder problem
as shown in Fig. 28. Note that we discretize the shell with solid NURBS elements. Due to symmetry, only 1/8 of
the model is analyzed. A tri-quadratic NURBS mesh (p = q = r = 2) was used for the computation. Details can be
67
y
ex
D
E
eq
v1
nu
eq
v2
R
eq
v3
l
eq
v4
ex
l
sy
e
R
C
theta
A
x
B
sy
r
Figure 25: Infinite plate with a circular hole under constant in-plane tension: quarter model. Boundary conditions
T
T
include: uy = 0 (AB), ux = 0 (CD), t = (−σxx , −σxy ) (AE), t = (σxy , σyy ) (ED).
Figure 26: Plate with a hole: coarse mesh of 2 bi-quadratic elements (left) and refined mesh (right).
68
Figure 27: Plate with a hole: distribution of numerical σxx obtained with a 32 × 16 quadratic mesh having 4488 dofs.
found in the file igaPinchedCylinder.m. Fig. 29 shows the mesh and the contour plot of the displacement in the
point load direction. Post-processing is done in Paraview and we refer to Section 6.4 for details. We recognise that the
problem under consideration is a shell like structure that would be more accurately modelled using appropriate shell
elements, but the example is merely intended to illustrate the ability of the method to analyse 3D geometries.
Figure 28: Pinched cylinder. Problem description and data [137].
69
Figure 29: Pinched cylinder: (a) mesh of one octant of the cylinder and (b) contour plot of the displacement in the
direction of the point load.
8.2. Two and three dimensional fracture mechanics
8.2.1. Edge cracked plate in tension
A plate of dimension b × 2h is loaded by a tensile stress σ = 1.0 along the top edge and bottom edge as shown in
Fig. 30. In the computation, the displacement along the y-axis is fixed at the bottom edge and the bottom left corner
is fixed in both x and y directions. The material parameters are E = 103 and ν = 0.3. A plane strain condition is
assumed. The reference mode I stress intensity factor (SIF) for this problem is given in [138] and is calculated as
KI = F
a √
σ πa,
b
(95)
where a is the crack length, b is the plate width, and F (a/b) is an empirical function. For a/b ≤ 0.6, function F is
given by
F
a
b
= 1.12 − 0.23
a
b
+ 10.55
a 2
b
− 21.72
a 3
b
+ 30.39
a 4
b
.
(96)
In the present implementation, this problem is solved using both XFEM and an extended isogeometric formulation.
The SIF is computed using an interaction integral. We refer to [70] for details.
We first verify the implementation of the XIGA code by comparing the XIGA result with the XFEM result for the
case of a = 0.45, b = 1 and h = 1. The XFEM and XIGA meshes are given in Fig. 31. Both meshes have the same
uniform distribution of nodes/control points. For the XFEM mesh, bilinear Q4 elements are used. For the XIGA mesh,
70
h
si
h
a
b
si
Figure 30: Edge cracked plate in tension: geometry and loading.
cubic (p = q = 3) B-spline basis functions are adopted. Fig. 32 shows the contour plots of the vertical displacement
obtained with XFEM and XIGA.
We now consider the computation of the mode I SIF for a crack of length a = 0.3. The reference SIF for this
problem is KIref = 1.6118. Both a linear and a cubic B-spline basis are used for three different meshes. The results are
given in Table 1. It should be emphasized that in the computation of the interaction integral, we use bilinear Lagrange
shape functions i.e., shape functions of Q4 elements to compute the derivatives of the weight function. This guarantees
that the weight function takes a value of unity on an open set containing the crack tip and vanishes on an outer contour
as shown in Fig. 33. B-splines functions are not interpolatory and therefore cannot be used to approximate the weight
functions.
mesh
disp. dofs
KI (linear)
Error (%)
KI (cubic)
Error (%)
9 × 18
18 × 36
36 × 72
324
1296
5184
1.4997
1.5823
1.5968
6.96
1.83
0.93
1.5560
1.6150
1.6117
3.46
0.20
0.01
Table 1: Edge cracked plate: SIFs results. The reference SIF is KIref = 1.6118. Note that linear NURBS are equivalent
to the conventional bilinear finite elements.
71
(a) XFEM mesh (Q4)
(b) XIGA
Bsplines)
mesh
(cubic
(c) XIGA mesh, enriched nodes
Figure 31: Edge cracked plate: XFEM and XIGA meshes. Both have the same number of displacement dofs of 1296.
The thick line denotes the crack. Square nodes denote tip enriched nodes whereas star nodes represent Heaviside
enriched nodes.
(a) XFEM
(b) XIGA
Figure 32: Edge cracked plate: uy contour plots on deformed configuration (enlargement factor of 30 used).
72
Figure 33: Distribution of weight function used in the computation of the interaction integral. Four-noded quadrilateral
elements with bilinear Lagrange shape functions are used to interpolate the weight function.
8.2.2. Three-dimensional mode I fracture problem
This example aims to show the capability of MIGFEM for solving three-dimensional (3D) fracture problems. For
3D cracks, the polar coordinates in the branch functions given in Eq. (83) are defined in terms of the level sets as [139]
r=
p
ϕ(ξ, η, ζ)2 + ψ(ξ, η, ζ)2 ,
θ = atan
ϕ(ξ, η, ζ)
ψ(ξ, η, ζ)
,
(97)
where the level set field Φ = (ϕ, ψ) is interpolated as
Φ(ξ, η, ζ) =
X
RI (ξ, η, ζ)ΦI .
(98)
I
We refer to [139] for details concerning the derivatives of the branch functions with respect to the parametric coordinates
(ξ, η, ζ).
73
The mode I 3D fracture problem we are solving is given in Fig. 34. The exact displacement field is given by
ux (r, θ) =
θ
θ
2(1 + υ) KI √
√
r cos
2 − 2υ − cos2
2
2
2π E
(99)
uy (r, θ) = 0
uz (r, θ) =
2(1 + υ) KI √
θ
θ
√
r sin
2 − 2υ − cos2
,
2
2
2π E
√
where KI = σ πa is the stress intensity factor, υ is Poisson’s ratio and E is Young’s modulus. In our example, a = 100
2
mm; E = 107 N/mm , υ = 0.3, σ = 104 N/mm2 . On the bottom, right and top surfaces, essential BCs taken from Eq.
(99) are imposed using the penalty method. A penalty parameter of 1e10 was used. We note that this problem can
be more effectively solved with two-dimensional elements. This example however aims at presenting how 3D extended
IGA can be implemented. Furthermore, it also illustrates how Dirichlet BCs are enforced on surfaces rather than the
usual case of line boundaries. To this end, a two-dimensional NURBS mesh for a given surface is generated from the
set of control points that define this surface (see the file surfaceMesh.m).
10
z
y
x
10
Figure 34: Three-dimensional mode I fracture problem: infinite plate with a center planar crack. The plate thickness
is 2, the crack length is 5 and the crack width is 2. The crack is located in the mid-plane of the plate.
The problem is first solved using a linear B-spline basis. A mesh of 9 × 9 × 1 elements is used. The mesh, enriched
74
nodes and comparison of the numerical deformed configuration against the exact profile are given in Fig. 35. Next, a
mesh of 7 × 7 × 2 elements is used where, in the through-thickness direction, there are two linear elements (q = 1) and
for the two other directions, a quadratic basis (p = r = 2) is used. The result is given in Fig. 36 and we note a good
qualitative agreement between Figs. 35 and 36. We did not perform a SIF computation for this problem as 3D SIF
computation is not yet implemented at present.
Figure 35: Three-dimensional mode I fracture problem: mesh of linear B-spline elements and enriched control points
(left); numerical deformed configuration (magnification factor of 40) superimposed on the exact deformed configuration
(right).
8.3. Structural mechanics
We consider the pull out of an open-ended cylindrical shell, see Fig. 37, as one of the most common benchmarks
for geometrically nonlinear thin shell problems [140, 141]. Due to symmetry only 1/8 of the model is studied. The
1/8 model can be exactly described using one single quadratic-linear NURBS surface as shown in Fig.38. Refined
meshes are then obtained from this initial mesh by using the k-refinement. The analysis was performed using a mesh of
8 × 8 bi-quartic NURBS elements (144 control points). The enforcement of symmetry BCs is achieved by constraining
two row of control points as shown in the same figure. These constraints–the so-called multipoint linear homogeneous
constraints are handled using the master-slave method [142] for the penalty method as described in Section 7.2 could
endanger the convergence of the Newton-Raphson solver. Note that for this particular nonlinear problem, of which the
M file is igaGNLFreeEndCylinderShell.m, the analysis was performed with a C++ implementation [143].
75
Figure 36: Three-dimensional mode I fracture problem: mesh of quadratic B-spline elements and enriched control points
(left); numerical deformed configuration (magnification factor of 40) superimposed on the exact deformed configuration
(right).
We adopt the Kirchhoff-Love thin shell model, see e.g., [144, 141] for details that involves only displacement dofs.
The maximum applied load 0.25Pmax is applied in 80 equal increments and for each increment the full Newton-Raphson
method is used to solve for the displacements. Fig. 39 shows the deformed configuration of the shell and in Fig.40 the
load versus the z-displacement at point A is plotted together with the result reported in [140].
Figure 37: The open-end cylindrical shell subjected to radial pulling forces: problem setup [140].
76
ux
sym2
sym1
f
uxxx
sym3
uxx
ui
ui1
free
Figure 38: Open ended cylindrical shell: 1/8 of the model with one bi-quadratic NURBS surface. The control points
and weights are given in file freeEndsGNLCylinderShellData.m. For symmetry edges, in order to satisfy the
symmetry condition i.e., zero rotation, we constraint two rows of control points in the sense that the displacements of
the control points right next to the control points locating on the symmetry edges (uI1 ) matches uI .
Figure 39: The open-ended cylindrical shell subjected to radial pulling forces: deformed shape without scaling.
9. Conclusion
We presented a Matlab R implementation for one, two and three-dimensional isogeometric finite element analysis for
solid and structural mechanics. This paper is addressed to students or researchers who wish to learn the concepts of IGA
in a clear and concise manner and is especially suited to those with solid mechanics applications in mind. NURBS are
77
40
IGA
reference
35
P x 1000
30
25
20
15
10
5
0
0
0.5
1
1.5
2
2.5
3
displacement u
Figure 40: The open-ended cylindrical shell subjected to radial pulling forces: load-displacement responses.
used throughout, where the underlying construction of the basis functions is detailed along with associated refinement
algorithms essential for numerical analysis. Differences with conventional FE implementations are made clear with
the use of Matlab R source code to illustrate isogeometric FE concepts explicitly. In addition, the implementation of
an isogeometric XFEM formulation for both two-dimensional and three-dimensional problems is described allowing
for simple linear elastic fracture analysis to be performed directly from CAD data. The benefits of isogeometric
analysis are made evident. Although not presented, the code supports multi-patch analysis in which compatibility
between connecting patches is not required [145, 146]. Geometrical nonlinearities for solid elements under the Total
Lagrange framework are provided as well. In addition, PUM enrichment for holes and inclusions is provided along
with implementation of the least squares method for imposing essential boundary conditions. Mass matrices are
implemented for almost every elements and popular time integration schemes such as Newmark and central difference
explicit are available so that transient analysis on CAD objects can be performed. The code is available for download
from https://sourceforge.net/projects/cmcodes/ for Linux and Mac OS machines. A Windows version can be
found at http://www.softpedia.com/get/Science-CAD/igafem.shtml.
The preliminary concepts and implementation details of isogeometric analysis have been described, but many
challenges remain in the field. These include the creation of suitable volume discretisations from given CAD boundary
representations, efficient integration schemes and suitable error estimators.
78
Acknowledgements
The authors would like to acknowledge the partial financial support of the Framework Programme 7 Initial Training
Network Funding under grant number 289361 “Integrating Numerical Simulation and Geometric Design Technology”.
Stéphane Bordas also thanks partial funding for his time provided by 1) the EPSRC under grant EP/G042705/1
Increased Reliability for Industrially Relevant Automatic Crack Growth Simulation with the eXtended Finite Element
Method and 2) the European Research Council Starting Independent Research Grant (ERC Stg grant agreement No.
279578) entitled “Towards real time multiscale simulation of cutting in non-linear materials with applications to surgical
simulation and computer guided surgery”. The first author would like to express his gratitude towards Professor L.J.
Sluys at Delft University of Technology, The Netherlands for his support to VPN during the PhD period and the
Framework Programme 7 Initial Training Network Funding. The authors acknowledge the comments of Dr. Robert
Simpson at Cardiff University that has improved the paper.
Appendix A. Knot vector conventions
The definition of an open knot vector as Ξ = {ξ1 , ξ2 , . . . , ξn+p+1 } where the first and last knot values are repeated
p + 1 times is sometimes simplified to remove redundant information. It is found that if a knot vector is open, the first
and last knot values play no role in the definition of the curve and can be removed entirely. Knot vectors are then
defined as Ξ = {ξ1 , ξ2 , . . . , ξn+p−1 } where the first and last knot vectors are repeated p times. This notation is used
frequently through the CAGD community.
Furthermore, B-spline algorithms can be written entirely in terms of knot interval vectors [147], defined as
∆Ξ = {∆ξ1 , ∆ξ2 , . . . , ∆ξn+p−2 }
where
∆ξi = ξi+1 − ξi
(A.1)
using the previous simplified knot vector notation. T-splines are based entirely on knot interval vector notation.
Appendix B. C k and Gk continuity
Two types of continuity are commonplace in the CAGD community, referred to as C and G continuity. From a
numerical analysis perspective, C continuity refers to the traditional definition of continuity in which, given a parametric
79
function f : R → R2 which defines a 2D curve through a parameter t, the parameterisation is said to be C k continuous
at point a if dn f /dtn is continuous (single-valued) for n = 0, 1, . . . , k at t = a. Intuitively, for a curve to be C 0 , it must
contain no discontinuities with t single-valued throughout its domain. Letting df /dt represent the ‘speed’ of the curve,
C 1 continuity implies that the tangent vector is equal in magnitude and direction while also assuming C 0 continuity.
C 2 continuity assumes both C 0 and C 1 continuity and requires that the ‘acceleration’ of the curve is equal in magnitude
and direction.
G continuity (sometimes referred to as ‘visual’ continuity) is independent of the parameter t. For a curve to be G0
continuous at a point it must contain no discontinuities, but the parameter value may have multiple values. For G1
continuity, the tangent of the curve must be equal in its direction, but its magnitude may differ. And G2 continuity
implies that the direction of acceleration is continuous, but its magnitude may differ.
If desired, a Gk curve can be made C k -continuous through an appropriate ‘reparameterisation’.
Appendix C. Homogenous and non-homogenous coordinates
It is found that if rational basis functions are used, then it is often more appropriate to work in terms of homogeneous
coordinates, since they can be applied straightforwardly to non-rational basis function algorithms. For example, in the
case of NURBS, the use of homogenous coordinates allows B-spline algorithms to be used directly.
Given a physical coordinate PA = (xA , yA , zA ) with weighting wA , the corresponding four-dimensional homogeneous
coordinate P̃A is written as
P̃A = (xA wA , yA wA , zA wA , wA )
= (x̃A , ỹA , z̃A , wA ).
(C.1)
(C.2)
To convert back to non-homogenous coordinates, it is a simple case of dividing x̃A , ỹA , z̃A through by wA .
Appendix D. Bézier extraction
Bézier extraction [130, 52] provides a tool to facilitate the incorporation of NURBS and even T-splines in to any
FE codes. The idea is based on the fact that any B-splines basis can be written as a linear combination of Bernstein
80
polynomials. Written for element e, the shape functions read
Ne (ξ) = Ce B(ξ̃)
(D.1)
where Ce denotes the elemental Bézier extraction operator and B are the Bernstein polynomials which are defined
on the parent element Ω̃. Index space, knot vectors are all embedded in the Bézier extrators which are computed in
a pre-processing step. Therefore, existing FE solvers can use NURBS/T-splines straightforwardly. We refer to folder
bezier-extraction for examples.
References
[1] L. A. Piegl and W. Tiller. The NURBS Book. Springer, 1996.
[2] D. F. Rogers. An Introduction to NURBS with Historical Perspective. Academic Press, 2001.
[3] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact
geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194(39-41):4135–
4195, 2005.
[4] J. A. Cottrell, T. J.R. Hughes, and Y. Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA.
Wiley, 2009.
[5] P. Kagan, A. Fischer, and P. Z. Bar-Yoseph. New B-Spline Finite Element approach for geometrical design and
mechanical analysis. International Journal for Numerical Methods in Engineering, 41(3):435–458, 1998.
[6] P. Kagan and A. Fischer. Integrated mechanically based CAE system using B-Spline finite elements. ComputerAided Design, 32(8–9):539 – 552, 2000.
[7] F. Cirak, M. Ortiz, and P. Schröder. Subdivision surfaces: a new paradigm for thin-shell finite-element analysis.
International Journal for Numerical Methods in Engineering, 47(12):2039–2072, 2000.
[8] İ. Temizer, P. Wriggers, and T.J.R. Hughes. Contact treatment in isogeometric analysis with NURBS. Computer
Methods in Applied Mechanics and Engineering, 200(9-12):1100–1112, 2011.
81
[9] L. Jia. Isogeometric contact analysis: Geometric basis and formulation for frictionless contact. Computer Methods
in Applied Mechanics and Engineering, 200(5-8):726–741, 2011.
[10] İ. Temizer, P. Wriggers, and T.J.R. Hughes. Three-Dimensional Mortar-Based frictional contact treatment in
isogeometric analysis with NURBS. Computer Methods in Applied Mechanics and Engineering, 209–212:115–128,
2012.
[11] L. De Lorenzis, İ. Temizer, P. Wriggers, and G. Zavarise. A large deformation frictional contact formulation using
NURBS-bases isogeometric analysis. International Journal for Numerical Methods in Engineering, 87(13):1278–
1300, 2011.
[12] M.E. Matzen, T. Cichosz, and M. Bischoff. A point to segment contact formulation for isogeometric, NURBS
based finite elements. Computer Methods in Applied Mechanics and Engineering, 255:27 – 39, 2013.
[13] W. A. Wall, M. A. Frenzel, and C. Cyron. Isogeometric structural shape optimization. Computer Methods in
Applied Mechanics and Engineering, 197(33-40):2976–2988, 2008.
[14] N. D. Manh, A. Evgrafov, A. R. Gersborg, and J. Gravesen. Isogeometric shape optimization of vibrating
membranes. Computer Methods in Applied Mechanics and Engineering, 200(13-16):1343–1353, 2011.
[15] X. Qian and O. Sigmund. Isogeometric shape optimization of photonic crystals via Coons patches. Computer
Methods in Applied Mechanics and Engineering, 200(25-28):2237–2255, 2011.
[16] X. Qian. Full analytical sensitivities in NURBS based isogeometric shape optimization. Computer Methods in
Applied Mechanics and Engineering, 199(29-32):2059–2071, 2010.
[17] R.N. Simpson, S.P.A. Bordas, J. Trevelyan, and T. Rabczuk. A two-dimensional isogeometric boundary element
method for elastostatic analysis. Computer Methods in Applied Mechanics and Engineering, 209–212:87–100,
2012.
[18] M.A. Scott, R.N. Simpson, J.A. Evans, S. Lipton, S.P.A. Bordas, T.J.R. Hughes, and T.W. Sederberg. Isogeometric boundary element analysis using unstructured T-splines. Computer Methods in Applied Mechanics and
Engineering, 254:197 – 221, 2013.
82
[19] D.J. Benson, Y. Bazilevs, M.C. Hsu, and T.J.R. Hughes. Isogeometric shell analysis: The Reissner–Mindlin shell.
Computer Methods in Applied Mechanics and Engineering, 199(5-8):276–289, 2010.
[20] J. Kiendl, K.-U. Bletzinger, J. Linhard, and R. Wüchner. Isogeometric shell analysis with Kirchhoff-Love elements.
Computer Methods in Applied Mechanics and Engineering, 198(49-52):3902–3914, 2009.
[21] D.J. Benson, Y. Bazilevs, M.-C. Hsu, and T.J.R. Hughes. A large deformation, rotation-free, isogeometric shell.
Computer Methods in Applied Mechanics and Engineering, 200(13-16):1367–1378, 2011.
[22] L. Beirão da Veiga, A. Buffa, C. Lovadina, M. Martinelli, and G. Sangalli. An isogeometric method for the
Reissner-Mindlin plate bending problem. Computer Methods in Applied Mechanics and Engineering, 209–212:45–
53, 2012.
[23] T. K. Uhm and S. K. Youn. T-spline finite element method for the analysis of shell structures. International
Journal for Numerical Methods in Engineering, 80(4):507–536, 2009.
[24] R. Echter, B. Oesterle, and M. Bischoff. A hierarchic family of isogeometric shell finite elements. Computer
Methods in Applied Mechanics and Engineering, 254:170 – 180, 2013.
[25] D.J. Benson, S. Hartmann, Y. Bazilevs, M.-C. Hsu, and T.J.R. Hughes. Blended isogeometric shells. Computer
Methods in Applied Mechanics and Engineering, 255:133 – 146, 2013.
[26] J. Kiendl, Y. Bazilevs, M.-C. Hsu, R. Wüchner, and K.-U. Bletzinger. The bending strip method for isogeometric
analysis of Kirchhoff-Love shell structures comprised of multiple patches. Computer Methods in Applied Mechanics
and Engineering, 199(37-40):2403–2416, 2010.
[27] Lu J. Circular element: Isogeometric elements of smooth boundary. Computer Methods in Applied Mechanics
and Engineering, 198(30-32):2391–2402, 2009.
[28] J. Lu and X. Zhou. Cylindrical element: Isogeometric model of continuum rod. Computer Methods in Applied
Mechanics and Engineering, 200(1-4):233–241, 2011.
[29] H. Gomez, T.J.R. Hughes, X. Nogueira, and V. M. Calo. Isogeometric analysis of the isothermal Navier-StokesKorteweg equations. Computer Methods in Applied Mechanics and Engineering, 199(25-28):1828–1840, 2010.
83
[30] P. N. Nielsen, A. R. Gersborg, J. Gravesen, and N. L. Pedersen. Discretizations in isogeometric analysis of
Navier-Stokes flow. Computer Methods in Applied Mechanics and Engineering, 200(45-46):3242–3253, 2011.
[31] Y. Bazilevs and I. Akkerman. Large eddy simulation of turbulent Taylor-Couette flow using isogeometric analysis
and the residual-based variational multiscale method. Journal of Computational Physics, 229(9):3402–3414, 2010.
[32] Y. Bazilevs, V. M. Calo, T. J. R. Hughes, and Y. Zhang. Isogeometric fluid-structure interaction: theory,
algorithms, and computations. Computational Mechanics, 43:3–37, 2008.
[33] Y. Bazilevs, J.R. Gohean, T.J.R. Hughes, R.D. Moser, and Y. Zhang. Patient-specific isogeometric fluid-structure
interaction analysis of thoracic aortic blood flow due to implantation of the Jarvik 2000 left ventricular assist
device. Computer Methods in Applied Mechanics and Engineering, 198(45-46):3534–3550, 2009.
[34] H. Gómez, V. M. Calo, Y. Bazilevs, and T.J.R. Hughes. Isogeometric analysis of the Cahn-Hilliard phase-field
model. Computer Methods in Applied Mechanics and Engineering, 197(49-50):4333–4352, 2008.
[35] C. V. Verhoosel, M. A. Scott, T. J. R. Hughes, and R. de Borst. An isogeometric analysis approach to gradient
damage models. International Journal for Numerical Methods in Engineering, 86(1):115–134, 2011.
[36] P. Fischer, M. Klassen, J. Mergheim, P. Steinmann, and R. Müller. Isogeometric analysis of 2D gradient elasticity.
Computational Mechanics, 47:325–334, 2010.
[37] A. Masud and R. Kannan. B-splines and NURBS based finite element methods for Kohn-Sham equations.
Computer Methods in Applied Mechanics and Engineering, 241-244:112 – 127, 2012.
[38] J.A. Cottrell, A. Reali, Y. Bazilevs, and T.J.R. Hughes. Isogeometric analysis of structural vibrations. Computer
Methods in Applied Mechanics and Engineering, 195(41-43):5257–5296, 2006.
[39] T.J.R. Hughes, A. Reali, and G. Sangalli. Duality and unified analysis of discrete approximations in structural
dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS. Computer
Methods in Applied Mechanics and Engineering, 197(49–50):4104 – 4124, 2008.
[40] C. H. Thai, H. Nguyen-Xuan, N. Nguyen-Thanh, T-H. Le, T. Nguyen-Thoi, and T. Rabczuk. Static, free vibration,
84
and buckling analysis of laminated composite Reissner-Mindlin plates using NURBS-based isogeometric approach.
International Journal for Numerical Methods in Engineering, 91(6), 2012.
[41] D. Wang, W. Liu, and H. Zhang. Novel higher order mass matrices for isogeometric structural vibration analysis.
Computer Methods in Applied Mechanics and Engineering, pages –, 2013.
[42] J. A. Evans, Y. Bazilevs, I. Babuška, and T.J.R. Hughes. n-Widths, sup-infs, and optimality ratios for the
k-version of the isogeometric finite element method. Computer Methods in Applied Mechanics and Engineering,
198(21-26):1726–1741, 2009.
[43] A. Shaw and D. Roy. Nurbs-based parametric mesh-free methods. Computer Methods in Applied Mechanics and
Engineering, 197(17–18):1541 – 1567, 2008.
[44] Hyun-Jung Kim and Sung-Kie Youn. Spline-based meshfree method. International Journal for Numerical Methods
in Engineering, 92(9):802–834, 2012.
[45] D. Großmann, B. Jüttler, H. Schlusnus, J. Barner, and Anh-Vu Vuong. Isogeometric simulation of turbine blades
for aircraft engines. Computer Aided Geometric Design, 29(7):519 – 531, 2012.
[46] Thomas Elguedj, Julien Réthoré, and Aurélien Buteri. Isogeometric analysis for strain field measurements.
Computer Methods in Applied Mechanics and Engineering, 200(1-4):40–56, 2011.
[47] T. W. Sederberg, J. Zheng, A. Bakenov, and A. Nasri. T-splines and T-NURCCs. ACM Transactions on Graphics,
22:477–484, 2003.
[48] A. Buffa, D. Cho, and G. Sangalli. Linear independence of the T-spline blending functions associated with some
particular T-meshes. Computer Methods in Applied Mechanics and Engineering, 1437-1445:199, 2010.
[49] M.A. Scott, X. Li, T.W. Sederberg, and T.J.R. Hughes. Local refinement of analysis suitable T-splines. Computer
Methods in Applied Mechanics and Engineering, 213-216:206–222, 2012.
[50] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J.A. Evans, T.J.R. Hughes, S. Lipton, M.A. Scott, and T.W. Sederberg.
Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering, 199(5-8):229–
263, 2010.
85
[51] M. R. Dörfel, B. Jüttler, and B. Simeon. Adaptive isogeometric analysis by local h-refinement with T-splines.
Computer Methods in Applied Mechanics and Engineering, 199(5-8):264–275, 2010.
[52] M. A. Scott, M. J. Borden, C. V. Verhoosel, T. W. Sederberg, and T. J. R. Hughes. Isogeometric finite element data structures based on Bézier extraction of T-splines. International Journal for Numerical Methods in
Engineering, 88(2):126–156, 2011.
[53] H. Wang, Y. He, X. Li, X. Gu, and H. Qin. Polycube splines. Computer-Aided Design, 40:721–733, 2008.
[54] J. Deng, F. Chen, X. Li, C. Hu, W. Tong, Z. Yang, and Y. Feng. Polynomial splines over hierarchical T-meshes.
Graphical Models, 70:76–86, 2008.
[55] T. Dokken and V. Skytt. Locally refined splines. In Proceedings of IV European Conference On Computational
Mechanics. Solids, Structures and Coupled Problems in Engineering, Paris, France, 16-21 May 2010.
[56] N. Nguyen-Thanh, H. Nguyen-Xuan, S.P.A. Bordas, and T. Rabczuk. Isogeometric analysis using polynomial
splines over hierarchical T-meshes for two-dimensional elastic solids. Computer Methods in Applied Mechanics
and Engineering, 200(21-22):1892–1908, 2011.
[57] N. Nguyen-Thanh, J. Kiendl, H. Nguyen-Xuan, R. Wüchner, K.U. Bletzinger, Y. Bazilevs, and T. Rabczuk.
Rotation free isogeometric thin shell analysis using PHT-splines. Computer Methods in Applied Mechanics and
Engineering, 200(47-48):3410–3424, 2011.
[58] A.-V. Vuong, C. Giannelli, B. Jüttler, and B. Simeon. A hierarchical approach to adaptive local refinement in
isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 200(49-52):3554–3567, 2011.
[59] D. Schillinger, L. Dedé, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, and T. J. R. Hughes. An Isogeometric Design-through-analysis Methodology based on Adaptive Hierarchical Refinement of NURBS, Immersed
Boundary Methods, and T-spline CAD Surfaces. Computer Methods in Applied Mechanics and Engineering,
242-252:116–150, 2012.
[60] P.B. Bornemann and F. Cirak. A subdivision-based implementation of the hierarchical b-spline finite element
method. Computer Methods in Applied Mechanics and Engineering, 253:584 – 598, 2013.
86
[61] H. Speleers, C. Manni, F. Pelosi, and M. L. Sampoli. Isogeometric analysis with Powell-Sabin splines for advectiondiffusion-reaction problems. Computer Methods in Applied Mechanics and Engineering, 221–222:132 – 148, 2012.
[62] S. K. Kleiss, B. Jüttler, and W. Zulehner. Enhancing isogeometric analysis by a finite element-based local
refinement strategy. Computer Methods in Applied Mechanics and Engineering, 213–216:168 – 182, 2012.
[63] D. Burkhart, B. Hamann, and G. Umlauf. Iso-geometric Finite Element Analysis Based on Catmull-Clark :
subdivision Solids. Computer Graphics Forum, 29(5):1575–1584, 2010.
[64] A. Wawrzinek, K. Hildebrandt, and K. Polthier. Koiter’s thin shells on catmull-clark limit surfaces. In Proceedings
of the Vision, Modeling, and Visualization Workshop, Berlin, Germany, 4-6 October 2011.
[65] Hyun-Jung Kim, Yu-Deok Seo, and Sung-Kie Youn. Isogeometric analysis with trimming technique for problems
of arbitrary complex topology. Computer Methods in Applied Mechanics and Engineering, 199(45-48):2796–2812,
2010.
[66] R. Sevilla, S. Fernández-Méndez, and A. Huerta. 3D NURBS-enhanced finite element method (NEFEM). International Journal for Numerical Methods in Engineering, 88(2):103–125, 2011.
[67] R. Sevilla, S. Fernández-Méndez, and A. Huerta. NURBS-enhanced finite element method (NEFEM). International Journal for Numerical Methods in Engineering, 76(1):56–83, 2008.
[68] R. Schmidt, R. Wuchner, and K.U. Bletzinger. Isogeometric analysis of trimmed NURBS geometries. Computer
Methods in Applied Mechanics and Engineering, 241–244:93 – 111, 2012.
[69] C. V. Verhoosel, M. A. Scott, R. de Borst, and T. J. R. Hughes. An isogeometric approach to cohesive zone
modeling. International Journal for Numerical Methods in Engineering, 87(15):336–360, 2011.
[70] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46(1):131–150, 1999.
[71] E. De Luycker, D. J. Benson, T. Belytschko, Y. Bazilevs, and M. C. Hsu. X-FEM in isogeometric analysis for
linear fracture mechanics. International Journal for Numerical Methods in Engineering, 87(6):541–565, 2011.
87
[72] S. S. Ghorashi, N. Valizadeh, and S. Mohammadi. Extended isogeometric analysis for simulation of stationary
and propagating cracks. International Journal for Numerical Methods in Engineering, 2012. In Press.
[73] A. Tambat and G. Subbarayan. Isogeometric enriched field approximations. Computer Methods in Applied
Mechanics and Engineering, 245–246:1 – 21, 2012.
[74] M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J.R. Hughes, and C. M. Landis. A phase-field description of
dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 217–220:77 – 95, 2012.
[75] V. P. Nguyen and H. Nguyen-Xuan. High-order B-splines based finite elements for delamination analysis of
laminated composites. Composite Structures, 102:261–275, 2013.
[76] V. P. Nguyen, P. Kerfriden, and S. Bordas. Isogeometric cohesive elements for two and three dimensional
composite delamination analysis. Composites Science and Technology, 2013. http://arxiv.org/abs/1305.2738.
[77] E. Cervera and J. Trevelyan. Evolutionary structural optimisation based on boundary representation of NURBS.
Part I: 2D algorithms. Computers and Structures, 83:1902–1916, 2005.
[78] R.A.K. Sanches, P.B. Bornemann, and F. Cirak. Immersed B-spline (i-spline) finite element method for geometrically complex domains. Computer Methods in Applied Mechanics and Engineering, 200(13–16):1432 – 1445,
2011.
[79] E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, and A. Düster. Geometric modeling, isogeometric analysis
and the finite cell method. Computer Methods in Applied Mechanics and Engineering, 249–252(0):104 – 115, 2012.
[80] M. Moumnassi, S. Belouettar, E. Béchet, S.P.A. Bordas, D. Quoirin, and M. Potier-Ferry. Finite element analysis
on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces. Computer
Methods in Applied Mechanics and Engineering, 200(5–8):774 – 796, 2011.
[81] G. Legrain. A nurbs enhanced extended finite element approach for unfitted cad analysis. Computational Mechanics, pages 1–17, 2013.
[82] F. Auricchio, L. da Veiga Beirão, A. Buffa, C. Lovadina, A. Reali, and G. Sangalli. A fully locking-free isogeometric
88
approach for plane linear elasticity problems: A stream function formulation. Computer Methods in Applied
Mechanics and Engineering, 197(1-4):160–172, 2007.
[83] T. Elguedj, Y. Bazilevs, V.M. Calo, and T.J.R. Hughes. B-bar and f-bar projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements. Computer Methods
in Applied Mechanics and Engineering, 197(33-40):2732–2762, 2008.
[84] R. Echter and M. Bischoff. Numerical efficiency, locking and unlocking of NURBS finite elements. Computer
Methods in Applied Mechanics and Engineering, 199(5–8):374 – 382, 2010.
[85] T. Elguedj, Y. Bazilevs, V.M. Calo, and T.J.R. Hughes. F-bar projection method for finite deformation elasticity
and plasticity using nurbs based isogeometric analysis. International Journal of Material Forming, 1(1):1091–
1094, 2008.
[86] Rui P. R. Cardoso and J. M. A. Cesar de Sa. The enhanced assumed strain method for the isogeometric analysis
of nearly incompressible deformation of solids. International Journal for Numerical Methods in Engineering,
92(1):56–78, 2012.
[87] R. L. Taylor. Isogeometric analysis of nearly incompressible solids. International Journal for Numerical Methods
in Engineering, 87(1-5):273–288, 2011.
[88] S. Lipton, J.A. Evans, Y. Bazilevs, T. Elguedj, and T.J.R. Hughes. Robustness of isogeometric structural
discretizations under severe mesh distortion. Computer Methods in Applied Mechanics and Engineering, 199(5–
8):357 – 373, 2010.
[89] N. Collier, D. Pardo, L. Dalcin, M. Paszynski, and V.M. Calo. The cost of continuity: A study of the performance
of isogeometric finite elements using direct solvers. Computer Methods in Applied Mechanics and Engineering,
213–216:353 – 361, 2012.
[90] F. Auricchio, L. Da Veiga, T.J.R. Hughes, A. Reali, and G. Sangalli. Isogeometric collocation methods. Mathematical Models and Methods in Applied Sciences, 20(11):2075–2107, 2010.
[91] F. Auricchio, L. Beiro da Veiga, T.J.R. Hughes, A. Reali, and G. Sangalli. Isogeometric collocation for elastostatics
and explicit dynamics. Computer Methods in Applied Mechanics and Engineering, 249–252:2 – 14, 2012.
89
[92] T. Takacs and B. Jüttler. Existence of stiffness matrix integrals for singularly parameterized domains in isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 200(49-52):3568–3582, 2011.
[93] Gang Xu, Bernard Mourrain, Régis Duvigneau, and André Galligo. Parameterization of computational domain
in isogeometric analysis: Methods and comparison. Computer Methods in Applied Mechanics and Engineering,
200(23-24):2021–2031, 2011.
[94] E. Cohen, T. Martin, R.M. Kirby, T. Lyche, and R.F. Riesenfeld. Analysis-aware modeling: Understanding quality
considerations in modeling for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering,
199(5-8):334–356, 2010.
[95] R. Schmidt, J. Kiendl, K.-U. Bletzinger, and R. Wüchner. Realization of an integrated structural design process: analysis-suitable geometric modelling and isogeometric analysis. Computing and Visualization in Science,
13(7):315–330, 2010.
[96] Xianlian Zhou and Jia Lu. Nurbs-based galerkin method and application to skeletal muscle modeling. In
Proceedings of the 2005 ACM symposium on Solid and physical modeling, SPM ’05, pages 71–78, New York, NY,
USA, 2005. ACM.
[97] Y. Zhang, Y. Bazilevs, S. Goswami, C. L. Bajaj, and T.J.R. Hughes. Patient-specific vascular NURBS modeling
for isogeometric analysis of blood flow. Computer Methods in Applied Mechanics and Engineering, 196(2930):2943–2959, 2007.
[98] Y. Zhang, W. Wang, and T.J.R. Hughes. Solid T-spline construction from boundary representations for genus-zero
geometry. Computer Methods in Applied Mechanics and Engineering, 249–252:185 – 197, 2012.
[99] Y. Zhang, W. Wang, and T.J.R. Hughes. Conformal solid T-spline construction from boundary T-spline representations. Computational Mechanics, pages 1–9, 2012.
[100] T. Martin, E. Cohen, and R. M. Kirby. Volumetric parameterization and trivariate B-spline fitting using harmonic
functions. Computer Aided Geometric Design, 26(6):648–664, 2009.
90
[101] M. Aigner, C. Heinrich, B. Jüttler, E. Pilgerstorfer, B. Simeon, and A.-V. Vuong. Swept volume parameterization
for isogeometric analysis. In E.R. Hancock, R.R. Martin, and M.A. Sabin, editors, Mathematics of Surfaces XIII,
volume 5654 of Lecture Notes in Computer Science, pages 19–44. Springer Berlin Heidelberg, 2009.
[102] J.M. Escobar, J.M. Cascn, E. Rodrguez, and R. Montenegro. A new approach to solid modeling with trivariate Tsplines based on mesh optimization. Computer Methods in Applied Mechanics and Engineering, 200(45-46):3210–
3222, 2011.
[103] M.J. Peake, J. Trevelyan, and G. Coates. Extended isogeometric boundary element method (xibem) for twodimensional helmholtz problems. Computer Methods in Applied Mechanics and Engineering, 259:93 – 102, 2013.
[104] R.N. Simpson, S.P.A. Bordas, H. Lian, and J. Trevelyan. An isogeometric boundary element method for elastostatic analysis: 2D implementation aspects. Computers & Structures, 118:2 – 12, 2013.
[105] T.J.R. Hughes, A. Reali, and G. Sangalli. Efficient quadrature for NURBS-based isogeometric analysis. Computer
Methods in Applied Mechanics and Engineering, 199(5-8):301–313, 2010.
[106] F. Auricchio, F. Calabro, T.J.R. Hughes, A. Reali, and G. Sangalli. A simple algorithm for obtaining nearly
optimal quadrature rules for NURBS-based isogeometric analysis. Computer Methods in Applied Mechanics and
Engineering, 249–252:15 – 27, 2012.
[107] A.-V. Vuong, Ch. Heinrich, and B. Simeon. ISOGAT: a 2D tutorial MATLAB code for Isogeometric Analysis.
Computer Aided Geometric Design, 27(8):644–655, 2010.
[108] C. de Falco, A. Reali, and R. Vázquez. GeoPDEs: a research tool for Isogeometric Analysis of PDEs. Advances
in Engineering Software, 42(12):1020–1034, 2011.
[109] D. Rypl and B. Patzák. From the finite element analysis to the isogeometric analysis in an object oriented
computing environment. Advances in Engineering Software, 44:116–125, 2012.
[110] D. J. Benson, Y. Bazilevs, E. De Luycker, M.-C. Hsu, M. Scott, T. J. R. Hughes, and T. Belytschko. A generalized
finite element formulation for arbitrary basis functions: From isogeometric analysis to XFEM. International
Journal for Numerical Methods in Engineering, 83(6):765–785, 2010.
91
[111] R.L. Taylor. A finite element analysis program, programmer manual. http://www.ce.berkeley.edu/rlt/feap/,
2001.
[112] L. Dalcin. PetIGA: A framework for high performance Isogeometric Analysis (using Petsc). https://bitbucket.
org/dalcinl/petiga, 2011.
[113] A. Ratnani and E. Sonnendrucker. Isogeometric analysis in reduced magnetohydrodynamics. Computational
Science & Discovery, 5(014007), 2012.
[114] A. Ratnani. Pigasus, Isogeometric Analysis simulations in Python. INRIA report, 2012.
[115] Rhino. CAD modeling and design toolkit. www.rhino3d.com.
[116] T.J.R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, Mineola, NY, 2000.
[117] B. Nayroles, G. Touzot, and P. Villon. Generalizing the finite element method: Diffuse approximation and diffuse
elements. Computational Mechanics, 10(5):307–318, 1992.
[118] T. Belytschko, Y. Y. Lu, and L. Gu. Element-free galerkin methods. International Journal for Numerical Methods
in Engineering, 37(2):229–256, 1994.
[119] V. P. Nguyen, T. Rabczuk, S. Bordas, and M. Duflot. Meshless methods: A review and computer implementation
aspects. Mathematics and Computers in Simulation, 79(3):763–813, 2008.
[120] J. Chessa.
Programming the finite element method with Matlab.
Northwestern University,
http://www.tam.northwestern.edu/jfc795/Matlab/, 2002.
[121] D. Wang and J. Xuan. An improved NURBS-based isogeometric analysis with enhanced treatment of essential
boundary conditions. Computer Methods in Applied Mechanics and Engineering, 199(37-40):2425–2436, 2010.
[122] T. Rabczuk, R. Gracie, J. H. Song, and T. Belytschko. Immersed particle method for fluid-structure interaction.
International Journal for Numerical Methods in Engineering, 81(1):48–71, 2010.
[123] A. Embar, J. Dolbow, and I. Harari.
Imposing dirichlet boundary conditions with nitsche’s method and
splinebased finite elements. International Journal for Numerical Methods in Engineering, 83(7):877–898, 2010.
92
[124] T. Belytschko, Y. Y. Lu, and L. Gu. Element-free Galerkin methods. International Journal for Numerical
Methods in Engineering, 37:229–256, 1994.
[125] M. Arroyo and M. Ortiz. Local maximum-entropy approximation schemes: a seamless bridge between finite
elements and meshfree methods. International Journal for Numerical Methods in Engineering, 65(13):2167–2202,
2006.
[126] T. P. Fries and T. Belytschko. The extended/generalized finite element method: An overview of the method and
its applications. International Journal for Numerical Methods in Engineering, 84(3):253–304, 2010.
[127] S. Bordas, P. V. Nguyen, C. Dunant, A. Guidoum, and H. Nguyen-Dang. An extended finite element library.
International Journal for Numerical Methods in Engineering, 71(6):703–732, 2007.
[128] P. Laborde, J. Pommier, Y. Renard, and M. Salaun. High order extended finite element method for cracked
domains. International Journal for Numerical Methods in Engineering, 190(47):6183–6200, 2004.
[129] A. Henderson. ParaView Guide, A Parallel Visualization Application. Kitware Inc., 2007.
[130] M. J. Borden, M. A. Scott, J. A. Evans, and T. J. R. Hughes. Isogeometric finite element data structures based on
Bézier extraction of NURBS. International Journal for Numerical Methods in Engineering, 87(15):15–47, 2011.
[131] N. Sukumar, D. L. Chopp, N. Moës, and T. Belytschko. Modelling holes and inclusions by level sets in the
extended finite element method. Computer Methods in Applied Mechanics and Engineering, 190:6183–6200,
2000.
[132] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering, 83(10):1273–1311, 2010.
[133] N. Moës, C. Stolz, P.-E. Bernard, and N. Chevaugeon. A level set based model for damage growth: The thick
level set approach. International Journal for Numerical Methods in Engineering, 86(3):358–380, 2011.
93
[134] S. Natarajan, D. R. Mahapatra, and S. P. A. Bordas. Integrating strong and weak discontinuities without integration subcells and example applications in an XFEM/GFEM framework. International Journal for Numerical
Methods in Engineering, 83(3):269–294, 2010.
[135] T. Elguedj, A. Gravouil, and A. Combescure. Appropriate extended functions for x-fem simulation of plastic
fracture mechanics. Computer Methods in Applied Mechanics and Engineering, 195(7–8):501 – 515, 2006.
[136] J. N. Reddy. An Introduction to the Finite Element Method. McGraw-Hill Education, 3rd edition, 2005.
[137] C. A. Felippa. Advanced Finite Element Methods (lecture notes). http://www.colorado.edu/engineering/
CAS/courses.d/AFEM.d/, 2001.
[138] H. P. Tada, P. C. Paris, and G. R. Irwin. The Stress Analysis of Cracks Handbook. Del Research Corporation,
St. Louis, MO, 1985.
[139] J. Shi, D. Chopp, J. Lua, N. Sukumar, and T. Belytschko. Abaqus implementation of extended finite element
method using a level set representation of three-dimensional fatigue crack growth and life predictions. Engineering
Fracture Mechanics, 77:2840–2863, 2010.
[140] K.Y. Sze, X.H. Liu, and S.H. Lo. Popular benchmark problems for geometric nonlinear analysis of shells. Finite
Elements in Analysis and Design, 40(11):1551 – 1569, 2004.
[141] D. Millán, A. Rosolen, and M. Arroyo. Nonlinear manifold learning for meshfree finite deformation thin-shell
analysis. International Journal for Numerical Methods in Engineering, 93(7):685–713, 2013.
[142] C. A. Felippa. Introduction to Finite Element Methods. Multifreedom constraints II (lecture notes). http:
//www.colorado.edu/engineering/CAS/courses.d/IFEM.d/, 2001.
[143] E. J. Lingen and M. Stroeven. Jem/Jive-a C++ numerical toolkit for solving partial differential equations.
http://www.habanera.nl/.
[144] T.D. Nguyen and G. N. Wells. Geometrically nonlinear formulation for thin shells without rotation degrees of
freedom. Computer Methods in Applied Mechanics and Engineering, 197(33–40):2778–2788, 2008.
94
[145] V. P. Nguyen, P. Kerfriden, M. Brino, S.P.A. Bordas, and E. Bonisoli. Nitsche’s method for two and three
dimensional nurbs patch coupling. Computational Mechanics, 2013. submitted.
[146] V. P. Nguyen, P. Kerfriden, S. Claus, and S.P.A. Bordas. A Nitsche’s method for mixed dimensional analysis:
conforming and non-conforming solid-structure coupling. 2013. In preparation.
[147] T. W. Sederberg.
Computer Aided Geometric Design:
educationportal.html, 2010.
95
Course Notes.
http://www.tsplines.com/