Breitenberger IBRA

Analysis in Computer Aided Design:
Nonlinear Isogeometric B-Rep Analysis of Shell Structures
M. Breitenbergera,∗, A. Apostolatosa , B. Philippa , R. Wüchnera , K.-U. Bletzingera

a Lehrstuhl für Statik, Technische Universität München, Arcisstr. 21, 80333 München, Germany
Abstract
A new concept called analysis in computer aided design (AiCAD) is proposed for design-trough-analysis workflow.
This concept uses non-uniform rational B-Splines (NURBS)-based B-Rep models for the entire workflow. Such
models consist of trimmed NURBS surfaces and are considered standard in the industry, especially for modeling free-
form geometries. The newly developed isogeometric B-Rep analysis (IBRA) used in AiCAD is also presented. IBRA
can be considered as a generalization of isogeometric analysis (IGA) that uses the boundary representation (B-Rep)
of the design model in addition to the same basis functions as in IGA for approximating the solution fields. IBRA
provides the framework for creating a direct and complete analysis model from computer aided design (CAD) in a
consistent finite-element-like manner. Thus, IBRA allows analyzing a CAD model without remodeling and meshing,
even for complex geometries.
For the numerical integration of trimmed surfaces, the concept of nested Jacobian approach (NEJA) with NURBS
surfaces is introduced. In addition, for enforcing the different types of boundary conditions or mechanical entities, a
new finite element type called isogeometric B-Rep element is introduced. Elements of this type permit enforcing, e.g.,
coupling or Dirichlet boundary conditions. A corresponding formulation based on a penalty approach is presented as
well.
The proposed workflow is realized exemplarily for surface modeling and the geometrical nonlinear analysis of
shell structures. The differences between the standard analysis procedure and the AiCAD workflow are explained in
detail. Various numerical examples confirm the accuracy, flexibility, and robustness of the proposed IBRA concept,
thus highlighting its advantages for the realization of design-through-analysis workflow with a uniform geometry
representation.
Keywords:
Analysis in Computer Aided Design, AiCAD, CAD, integration of design and analysis, finite elements, boundary
representation, isogeometric B-Rep analysis, IBRA, IGA, trimmed surfaces, shell structures, nonlinear kinematics,
Nested Jacobian Approach, NEJA, isogeometric B-Rep elements, multi-patches, non-watertightness, non-matching
coupling
1. Introduction
One of the bottlenecks in an efficient virtual product development process lies at the interface between design and
analysis, and it is caused by the use of different geometry representations. In industry the standard output of geometri-
cal modeling in computer aided design (CAD) (see Section 2.1) is a non-uniform rational B-Splines (NURBS)-based
B-Rep model (see Section 3), especially for free-form geometries. The predominant method used in the industry for
solving structural problems is classical finite element analysis (FEA). This method typically uses linear respectively
quadratic polynomials defined over non-overlapping subdomains (the elements) for geometry representation. Thus,
for standard design-through-analysis workflow, a complex geometry transformation called meshing is necessary (see
Section 2.2). The analysis is usually performed using a computer aided engineering (CAE) system.
∗ Correspondingauthor
Email address: [email protected] (M. Breitenberger)
The final publication is available at Comput. Methods Appl. Mech. Engrg. via http://dx.doi.org/10.1016/j.cma.2014.09.033
For design cycles that require the incorporation of the analysis result in the design model, reconstruction of the
geometry model to ensure that it is CAD-suited is necessary. This process is error-prone and, generally, only approx-
imative.
To avoid the meshing, a few alternatives for the integration of design and analysis have been suggested in the past
few years. One of these approaches is the finite cell method (FCM) [1, 2, 3, 4, 5]. This approach has high potential,
especially for three-dimensional structures, and it is being developed continuously by Rank and his co-workers. Other
approaches include the Kantorovich method [6], implemented in Scan&SolveTM [7], or isogeometric analysis (IGA)
introduced by Hughes et al. [8, 9].
IGA entails the use of basis functions for representing CAD geometry as well as for approximating solution
fields. Usually, IGA is based on NURBS [8, 9] because NURBS represent the standard for geometry description
in contemporary CAD systems. NURBS have high continuity, allow for refinement on the identical geometry and
raise the degree of the underlying polynomial functions [10, 11]. NURBS basis functions show very good properties
for analysis. In the past few years, many IGA researchers have demonstrated the advantages of these properties
for analysis purposes, such as robustness and excellent approximation power. IGA has been applied successfully in
many fields of engineering such as structural mechanics [8, 9, 12, 13, 14, 15, 16, 17, 18], fluid mechanics [19], fluid-
structure interaction [20, 21, 22, 23], contact mechanics [24, 25, 26], and optimization [27, 28, 29]. The use of the
basis functions from design also for analysis simplifies the integration of these two tasks greatly. However, a general
and consistent concept for creating direct and complete analysis models from CAD still is not available.
Recent publications have dealt with this shortcoming of IGA and suggested solutions pertaining to specific fields.
Schillinger et al. [3] presented a design-through-analysis workflow that is well suited for solving three-dimensional
problems. Another approach for integrating design and analysis is the use of T-Splines [30, 31] or subdivision surfaces
for geometric modeling and analysis [32, 33, 34, 35]. However, assuming a NURBS-based B-Rep model as standard
output from CAD, these latter approaches also require geometric transformation.
The goal of the design-trough-analysis workflow proposed in the present publication, named analysis in computer
aided design (AiCAD), is to use NURBS-based B-Rep models for the entire workflow, i.e. for geometric modeling,
model analysis and postprocessing (see Section 2). Thus, AiCAD also allows for unification of the design and analysis
models (see Section 2.3), and AiCAD can be fully integrated into an existing CAD system (Section 6). For realizing
the proposed AiCAD workflow, we present the newly developed isogeometric B-Rep analysis (IBRA) method.
IBRA is a general and consequent extension of IGA, which in addition to the basis functions also uses the B-Rep
description of the design model. This is expressed by the extension ”B-Rep” which is the abbreviation of boundary
representation, a well known term in the CAD community [36, 37]. The goal of IBRA is to directly use the design
model, coming from a CAD system, for analysis, without meshing or re-parametrization. Hence, IBRA allows for
the creation of a direct and complete analysis model from the CAD model in a consistent manner, analogous to well
established finite-element-like workflows.
The focus of this publication is on shell structures, so IBRA is explained in detail for shell problems, whereas the
IBRA concept is outlined briefly using one and three-dimensional problems. IBRA for shell structures uses NURBS-
based surface B-Rep models, which consist of trimmed NURBS surfaces for geometry representation. In the past few
years, some investigations on trimmed NURBS analysis have been conducted. Kim et al. [38] presented an approach
to handle membranes, Wang et al. [39] presented a technique for analyzing compound B-Spline surfaces for linear
shell problems, and Schmidt et al. [14] presented an approach for handling shell structures including the treatment of
trimmed multi-patch structures.
The proposed concept is different in the sense that the treatment of trimmed NURBS surfaces is embedded within
a general concept, which finally allows for the direct use of CAD models for analysis. Aspects relevant to practice
such as the application of arbitrarily located loads, local refinement, enforcing weak boundary conditions along trim-
ming curves, dealing with non-matching multi-patches including gaps and overlaps, and the numerical integration
of trimmed surfaces with the nested Jacobian approach (NEJA) using NURBS surfaces are discussed in the present
publication. To analyze NURBS-based B-Rep models in a finite-element-like manner, we introduce a new type of
element called isogeometric B-Rep element.
Since the present publication tries to unify the fields of design and analysis, both fields are reviewed briefly as the
bases for the new proposed methods and techniques.
Thus far, AiCAD workflow is realized within the CAD program Rhinoceros 4, 5 R
[40] and within the integrated
CAD/CAE/CAM program Siemens NX 8.5
R
[41]. This demonstrates the general applicability of the presented
2
concept. Specific implementation aspects are discussed in Section 6.
In order to show the accuracy, flexibility and robustness of IBRA, a few well selected benchmark examples with
increasing complexity (ranging from linear plate problems up to complex nonlinear shell problems) are discussed. An
industrial problem is solved as well, to demonstrate the potential of the proposed method.
The present publication is organized as follows: In Section 2, the new AiCAD workflow is compared with the
standard analysis workflow. In Section 3, some fundamentals of NURBS-based B-Rep models, including an explana-
tion of trimmed NURBS surfaces, are discussed. Some general aspects of the analysis, especially of IGA, are reviewed
in Section 4. In Section 5, the new isogeometric B-Rep analysis is presented, with a special focus on nonlinear shell
analysis. Here as well, the new type of finite element called, the isogeometric B-Rep element is introduced in addition
to a formulation that allows for the coupling of trimmed NURBS surfaces and the enforcing of weak boundary con-
ditions. In Section 6, some aspects concerning the implementation of AiCAD and IBRA are detailed. In Section 7,
a suite of selected benchmark examples with increasing complexity is presented to show the accuracy, flexibility and
robustness of the proposed methods. In the last section, a short conclusion and outlook for future research are given.
2. Design and Analysis Modeling
The design and analysis models are essential components of a design-through-analysis workflow. Thus, in the
following section, a few important aspects of design and analysis modeling are summarized, with special attention on
the representation of the two ”different” models.
In general, a model substitutes an object in a more convenient form that is easier to use and analyze [37]. In the
age of computers and global networks, virtual models play the dominant role in industry and are considered as such
in this publication. A virtual model is a collection of selected information about an object, where the selection and
the level of detail depend on the model’s purpose. In the following list, some important types of models with different
purposes are listed:
– design model (geometry description)
– structural analysis model (geometry description, material properties, boundary conditions)
– rendering model (geometry description, reflection properties, colors, textures)
– cost model (materials, volumes, prices)
In the past few years, industry has strongly trended toward incorporating all information about an object in one
common model for avoiding remodeling and inconsistencies. The problem in merging a structural analysis model
with common uniform model is that within the standard analysis (CAD-CAE) workflow, which is based on classical
finite element analysis (FEA), the design (CAD) and analysis (CAE) models use different geometry representations
(descriptions). Thus, in this publication, an alternative approach called analysis in computer aided design (AiCAD)
workflow is introduced for unifying the design and analysis models. In order to unify both models within the AiCAD
workflow, the analysis model adapts the geometry description of the design model and classical finite element analysis
is replaced with isogeometric B-Rep analysis, a new analysis technique introduced in Section 5.
In the following, the design and analysis modeling processes for standard analysis and the AiCAD workflow are
explained briefly. The different workflows are illustrated in Figure 1 in an abstract manner and are exemplified with
the shell analysis of an oil sump (Figure 2).
2.1. Design Modeling

Design or geometric modeling is used for defining the shape and other geometric characteristics of an object [37].
Depending on the application, several representations and techniques are used because each concept and technique
has its own limitations [42].
Nowadays, geometrical modeling is usually computer based and is performed using computer aided design (CAD)
systems. Thus, in this publication, geometric modeling with CAD systems is assumed. In contemporary CAD sys-
tems, the user can choose between different modeling techniques, which use their own representation e.g. cell decom-
position, sweep representation, constructive solid geometry (CSG), and boundary representation (B-Rep). All of these
representations use a varying different, for its purpose optimal, data structures [37, 36, 43]. Owing to the limitations
of each representation, a transformation to other representations is not always possible [43]. Indeed, a CSG model,
3
defined as a composition of simpler solids (primitives) using Boolean operations [37] can easily be converted into a
B-Rep model (see Figure 3), but in general it is not possible to derive a CSG model from a B-Rep model.
In this publication, the NURBS-based B-Rep model (see Section 3) is assumed as the output of geometric mod-
eling because, this representation is very powerful, included in the standard exchange format IGES [44], and is the
industry standard, especially for free-form geometries. Thus, the NURBS-based B-Rep model is used as the starting
point for both analysis workflows.
2.2. Analysis Modeling – Standard Analysis (CAD-CAE) Workflow

In this publication, analysis modeling refers to the preparation of a model for analysis itself. In the industry,
structural problems are solved predominantly with classical finite element analysis (FEA), which typically uses linear
or quadratic polynomials defined over non-overlapping subdomains (the elements) for describing the geometry of the
analyzed model. Thus, within the standard analysis workflow, illustrated in Figure 1 and Figure 2 , the NURBS-
based B-Rep model needs to be converted into a finite element mesh. This change of geometry description is called
meshing. Note that in almost all cases, mesh generation is just an approximation of the geometry of the design model.
In addition, meshing is a highly complex task in itself (see [45, 46, 47]) because tolerances, non-watertightness [48],
and ”dirty” geometries need to be considered and the meshes have to fulfill certain criteria [49, 47]. Moreover, for
complex models, manual corrections are often necessary, which take a lot of time and require experience. Thus, for
industrial applications 80% of the overall cost of engineering design is for meshing [9].
Another step to get the analysis model is the introduction of analysis-related data, e.g., material properties and
boundary conditions. They can either be applied directly to the finite element mesh, e.g., by selecting the elements or
nodes (option 1 in Figure 1) or to the CAD model (option 2). In option 2, properties are assigned as attributes to the
B-Rep entities (e.g. trimming curves), and after meshing they are transferred to the finite element mesh. This can be
done automatically, and option 2 is thus more convenient for the user, especially when re-meshing of the CAD model
is necessary. The analysis itself is performed using a finite element program.
The steps in the standard analysis workflow can be performed by using either different programs (CAD and FE
program) or by using integrated CAD/CAE system.
In the standard analysis workflow, the analysis model consists of a geometry description, finite element mesh, and
analysis-related data assigned to mesh entities (see Figure 2).
2.3. Analysis Modeling – Analysis in Computer Aided Design (AiCAD) Workflow

In contrast to the standard analysis (CAD-CAE) workflow, in AiCAD the geometry representation of the design
and analysis models is the same. For the analysis, the design model just needs to be refined (see Section 3.2.5), which
does not change the representation of the model but only its description. The refinement is necessary to approximate
the solution with adequate accuracy and can be assigned as a property to the B-Rep entities, e.g., trimmed surfaces
and, trimming curves. In AiCAD, not only the introduction of analysis-related data can be achieved in the CAD
environment (option 2 in the standard analysis workflow) but also the postprocessing of the analysis (see Figure 2).
Here, of course, the CAD system needs to be extended to include the desired postprocessing functionalities, e.g.,
visualization of the stresses on the CAD model. This is usually done by using user interface software (see Section
6.2).
The difference between standard analysis and the AiCAD workflow lies in the representation of the analysis model
and its solution. In classical FEA, the (deformed and undeformed) geometries are represented by a finite element
mesh, whereas in AiCAD, the geometries are represented by a NURBS-based B-Rep model (see Figure 2).
Note that for both workflows, the overall model topology is very important for creating a suitable analysis model.
3. Design – NURBS based B-Rep Model
The proposed AiCAD workflow with its isogeometric B-Rep analysis (IBRA) uses NURBS-based B-Rep models,
commonly used in contemporary CAD systems, for geometry representation (see Section 2.1). Therefore, in the
following text, the NURBS-based B-Rep model composed of trimmed NURBS surfaces, is explained in detail.
4
Standard Analysis Workflow Definition of analysis-
based on classical finite element analysis related data e.g.
option 2 - boundary conditions
- material properties
- ...
CAD geometry kernel option 1 FE geometry kernel
meshing finite element
design model analysis model
mesh
B-Rep
topology of the model
tolerances
non-watertightness
dirty geometries
...
isogeometric refinement
design model analysis model
CAD geometry kernel Definition of analysis- CAD geometry kernel

related data e.g.
- boundary conditions
- material properties Analysis in Computer Aided Design Workflow
- ... based on isogeometric B-Rep analysis
Figure 1: Difference between standard analysis (CAD-CAE) workflow based on classical finite element analysis (FEA) and proposed analysis in
computer aided design (AiCAD) workflow based on isogeometric B-Rep analysis (IBRA)
5
Standard Analysis Workflow Analysis in Computer Aided Design Workflow
based on classical finite element analysis based on isogeometric B-Rep analysis
geometrical modeling
DESIGN
CAD geometry kernel

design model
(NURBS based B-Rep model)
analysis modeling
MESHING REFINING
(DIFFERENT geometry (SAME geometry
representations) design model with analysis-related representation)
data
analysis modeling (preprocessing)
analysis modeling (preprocessing)

ANALYSIS
finite element mesh refined design model
CAD geometry kernel

FE geometry kernel
standard analysis model AiCAD analysis model
classical finite element analysis isogeometric B-Rep analysis

postprocessing
postprocessing
approximated solution represented approximated solution represented

by the finite element mesh by the refined design model
Figure 2: Comparison of standard analysis (CAD-CAE) and proposed analysis in computer aided design (AiCAD) workflow exemplified with
shell analysis of oil sump. Both workflows are realized within the integrated CAD/CAE/CAM system Siemens NX 8.5 R
. The classical finite
element analysis is performed using NASTRAN NX R
[41], whereas the isogeometric B-Rep analysis is performed using the finite element code
CARAT++, developed at the authors’ chair. Within the AiCAD workflow, for transferring data between Siemens NX 8.5 R
and CARAT++, the
interface software TEDA (Tool to Enhance Design by Analysis) is used. TEDA, too, is developed at the authors’ chair.
6
U
b) CSG model: composition of simpler solids

(primitives) using Boolean operations
a) Object
c) NURBS based B-Rep model: composition of

timmed NURBS surfaces
Figure 3: Different geometry representations of same object (after [50])
7
3.1. Boundary Representation (B-Rep)
A boundary representation describes an object by its ”skin” between the ”model” and ”non-model” [43]. The
skin of a three dimensional object is composed of a set of adjacent bounded surface elements called faces. Faces
are bounded by sets of edges, which are curves that lie on the surface of the faces on either side of the edge. The
points, where several faces meet, are called vertices [43, 36]. The main topology entities mentioned above are faces
(F), edges (E), and vertices (V), together with their geometric forms: surfaces (S), curves (C), and points (P) (see
geometry space and topology in Figure 9). A boundary representation of an object consists of two parts: structure
(topology) and shape (geometry) [43]. Depending on the purpose, the face-edge-vertex data kernel is augmented by
additional elements such as shells and/or loops (see topology in Figure 9).
3.2. Non-Uniform Rational B-Splines (NURBS)

NURBS are piece-wise rational basis functions that can be used for parametrizing curves in general as well as for
surfaces and volumes within tensor product construction. They can be used for modeling a large variety of shapes
[11, 51]. The term NURBS stands for Non-Uniform Rational B-Splines and indicates that NURBS are a generalization
of B-Splines. Therefore, a short introduction to B-Splines is given first. For more details about NURBS the reader is
referred to [10, 51]. B-Splines are non-interpolating, piece wise defined polynomial curves. They are defined by the
following:
– set of control points, Pi , i = 1, ..., n
– polynomial degree p
– knot vector Ξ = [ξ1 , ξ2 , ...ξn+p+1 ]
– set of basis functions Ni , i = 1, ..., n
3.2.1. Knot Vector

The knot vector Ξ is a set of parametric coordinates ξi arranged in non-descending order, this set divides the B-
Spline into sections. If all knots are spaced equally, the knot vector is called uniform. NURBS can have non-uniform
knot vectors and are therefore more flexible. Knot values that appear more than once are called multiple knots. The
intervals between two consecutive, distinct knots are called non-zero knot spans. In case the first and last entries in
knot vectors have a multiplicity of p + 1 they are called open knot vectors. They are assumed in the remainder of this
paper.
3.2.2. Basis Functions

B-Spline basis functions Ni are defined by the knot vector Ξ and the polynomial degree p. They can be computed
using the Cox-deBoor recursion formula [10, 51, 52]. The computation starts with p = 0:

1, ξi 6 ξ < ξi+1

Ni,0 (ξ) = 

(1)
0, otherwise

For p > 1 it is
ξ − ξi ξi+p+1 − ξ
Ni,p (ξ) = Ni,p−1 (ξ) + Ni+1,p−1 (ξ) (2)
ξi+p − ξi ξi+p+1 − ξi+1
The basis functions are C ∞ continuous inside a knot span and C p−1 continuous across single knots. At knots with
multiplicity k the continuity of the basis functions is reduced to C p−k . The following list contains some important
properties of the B-Spline basis functions:
– local support, i.e. a basis function Ni,p (ξ) is non-zero only in the interval [ξi , ξi+p+1 ].
– partition of unity, i.e. ni=1 Ni,p (ξ) = 1.
P
– non-negativity, i.e. Ni,p (ξ) > 0.
– linear independence, i.e. ni=1 αi Ni,p (ξ) = 0 ⇔ αi = 0, i = 1, 2, ..., n
P
NURBS geometries have an additional weight for each control point, which can lead to rational basis functions. The
weights are necessary for representing, for example, conical sections exactly. In the case that all weights are equal to
one, the basis functions reduce to simple piece-wise polynomials. The formulas for the NURBS basis functions are
given in Section 3.2.3 and Section 3.2.4.
8
P2 P4 P6
1.0
C(ξ5 )
N7,3
0.8 N1,3
C(ξ6 ) N4,3
N2,3 N3,3 N5,3 N6,3
P3 0.6
0.4
P7
P1 C(ξ7 )
C(ξ1−4 ) C(ξ8−11 ) 0.2
ξ1−4 ξ5 ξ6 ξ7 ξ8−11
0
P5 0 0.2 0.4 0.6 0.8 1.0
(a) B-Spline curve with its control point polygon (dashed lines) (b) Corresponding basis functions defined in parameter space
and knot positions on the curve
Figure 4: Cubic B-Spline curve with the knot vector Ξ = [0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1] and its basis functions
3.2.3. NURBS Curves

B-Spline curves of degree p are computed using linear combinations of the control point coordinates and the
corresponding basis functions. The formula for determining B-Spline curves is given by
n
X
C(ξ) = Ni,p (ξ)Pi (3)
i=1
NURBS curves have non-uniform knot vectors and the additional weight wi for each control point Pi must be consid-
ered. The formula for NURBS curves with their basis functions Ri,p is given as follows:
n n
X Ni,p (ξ)wi Pi X
C(ξ) = Pn = Ri,p (ξ)Pi (4)
i=1 k=1 Nk,p (ξ)wk i=1
A cubic B-Spline with an open knot vector is shown as an example Figure 4. Owing to the open knot vector, the
first and last control points (P1 and P7 ) are interpolated and the curve is tangential to the control point polygon at its
start and end.
3.2.4. NURBS Surfaces

B-Spline surfaces are constructed as tensor products of B-Spline basis functions in the parametric dimensions ξ
and η. They are defined by a net of n × m control points, two knot vectors Ξ and H, two polynomial degrees p and q
and the corresponding basis functions Ni,p (ξ) and M j,q (η). The formula for determining B-Spline surfaces S is
n X
X m
S(ξ, η) = Ni,p (ξ)M j,q (η)Pi j , (5)
i=1 j=1
and NURBS surfaces are defined by
n X
m n X m
X Ni,p (ξ)M j,q (η)wi j Pi j X
S(ξ, η) = Pn Pm = Ri j,pq (ξ, η)Pi j . (6)
i=1 j=1 k=1 l=1 Nk,p (ξ)Ml,q (η)wkl i=1 j=1
An example of a quadratic B-Spline surface is illustrated in Figure 5, which has the knot vectors Ξ = H =
[0, 0, 0, 0.5, 1, 1, 1]. The surface is divided into four parts because it has two non-zero knot spans in each direction.
9
η5−7 ξ5−7
η4 ξ4
η1−3 ξ1−3
(a) B-Spline surface with its control point net (dashed lines) and (b) Corresponding basis functions in parameter space
parameter curves (knots)
Figure 5: Quadratic B-Spline surface with open knot vectors Ξ = H = [0, 0, 0, 0.5, 1, 1, 1] and its basis functions
3.2.5. Geometry Refinement

For the later mentioned purpose of analysis the basis might need to be enriched to properly describe deformed
geometries (see also [8]). There are two ways of refining B-Spline and NURBS geometries. They are called degree
elevation and knot insertion [10]. Both refinement techniques increase the number of control points and thus impart
the capability to represent a greater number of shapes. The refinement itself does not change the initial shape of
the geometry. The continuity of a curve can be decreased by generating knots with multiplicity greater than one
(see Section 3.2.2). In the case of multiplicity equal to (p + 1), the curve gets interpolating at this knot and can be
split easily into two curves. An example of a similar case is shown in Figure 6, which illustrates the knot insertion
refinement of a B-Spline curve by introducing a C 0 continuity at ξ = 0.1. An example of degree elevation is shown in
Figure 7.
For surfaces and solids, for each parameter direction, refinement can be performed in the same manner and inde-
pendently.
3.3. Trimmed NURBS Surfaces

A trimmed NURBS surface is described by a NURBS surface (see Section 3.2.4) and a set of properly ordered
boundary (trimming) curves lying within the parameter space of the surface [53]. Thus, a trimmed surface is a partially
visible surface, defined by the trimmed domain D ∈ R2 (see parameter space in Figure 9) using a B-Rep description.
Trimmed surfaces simplify the representation of complex objects. Examples of trimmed NURBS surfaces are shown
in Figure 9 b) - d).
In general, trimming curves can be of any form; however, when dealing with NURBS entities, it is desirable to
represent these with NURBS. Assume that M boundary curves (see Section 3.2.3) are given by
nk
uk (ξ̃)
" # X
C̃k (ξ̃) = = Ri,l (ξ̃) P̃ki , k = 1, 2, ..., M, (7)
vk (ξ̃)
i=1
where l is the polynomial degree, ξ̃ is the curve parameter, uk and vk are parameters of the surface of the trimming
curve and P̃ki represents the control points of the trimming curve in the parameter space. The curves C̃k (ξ̃) are joined
properly to form outer and inner loops (see topology in Figure 9). The outer loops are oriented counter-clockwise,
whereas the inner loops are oriented clockwise. The curves Ck , which bound the trimmed surface, can be obtained by
mapping the C̃k curves onto the surface [54] (see also Figure 9).
10
P7 P9
N1,3 N4,3
P5 1.0
C(ξ5−7 ) C(ξ8 )
0.8 N10,3
P3 P4 C(ξ9 ) N5,3 N7,3
P6 0.6 N6,3 N8,3 N9,3
P2
0.4
P10
P1 C(ξ10 )
C(ξ1−4 ) C(ξ11−14 ) 0.2
ξ8 ξ9 ξ10 ξ11−14
0
0 ξ5−7 0.2 0.4 0.6 0.8 1.0
P8
(a) Refined B-Spline curve with its control point polygon (dashed (b) Corresponding refined basis functions
lines) and knot positions on the curve
Figure 6: Knot insertion exemplified by B-Spline curve from Figure 4 with knot vector Ξ = [0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1]. After inserting
three knots at ξ = 0.1 the knot vector becomes Ξ = [0, 0, 0, 0, 0.1, 0.1, 0.1, 0.25, 0.5, 0.75, 1, 1, 1, 1] and a C 0 continuity is introduced at the new
knot. The refinement does not change the curve’s shape.
1.0
P3 P6 N11,4
N1,4
P2 C(ξ6 ) P5 P10 0.8
P9 N2,4 N5,4 N7,4 N10,4
C(ξ7 )
P4 0.6 N3,4 N9,4
N4,4 N6,4 N8,4
P7 0.4
C(ξ8 ) P11
P1
C(ξ1−5 ) C(ξ9−13 ) 0.2
P8 ξ1−5 ξ6 ξ7 ξ8 ξ9−13
0
0 0.2 0.4 0.6 0.8 1.0
(a) Refined B-Spline curve with its control point polygon (dashed (b) Corresponding quartic basis functions
lines) and knot positions on the curve
Figure 7: Degree elevation exemplified by B-Spline curve from Figure 4 with knot vector Ξ = [0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1]. The polynomial
degree is increased by one to four. The refinement does not change the curve’s shape.
11
η5−7 ξ5−7
η4 ξ4
η1−3 ξ1−3
(a) Trimmed B-Spline surface with its control point net (dashed (b) Corresponding basis functions defined on trimmed domain
lines)
Figure 8: Trimmed B-Spline surface from Figure 10 and its basis functions defined on trimmed domain, which is described by a boundary
representation (B-Rep).
Ck (ξ̃) = S(uk (ξ̃), vk (ξ̃)), k = 1, 2, ..., M (8)

An example of a trimmed B-Spline surface with its basis functions, defined on the trimmed domain, is illustrated
in Figure 8. The corresponding untrimmed B-Spline surface with its basis functions is shown in Figure 5.
3.4. Trimming
Trimming refers to the operation of creating trimmed surfaces (see Section 3.3), e.g. by performing a Boolean or
blending [50] operation. Mathematically this results in a surface-to-surface intersection (SSI) problem [55, 56, 57,
50], where the trimming curves C̃ in the parameter space are determined.
Figure 10 shows an example of a trimming operation with three intersecting surfaces before and after the execution
of several Boolean operations. The trimming operation consists of the following steps:
– approximation of the intersection curve with the trimming curves C̃k (ξ̃) in the parameter space of each surface
using an SSI algorithm [55, 56, 57]
– determination of the B-Rep description (see Section 3.3 and Section 3.5) for each surface, including specifica-
tions of the trimmed domain (see also topology and parameter space in Figure 9)
– determination of the space curves Ck (ξ̃) by mapping the trimming curves C̃k (ξ̃) onto the corresponding surface
[54]
Note that for the example in Figure 10 the specification of the trimmed domain depends on the type of Boolean
operation. For solving the SSI problem and mapping the trimming curves onto the surfaces most CAD systems use
approximation algorithms, even if the mapping could be done exactly by using polynomials of high degrees (e.g.
p=27) [54]. Indeed, the use of approximation techniques requires the specification of tolerances but it reduces the
computational cost and increases the robustness against numerical instabilities [54].
3.4.1. Trimming Tolerances

Figure 11 shows two different representations of the trimmed B-Spline surface S(1)visible in Figure 9; these were
created by using two different trimming tolerances. The trimmed surface shown in the Figure 11a and Figure 11b
was created in the CAD program Rhinoceros 5 R
by using a trimming tolerance of 10−1 units, while the surface
described in the Figure 11c and Figure 11d was created using a tolerance of 10−8 units. It can be seen that the more
accurately trimmed surface uses a greater number of control points for the trimming (C̃3 and C̃7 ) and space (C3 and
C7 ) curves to approximate the intersection curves.
12
V Vertices
E Edges (direction is not considered)
C Space curves
C̃ Trimming curves
S Surfaces
D Trimmed domain
a) Surface B-Rep model
Geometry space
C(5) (ξ̃)
C(14) (ξ̃)
C (ξ̃)
(6) C(9) (ξ̃) C(15) (ξ̃)
S(1)
visible
(ξ, η)
C (ξ̃)
(4) C (10)
(ξ̃) C (11)
(ξ̃)
S(3)
visible
(ξ, η)
C (ξ̃)
(7)
S(2)
visible
(ξ, η)
C(3) (ξ̃) C(13) (ξ̃)
C(1) (ξ̃) C(8) (ξ̃)
C(12) (ξ̃)
C(2) (ξ̃)
Topology (abstract)
V6
E5 V7 E10
E6 V7,8 V10
Inner loop V5 V12
E7 E12
V1 E4 E7
E13 V11
E3 V4 E9
Face 1 Face 3
Outer loop Outer loop V3
E1 V3 Face 2 E11
E3
V8 Outer loop
E2
E8 V4
V9
V2
(9)
Parameter space C̃ (ξ̃) (15)
C̃ (ξ̃)
(5) D(2)
C̃ (ξ̃)
(12)
(4) C̃ (ξ̃) C̃
(14)
(ξ̃)
(6)
C̃ (ξ̃) C̃ (ξ̃)
(8)
D(1) C̃ (ξ̃) D(3)
(3) (10)
C̃ (ξ̃) C̃ (ξ̃) (13)
C̃ (ξ̃)
(7)
C̃ (ξ̃)
(1)
C̃ (ξ̃) C(2) (ξ̃) (11)
C̃ (ξ̃)
b) Trimmed B-Spline surface 1 c) Trimmed NURBS surface 2 d) Trimmed NURBS surface 3
Figure 9: NURBS-based surface B-Rep model
13
S(2) (ξ, η)
S(1) (ξ, η)
S(3) (ξ, η)
a) Original B-Spline surface 1 b) Original NURBS surfaces 2 and 3
c) Intersecting surfaces before trimming operations
Trimming operation
- approximation of the intersection curve in the parameter
space (C̃ curves) of each surface using an SSI algorithm
- determination of the B-Rep description for each surface
including the specification of the trimmed domain
- determination of the space curves C by mapping the
trimming curves C̃ onto the corresponding surface
d) Intersecting surfaces after trimming operations
Figure 10: Boolean operations including trimming of intersecting surfaces
14
C5
C̃5 C6
C7 C4
C̃4
C̃6
C̃7
C3
C1
C̃3
C̃1 C̃2
C2
(a) Trimming curves in the parameter (b) Corresponding trimming curves in the geometry space deter-
space determined by using a trimming mined by using a trimming tolerance of 10−1 units
tolerance of 10−1 units
C5
C̃5 C6
C7 C4
C̃4
C̃6
C̃7 C3
C1
C̃3
C̃1 C̃2 C2
(c) Trimming curves in the parameter (d) Corresponding trimming curves in the geometry space deter-
space determined by using a trimming mined by using a trimming tolerance of 10−8 units
tolerance of 10−8 units
Figure 11: Two different representations of trimmed B-Spline surface created in Figure 10 based on different trimming tolerances
In CAD the exact intersection curve of two surfaces is represented by the two space curves of the two intersecting
surfaces, where the two space curves are simply an approximation of the exact intersection curve. Thus, there can
exist gaps between the two space curves. Such gaps pertain to non-watertightness (see also [48]). Thus almost all
NURBS-based B-Rep models can be seen as non-watertight models. For example, in Figure 9 the space curves C3
and C7 of the surface S(1) (2)
visible as well as the corresponding space curves C8 of the surface Svisible and C12 of the surface
(3)
Svisible do not match exactly.
Usually, in a CAD system, the user can specify the desired (relative and/or absolute) tolerance, depending on the
application. In Rhinoceros 5 R
[40] for instance the default value for large models is 0.01 units and that for small
models is 0.001 units.
3.5. Topology
In addition the geometric description of surfaces and curves, the topology of the boundary representation needs
to be discussed. Depending on the application, different data structures can be used to describe the topology and
geometry of a boundary representation. Because no data structure is superior to all others, several descriptions can be
used, such as the winged-edge and face-edge-vertex representation [43].
The isogeometric B-Rep analysis is independent of the B-Rep data structure (topology), as long as the data struc-
ture is complete. In this publication the topology is split into two: one for trimmed surfaces and one for the B-Rep
model, which consists of trimmed surfaces.
15
3.5.1. Topology of Trimmed Surfaces
The topology of trimmed surfaces, illustrated in Figure 9 (topology), is already described in principle in Section
3.1 and Section 3.3. Each trimmed surface consists of one face and a set of properly joined edges (E), which define
inner (clockwise oriented) and outer (counter-clockwise oriented) loops.
The vertices (V) with their geometric form the points are used to describe edges. Here, the vertices resp. points
are defined as parameters ξ̃ in the corresponding trimming curve C̃. Thus the geometric description of an edge of each
surface is given by one trimming curve Ck (ξ̃) and two parameters ξ̃. The parameters can be computed using a point
inversion algorithm (see [58]).
For the definition of a trimmed surface, vertices are not really necessary because, usually, their parameters coincide
with the first and last knots of the trimming curves.
For the sake of discussion, in this publication, it is assumed that the orientation of the edges corresponds to the
one of the trimming curves. Thus, a trimmed NURBS surface can be described by defining correctly oriented loops
using trimming curves that are oriented accordingly (see parameter space in Figure 9).
3.5.2. Topology of Surface B-Rep Models

Boundary representations can be used for solid (see Figure 3) and surface models (see Figure 9). In the following
the explanation of topology is restricted to surface models because the focus of this publication is on shell analysis.
In the case of surface models, the object is represented by the ”skin” itself. Thus, such models consist of an
arbitrary number of faces connected through edges, such as face 1 and face 2 with its edge E7 (see topology in Figure
9). This means the topology of surface B-Rep models can be described easily by defining common edges as shown
in Figure 9 (topology). Here, common edges are described by two trimming curves of the adjacent faces and two
vertices, whereas one vertex can be described by the two parameters that are used for defining the start or end point of
an edge. Note that in the case the CAD program provides only the space coordinates of the vertices, the corresponding
parameters for the trimming curves can be computed easily using a point inversion algorithm [58].
4. Analysis
In this section, shell analysis is reviewed briefly as the basis for isogeometric B-Rep shell analysis. In this pub-
lication, the following conventions are used if not indicated differently: Latin indices such as i, j can take the values
1, 2, 3, Greek indices such as α, β can take the values 1, 2. Capital letters refer to the reference configuration and
small letters are used for the description in the current configuration. In addition, derivatives w.r.t. to a quantity r are
abbreviated by (·),r and vectors are expressed in the Cartesian coordinate system ei .
4.1. Analysis of Shells

Shells are thin-walled, curved structures, like car bodies, roof shells or cooling towers. For the analysis of such
structures, shell formulations are attractive because they concentrate on the mechanically relevant effects and thus
the computational effort can be reduced significantly. For solving shell problems two main theories are widely used:
the Reissner-Mindlin (RM) and the Kirchhoff-Love (KL) theory. In this publication the KL theory is used. The
assumptions of this theory are the following:
– the director remains straight
– the director remains perpendicular to the midsurface during the deformation
– no transverse shear deformation is taken into account
– ratio Rt of the characteristic radius R and the thickness t is larger than 20
4.1.1. Geometry Description

The body of a shell can be described by the position vector X resp. x in the reference resp. in the current
configuration as follows:
X(θ1 , θ2 , θ3 ) = Xsurf (θ1 , θ2 ) + θ3 A3 , x(θ1 , θ2 , θ3 ) = xsurf (θ1 , θ2 ) + θ3 a3 , (9)

where Xsurf is the midsurface and θ are the contravariant coordinates of the covariant base vectors Ai resp. ai . Here
i
A1 and A2 resp. a1 and a2 are defined by

16
∂Xsurf ∂xsurf
Aα = , aα = , (10)
∂θα ∂θα
and the third base vector A3 resp. a3 represents the surface normal, which can be computed as follows:
A1 × A2 a1 × a2
A3 = , a3 = (11)
|A1 × A2 | |a1 × a2 |
h i
The corresponding contravariant coordinate θ3 is defined in the range − 2t , + 2t . In the following, the thickness t is
assumed as being constant over the midsurface.
The displacement u between two configurations can be written as
u = usurf = xsurf − Xsurf . (12)
4.1.2. Equilibrium
For arriving at a solution to a structural problem, the internal and external forces need to be in equilibrium. This
can be expressed by the principle of virtual work, which is defined as the sum of the internal and external virtual work
and is given by
δW = δWint + δWext = 0. (13)

The internal and external virtual work in the reference configuration can be formulated as follows [59]:
Z
δWint = − S : δE dΩ
Z Ω Z (14)
δWext = ρ · δu dΩ + tΓ · δu dΓ,
Ω Γ
where Ω is the physical domain, Γ is the domain boundary, S is the 2nd Piola Kirchhoff (PK2) stress tensor, δE is
the energetically conjugated virtual Green-Lagrange strain tensor, ρ is the body force, tΓ is the stress vector at the
boundary, and δu is the virtual displacement. The equilibrium condition in Eq. (13) is fulfilled for an arbitrary virtual
displacement field δu [60]:
∂W
δW = δu = −Rδu = 0 (15)
∂u
Thus, the residual force R must vanish in the case of equilibrium. Here, Eq. (15) is a nonlinear expression that is
linearized in order to solve it with an iterative solution approach such as the Newton-Raphson method. This results in
the following linear expression:
K ∆u = R, (16)
where the tangential stiffness K and the residual force R are used for obtaining the displacement increment ∆u until
Eq. (15) is satisfied.
The equations Eq. (13) and Eq. (14), including appropriate boundary conditions, represent a typical boundary
value problem, which is within the context of structural analysis usually solved with the finite element method using
the same discretization for the geometry and the solution field (isoparametric concept). In a discretized form the
residual force R and the tangential stiffness K can be written as a vector and a matrix, respectively. The components
of the matrix and vector are given by
∂W ∂Wint ∂Wext
Rr = − =− − = Rint
r + Rr ,
ext
(17)
∂ur ∂ur ∂ur
17
∂Rr ∂W 2 ∂Wint
2
∂Wext
2
Krs = =− =− − = Krs
int
+ Krs
ext
, (18)
∂u s ∂ur ∂u s ∂ur ∂u s ∂ur ∂u s
where r, s are indices used for the discretization components. Here, different discretizations schemes can be applied.
Usually, classical finite element analysis (FEA) uses Lagrange polynomials, whereas IGA mostly uses NURBS basis
functions. The proposed IBRA uses the basis functions from IGA just within the trimmed domain (see Section 5).
4.2. Isogeometric Analysis

The idea behind IGA, introduced by Hughes et al. [8], is to use the basis functions used for presenting CAD
geometry also for approximating the solution fields. Usually, IGA is based on NURBS because they represent the
standard for geometry description in current CAD systems. IGA can be applied in many fields such as structural
mechanics [8, 9, 12, 13, 14, 15, 16, 17, 18], fluid-structure interaction [20, 21, 22, 23], contact mechanics [24, 25, 26]
and optimization [27, 28, 29].
For isogeometric shell analysis, NURBS-based element formulations for both theories, KL and RM, are available.
Kiendl et al. [61] developed the KL shell element, whereas Benson et al. [62] and Dornisch et al. [63] developed two
different RM shell element formulations. In addition, Echter et al. [64] developed a hierarchical NURBS-based shell
formulation. In the present paper the geometrically nonlinear formulated NURBS-based KL shell element of Kiendl
et al. [61, 65] is adapted. In the following text, the most important aspects of this shell formulation, necessary for the
proposed isogeometric B-Rep analysis for shell structures, are summarized.
4.2.1. NURBS-based Kirchhoff-Love (KL) Shell Element

The NURBS-based KL shell element is explained in detail in [65]. Considering the assumptions of the Kirchhoff-
Love theory and preintegration over the thickness the residual force vector components Eq. (17) and the stiffness
matrix Eq. (18) are expressed as follows:
∂ε ∂κ ∂ui ∂u j
Z ! XZ X
Rr = n: +m: dA − pi · dA − Fj · , (19)
A ∂ur ∂ur i
∂u r j
∂ur
A
∂n ∂ε ∂2 ε ∂m ∂κ ∂2 κ
Z ! !
Krs = : +n: + : +m: dA + Krs
ext
, (20)
A ∂u s ∂ur ∂ur ∂u s ∂u s ∂ur ∂ur ∂u s
where n, m are the internal membrane forces and moments, respectively. The terms δε and δκ are the corresponding
virtual change in the curvatures and strains within the shell. The variables pi represent the area loads, Fi are the
ext
point loads and Krs represent the external stiffness matrix components, which vanish in the case of deformation-
independent loads. The variable A refers to the surface domain of the shell structure in the reference configuration.
4.2.2. Numerical Integration of NURBS Surfaces

The shell formulation in Section 4.2.1 requires integration over the NURBS surface. Usually, this is done knot
span (element) wise by using a Gauss quadrature rule. Here, p+1 quadrature points in both directions per element
are recommended [8], where p represents the polynomial degree. In the present contribution, this approach is called
standard integration. The area of a knot span i j is given by (see also Figure 14)
Z Z ξend Z ηend Z
|A|i j = dAi j = J1 dξdη = J1 J3 dG, (21)
Ai j ξ start η start G
where the Jacobian J1 indicates the mapping from the geometry to the parameter space and is expressed as
J1 = |A1 × A2 | (22)
18
and the Jacobian J3 indicates the linear mapping from the parameter space to the Gaussian space G ∈ [−1, 1] × [−1, 1]
and is given by
∂ξ ∂η
J3 = (23)
∂ξG ∂ηG
with the parameters ξG and ηG in the Gaussian space. Using the quadrature points l, the area of the knot span i j can
be computed as follows:
ng
X
|A|i j ≈ J1 (ξl , ηl ) J3 (ξl , ηl ) wl (24)
l=0
where ng is the number of quadrature points for each knot span and wl are the quadrature weights.
The mapping J3 (ξl , ηl ) stays constant during nonlinear analysis. Thus, the following weighting factor is introduced:
w̃l = J3 (ξl , ηl ) wl (25)
The weighting factor must be computed just once in preprocessing and is assigned to its quadrature point. Considering
the weighting factor, Eq. (24) can be rewritten as follows:
ng
X
|A|i j ≈ J1 (ξl , ηl ) w̃l (26)
l=0
This is done for generalizing the numerical integration of trimmed NURBS surfaces in the context of IBRA (see
also Section 5.3). In addition, it is more convenient from the implementation viewpoint (see Section 6.3.2).
5. Isogeometric B-Rep Analysis
The proposed AiCAD workflow (see Section 2.3) uses NURBS-based B-Rep models (see Section 3) for geometry
representation. In order to analyze such models, isogeometric B-Rep analysis (IBRA), a new concept, is introduced
here. IBRA extends isogeometric analysis (IGA) [8] by adapting the basis functions as well as the B-Rep description
of the design model (see Section 3.1) for approximating the solution fields. The boundary description is used to
describe trimmed domains (see Section 3.3) in a consistent manner. The IBRA concept allows for analysis of an
entire CAD model in a consistent finite element-like way.
In this section, the numerical integration of trimmed domains resp. geometries is discussed (Section 5.3). A new
type of finite element, the so called isogeometric B-Rep elements or just B-Rep elements, is introduced (see Section
5.1). They are used for enforcing different types of boundary conditions (Neumann, Dirichlet, coupling). In addition,
the section covers the application of arbitrarily located loads (see Section 5.7) and local refinement (see Section 5.8)
within IBRA.
The focus of this publication is on IBRA for shell structures. Thus, this part is explained in detail. The idea of
solving one- and three-dimensional problems with IBRA is briefly discussed in Section 5.9.
5.1. Isogeometric B-Rep Elements

Isogeometric B-Rep elements are finite elements derived from B-Rep entities such as vertices, edges, and faces
(see Section 3.5), which are used to describe one-, two-, or three-dimensional trimmed domains.
They can be used for imposing analysis properties assigned to B-Rep entities within the design-through-analysis
process (see Section 2.3). The following list contains some examples of such properties
– Neumann boundary conditions (e.g. forces, moments)
– Dirichlet boundary conditions (e.g. displacements, rotations,...)
19
Physical space
Ω(1) Ω(2) Γ(2)

n
Γ(1)
c
Γ(1)
d Γ(2)
c
Parametric space
Γ̃(1)
c
Γ̃(2)
n
Ω̃(1)
Ω̃(2)
Γ̃(1)
d
Γ̃(2)
c
Figure 12: Description of a shell with arbitrarily defined boundary conditions within isogeometric B-Rep shell analysis
– Coupling boundary conditions

– Mechanically motivated entities (e.g. cables, beams,...)
– ...
B-Rep elements are defined by non-zero knot spans of embedded geometries such as trimming curves (see Section
3.3). Figure 17 shows examples of B-Rep edge elements. In Section 5.4, they are explained in detail. For B-Rep
elements one distinguishes between non-overlapping and overlapping B-Rep elements. Non-overlapping elements
extract the non-zero basis functions from only one underlying isogeometric element but overlapping elements extract
them from more than one underlying elements (see Figure 17).
The virtual work contributions of B-Rep entities can be external (e.g. Neumann boundary conditions) as well as
internal (e.g. cable element formulation). Therefore they are considered separately in the virtual work expression Eq.
(13) as sum over the different B-Rep contributions:
X
δW = δWint + δWext + δWB-Rep
k
=0 (27)
k
For solving Eq. (27) the terms in Eq. (17) and Eq. (18) need to be enriched by the components of the force vector
and the stiffness matrix of the B-Rep contributions:
∂Wint ∂Wext X δWB-Rep

k X
B-Rep k
Rr = − − − r + Rr +
= Rint ext
Rr (28)
∂ur ∂ur k
∂ur k
∂Wext X ∂WB-Rep
k
∂Wint X
B-Rep k
Krs = − − − = Krs
int
+ Krs
ext
+ Krs (29)
∂ur ∂u s ∂ur ∂u s k
∂ur ∂u s k
20
Svisible
θ2
θ3 A2
C
A3
A1
θ̃1 T3
T2
ωT2
T1
ωT1 θ1
θ 2)
θ 1,
urf (
1 )
(˜θ
Xs
p
x3 -Re
XB
e3
e2
x2
e1
x1
Figure 13: Geometry description of a shell subset Ω(a) in the context of IBRA
5.2. NURBS-based Isogeometric B-Rep Analysis of Shells

In the context of IBRA a shell structure is represented by a surface B-Rep model (see Figure 9), which can
be characterized by open and bounded subsets Ω(a) ∈ R3 and the corresponding boundaries Γ(a) ∈ R3 with a ∈ N
(see Figure 12). Here, each subset Ω(a) and its boundary Γ(a) are represented by a trimmed surface and its space
curves. The trimmed domains of the surfaces and their boundary (trimming) curves are denoted with Ω̃(a) ∈ R2 and
Γ̃(a) ∈ R2 , respectively, in the parameter space (see Figure 12). As shown in Figure 12 the boundaries Γ(a) resp. Γ̃(a)
are decomposed into subsets:
     
[ (a) k  [ (a) l  [ (a) m 
Γ (a)
=  Γ  ∪  Γn  ∪  Γc  ,
d
k l m
      (30)
[ (a) k  [ (a) l  [ (a) m 
Γ̃ (a)
=  Γ̃  ∪  Γ̃n  ∪  Γ̃c  ,
d
k l m
where the following indices are used for expressing the different types of boundary conditions:
– d for Dirichlet boundary conditions (displacement, rotation)
– n for Neumann boundary conditions
– c for coupling boundary conditions
Mechanically motivated formulations (e.g. cables) are no longer considered in the remainder of this paper because
they are not needed for the current discussion.
21
5.2.1. Geometry Description
Each subset of a shell Ω(a) can be described by the position vector of its midsurface Xsurf resp. xsurf (see also
Section 4.1.1). In the context of IBRA the corresponding boundaries are described by the position vector XB-Rep (θ˜1 )
resp. xB-Rep (θ˜1 ) as shown in Figure 13. Here, θ̃1 is the parameter of the trimming curves resp. space curves.
At each point XB-Rep a local orthogonal coordinate system T̃i , aligned with the trimming curve, can be determined
by
∂Xsurf ∂θ1 ∂θ2

T̃2 = = A1 + A2 ,
∂θ̃1 ∂θ̃1 ∂θ̃1
(31)
T̃3 = A1 × A2 ,
T̃1 = T̃2 × T̃3 ,
where T̃2 represents the tangent vector of the space curve, T̃3 is the not-normalized vector of the surface normal and
T̃1 is the vector perpendicular to T̃1 and T̃3 pointing outward from the surface (see Section 3.3 for the correct curve
direction). The length T i of the orthogonal vectors T̃i is defined by
p
T i = Ti · Ti (32)
Thus, the corresponding local Cartesian coordinate system Ti (see Figure 13) can be determined by
T̃i
Ti = (33)
Ti
For the sake implementation, it is more elegant to rewrite the vector T1 resp. t1 with the pseudo parameter θ̃2 by
using the triple product expansion as
∂θ1 ∂θ2 ∂θ1 ∂θ2

T̃1 = A1 + A2 , t̃1 = a1 + a2 , (34)
∂θ̃2 ∂θ̃2 ∂θ̃2 ∂θ̃2
∂θα
where the derivatives ∂θ̃2
can be determined in the reference and current configuration by
∂θ1 ∂θ1
= T̃2 · A2 , = t̃2 · a2
∂θ̃2 ∂θ̃2 (35)
∂θ2 ∂θ2
= −T̃2 · A1 , = −t̃2 · a1
∂θ̃2 ∂θ̃2
5.3. Numerical Integration of Trimmed Surface Domains

In the context of IBRA, the numerical integration of trimmed surfaces is necessary. Analogous to the numerical
integration of untrimmed NURBS surfaces (see Section 4.2.2), for trimmed NURBS surfaces, it is recommended to
integrate knot span (elements) wise as well. Thus, for numerical integration, the trimmed domain D ∈ R2 , which
describes a shell subset Ω̃(a) , is decomposed into subsets Di j ∈ R2 with i, j ∈ N
[[
D= Di j (36)
i j
The knot spans (elements) can be categorized as follows:

– void (no need for integration)
– full (standard integration)
– trimmed
22
Elements that require integration (full and trimmed) are called active elements, whereas void elements are calledinactive
elements. For the numerical integration of trimmed knot spans in the context of structural analysis, several techniques
have been proposed. H.-J. Kim et al. used NURBS-enhanced triangles [38], while Y.-W. Wang et al. used a cell based
approach [39] to parameterize trimmed knot spans for integration. In addition, a quadtree technique, usually used for
the finite cell method (FCM) [1], can be used. The question of the most efficient integration scheme is however still
open even if Nagy et al. [66] presented an integration technique that is optimal in the sense that the final quadrature
points and weights satisfy the moment fitting equations within the trimmed domain up to a predefined tolerance [66].
5.3.1. Nested Jacobian Approach

In this publication, NURBS surfaces are used for parametrizing trimmed knot spans to obtain a partially explicit
description of the trimmed domain. In the context of IBRA, techniques that parametrize the trimmed domain are
referred to as nested Jacobian approaches (NEJA), because they only need to consider an additional mapping (see
also Figure 14). In contrast to Eq. (21), the area of a trimmed knot span i j is given by
Z Z np Z
X ξ̂end Z η̂end np Z
X
|A|i j = dAi j = J1 dDi j = J1 J2 dξ̂dη̂ ≈ J1 J2 J3 dGk (37)
Ai j Di j k=0 ξ̂ start η̂ start k=0 Gk
where n p is the number of surfaces, used for parametrizing the trimmed knot span. The mappings J1 and J3 can be
computed as explained in Section 4.2.2. Here, J3 is the mapping from the parameter space D̂i j to the Gaussian space
G ∈ [−1, 1] × [−1, 1]
∂ξ̂ ∂η̂
J3 = (38)
∂ξG ∂ηG
J2 is the additional (nested) mapping from the parameter space Di j to the parametrization Ŝl of the trimmed element
domain D̂i j . In this publication NURBS surfaces Ŝl are used for parametrizing Di j
np
l
[
Di j = Ŝ (39)
l
The Jacobian J2 is defined by
∂Ŝ ∂Ŝ

J2 = × (40)
∂θ̂1 ∂θ̂2
For the numerical integration of subparameter surfaces the standard Gauss quadrature (see Section 4.2.2) is used,
in which the highest occurring polynomial degree from the trimming curves or surfaces in the trimmed knot span
controls the number of quadrature points.
Implementation aspects of parametrizing trimmed knot spans with NURBS surfaces are discussed in Section 6.
Analogously to the integration of untrimmed knot spans (see Section 4.2.2) the integration of trimmed knot spans
can be reduced to the expression (compare Eq. (26))
ng
X
|A|i j ≈ J1 (ξl , ηl ) w̃l (41)
l=0
where ng is the number of all quadrature point in a knot span and the weighting factors w̃l are given by
w̃i = J2 (ξi , ηi ) J3 (ξi , ηi ) wi . (42)
23
Geometry space
Svisible
S(ξ, η) = Ri j (ξ, η)Pi j

P P
i j Preprocessing
J1 Parameter space
η
J1
N (ξ)M j (η)wi j 1,1,1 R2

e3 Ri j (ξ, η) = P Pi
k l Nk (ξ)Ml (η)wkl
R3
e2 0.75
e1 Di j
0.50
J1
0.25
0,0,0 ξ
0,0,0 0.25 0.50 0.75 1,1,1
J2
Ŝ(ξ̂, η̂) = R̂i j (ξ̂, η̂)Q̂i j

P P
i j J2
Additional (nested) mapping for

trimmed knot spans using NURBS
surfaces
J3
Subpameter space
η̂ η̂
Standard mapping
η̂end R2 η̂end R2
for untrimmed
D̂î, ĵ knot spans
D̂î ĵ
η̂ start ξ̂ η̂ start ξ̂
ξ̂ start ξ̂end ξ̂ start ξ̂end
J3
J3
Gauss space
ηG
+1
ξG Pi j = Control points in geometry space

G
Q̂i j = Control points in parameter space
-1
-1 +1
Figure 14: Overview of spaces and mapping used for NEJA with NURBS surfaces
24
(a) Trimmed NURBS surface with 16 active elements (b) Corresponding location of the quadrature points within the
trimmed domain
Figure 15: Examples of quadrature point locations with respect to IBRA by using NEJA with NURBS surfaces
This can be done because the trimmed domains and, consequently, their mappings J2 and J3 remain unchanged during
nonlinear analysis. This means that, within IBRA the untrimmed and trimmed elements can be treated exactly in the
same manner by considering the weighting factor w̃i . This makes the implementation convenient (see Section 6).
Figure 15 shows an example of the locations of quadrature points for a trimmed surface by using NEJA with
NURBS surfaces.
5.4. B-Rep Edge Elements

A B-Rep edge element is a one-dimensional B-Rep element (see also Section 5.1) embedded into a surface or a
volume description. It is defined by a non-zero knot span of a trimming curve C̃. In this publication, the B-Rep edge
elements are embedded into NURBS surfaces and make use of the basis functions from these surfaces. An example
of basis functions for B-Rep edge elements is shown in Figure 16b. Usually, there is no available explicit description
of these basis functions, but this is not required because only the evaluation at the quadrature points is necessary.
Examples of B-Rep edge elements are shown in Figure 17.
Refinement of the trimming curves (see Section 3.2.5), which increases the number of B-Rep edge elements, influ-
ences only numerical integration because the B-Rep elements are embedded in an underlying geometry. Nevertheless,
refinement is useful for generating non-overlapping elements (see Figure 17), which are beneficial from the assembly
viewpoint. The length L of a B-Rep edge element is given by
∂Xsurf dθ + ∂Xsurf dθ dθ̃1 ,

1 2
Z Z Z !
|L| = dΓe = dXsurf dθ̃1 = (43)
dθ̃1 2 ∂θ1 dθ̃1 ∂θ2 dθ̃1 2
Γe θ̃1 θ̃1
where Γe is the boundary segment of the B-Rep element.

In the following text, a simple B-Rep edge element formulation is presented that allows for the consideration of
coupling and Dirichlet boundary conditions along trimming curves. It is based on a penalty approach and can be used
for geometrically nonlinear shell analysis. Other formulations, for example, those that consider the stress resultants at
the interface, are of course possible.
5.5. Coupling Boundary Conditions – Penalty Approach

In this section, a B-Rep element formulation for coupling two common boundary edges Γ(1)
c and Γc is introduced.
(2)
(α)
The formulation is based on a penalty approach, and it constrains the displacements uB-Rep as well as the rotation
25
C(ξ̃(5) )
C(ξ̃(6) )
6 5
C(ξ̃(4) )
7
C(ξ̃(7)
C(ξ̃(3) ) 3
C(ξ̃(1) )
4
C(ξ̃(2) ) 1
2
(a) B-Rep edges of trimmed B-Spline surface from Figure 10 (b) Basis functions of B-Rep edges extracted from basis functions
of underlying surface
Figure 16: Basis functions of B-Rep elements
Parameter space Geometry space
overlapping B-Rep edge element
non-overlapping B-Rep edge element
S(ξ(ξ̃h ), η(ξ̃h ))
B-Rep edge element defined by

the non-zero knot span [ξ̃h , ξ̃h+1 ]
C̃(ξ̃h ) C̃(ξ̃h+1 ) of the trimming curve C̃ S(ξ(ξ̃h+1 ), η(ξ̃h+1 ))
Figure 17: Examples of isogeometric B-Rep edge elements in parameter and geometry spaces, with arbitrarily chosen discretizations of trimming
curves. The figure shows an example of an overlapping and a non-overlapping B-Rep edge element, as well as the general definition of a B-Rep
edge element.
26
ω(α)
T2 around the local vector T2 to be the same for both boundaries. The virtual work contribution of this element
formulation can be written as follows:
coupling-Penalty disp-Penalty rot-Penalty

δWB-Rep = δWB-Rep + δWB-Rep (44)
disp-Penalty rot-Penalty
where δWB-Rep constrains the displacement and δWB-Rep constrains the rotation at the boundaries to be the
equal. Integration along the boundary is done on Γ(1)c . The integration itself is explained in Eq. (43). For practical
applications, the boundary curves Γ(α) of the two subdomains Ωα do not match exactly (see Section 3.4.1); thus,
for any coupling formulation, a strategy of connecting the master (1) with the slave (2) curve is required. For this
formulation, the minimal distance approach is used.
The virtual work term for constraining the displacement is given as follows:
Z (1)
disp-Penalty
δWB-Rep = −αdisp B-Rep − uB-Rep · δuB-Rep − δuB-Rep dΓc
u(1) (2) (2) (1)
(45)
Γ(1)
c
Thus, the components of the residual force vector and the tangential stiffness matrix are given by
Z
B-Rep disp-Penalty
Rr = αdisp u(1) − u(2) · u(1) (2)
,r − u,r c ,
dΓ(1)
(46)
Γ(1)
c
Z
B-Rep disp-Penalty
Krs = αdisp u(1) (2) (1) (2)
,s − u,s · u,r − u,r c ,
dΓ(1)
(47)
Γ(1)
c
where αdisp is the penalty factor which is selected manually. Guidelines for selecting an appropriate penalty factor are
given in Section 7. The virtual work term for constraining the rotation is given by
Z (1)
rot-Penalty
δWB-Rep = −αrot T2 ± ωT2 δωT2 ± δωT2 dΓc
ω(1) (2) (2) (1)
(48)
Γ(1)
c
and its residual force vector and the tangential stiffness matrix components are as follows:
Z (1)
B-Rep rot-Penalty
Rr = αrot T2 ± ωT2 ωT2 ,r ± ωT2 ,r dΓc ,
ω(1) (2) (2) (1)
(49)
Γ(1)
c
Z (1) (1) (1)

B-Rep rot-Penalty
Krs = αrot ω(1)
T2 ,s ± ωT2 ,s ωT2 ,r ± ωT2 ,r + ωT2 ± ωT2 ωT2 ,rs ± ωT2 ,rs dΓc ,
(2) (2) (2) (2) (1)
(50)
Γ(1)
c
where αrot is the penalty factor for the rotations, which is selected manually as well (see Section 7). The signs in Eq.
(47), Eq. (48) and Eq. (49) depend on how the tangent vectors of the two boundary curves are oriented with respect
to each other. The negative sign is used in the case that both tangent vectors T2 point in the same direction.
The term ωT2 expresses the rotation around the vector T2 (see Figure 13) [67], and it is given by
ωT2 = arcsin(ω · T2 ), (51)

27
where ω is the rotation vector that is defined as follows:
ω = T3 × w, (52)
with w being the displacement of the tip of the vector t3 , given by
w = t3 − T3 . (53)
Because the arcsin function is involved, the rotations at the shell boundaries are restricted to less than 90◦ , unless
angles over 90◦ are tracked separately during nonlinear shell analysis.
The first and second variations of the rotation ωT2 are given by
ω,r · T2
ωT2 ,r = p , (54)
1 − (ω · T2 )2
ω,rs · T2 (ω,r · T2 )(ω · T2 )(ω,s · T2 )

ωT2 ,rs = p + 3 , (55)
1 − (ω · T2 )2 1 − (ω · T2 )2 2
where the variations of the rotation vector ω and the displacement w are given by
ω,r = T3 × w,r
(56)
ω,rs = T3 × w,rs
w,r = t3,r
(57)
w,rs = t3,rs .
The variations of t3 correspond to those of a3 and can be determined directly or found in [65].
Since the penalty factor should be in the range of the Young’s modulus E (see also Section 7.6), a relative penalty
factor is used, which is defined as follows:
α
αrel = (58)
E
For patches with C 0 continuity inside the NURBS description, this B-Rep element formulation can be used for
preserving the relative angle within the patch as well. One simply needs to add parameter curves along the knots with
C 0 continuities.
5.6. Dirichlet Boundary Conditions – Penalty approach

The B-Rep element formulation for Dirichlet boundary conditions can be formally written in the same way, as the
coupling boundary conditions, with the exception that the displacements resp. rotations for one boundary are set to
zero. The virtual work contribution can then be written as follows:
Dirichlet-Penalty disp-Penalty rot-Penalty

δWB-Rep = δWB-Rep + δWB-Rep (59)
where the virtual work term for constraining the displacements is given by
Z (1)
disp-Penalty
δWB-Rep = −αdisp B-Rep · δuB-Rep dΓc ,
u(1) (1)
(60)
Γ(1)
c
28
p
Integration domain
(a) Aritrarily located load p on a B-Spline surface (b) Corresponding integration domain in parameter space
Figure 18: Arbitrarily located loads within IBRA
surface 1
coupling
surface 2 (refined)
coupling
surface 3 (refined)
Figure 19: Local refinement within IBRA. After trimming the surface, the desired regions are refined and the boundaries are coupled.
and that for constraining the rotations is given by
Z (1)
rot-Penalty
δWB-Rep = −αrot T2 δωT2 dΓc .
ω(1) (1)
(61)
Γ(1)
c
The corresponding variations can be computed in a straight forward manner, as shown in Section 5.5.
5.7. Considering Arbitrarily Located Loads

IBRA can be used for considering arbitrarily located loads (see Figure 18). We can define additional one-, two-
or three-dimensional integration domains using a B-Rep description (B-Rep elements) within the parameter space of
the underlying geometry. These integration domains (B-Rep elements) can then be used for applying loads such as
line, surface and volume loads. Integration over these domains can then be performed as explained in Section 5.3 and
Section 5.4.
5.8. Local Refinement

By considering coupling boundary conditions (see Section 5.5) within IBRA quasi arbitrary local refinement can
be realized easily. One simply needs to split the surfaces into different pieces, refine the desired trimmed regions
and couple the boundaries. An example of local refinement within IBRA is shown in Figure 19. Note that this local
refinement differs from those presented in previous publications [68, 69, 70, 71].
29
Trimmed domain
(a) Trimmed domain within NURBS solid cube (b) Trimmed domain within refined NURBS solid cube
Figure 20: Example of trimmed domain of solid model (see also Figure 3)
5.9. Outlook on Isogeometric B-Rep Analysis for Solids and Curves

Because the focus of this publication is on structural shell analysis in this section the idea of applying IBRA to
volumes and curves is introduced. However, the concept is the same as for surfaces: a NURBS-based B-Rep model is
used for the geometry representation throughout the analysis.
Usually, in CAD, a solid is represented by a closed surface B-Rep description, the outer shell (in the sense of a
B-Rep description [36]). To perform an analysis, the inner part of the solid needs to be discretized. In the context
of IBRA, the entire model or individual parts will be embedded into a cube with NURBS volume discretization (see
Figure 20). Inside the created cube, the solid model then defines the trimmed domain of the body.
For numerical integration of the newly created trimmed domain, different techniques can be used for example
the one used for the Finite Cell Method (FCM) [1]. For the application of weak boundary and coupling conditions,
two-dimensional B-Rep elements can be used. Possible formulations that can be used to this end, are given by Ruess
et al. [72, 73], where the formulations are used for the FCM.
Another important step is the representation of deformed geometry. For an IBRA-consistent representation of the
deformed geometry, the displacement field of the cube should be mapped to a refined surface B-Rep description of
the trimmed domain. For surface problems, this mapping would correspond to the mapping of a trimming curve C̃ to
a deformed space curve C (see Section 3.3). In contrast to space curve mapping, which is standard in CAD, surface
mapping is more challenging and requires further research.
Theoretically, IBRA could be used for curves as well because one can define a trimmed domain within the pa-
rameter space of a curve. Practically, this does not make sense because trimming of curves can always be done
exactly.
6. Implementation
In this section, a few important aspects of realizing an AiCAD workflow and extending an IGA code by IBRA are
explained.
The AiCAD workflow has the following three main components:
– design (CAD program)
– interface
– analysis (IBRA code)
Independent of the implementation of AiCAD, the workflow can be seen in an abstract form, as shown in Figure
21. In the following text, the various components and their tasks are explained briefly.
30
AiCAD
Analysis in Computer Aided Design
Interface
user interface e.g. textfiles
pre- and postprocessing provides data
obtains results
CAD IBRA
Computer Aided Design program Isogeometric B-Rep Analysis code
Figure 21: Software components of an AiCAD workflow
6.1. Design – CAD Program

The design component of an AiCAD workflow is realized in a CAD program, where the geometry kernel is used
for the following tasks:
– geometric modeling (see Section 2.1)
– pre- and postprocessing for IBRA
Geometric modeling is the intrinsic function of the CAD program, whereas pre-and postprocessing may be seen
as additional functionalities. The pre- and postprocessing in a CAD program require a user interface, which endows
the CAD software with missing functions, such as
– assigning analysis related data to B-Rep entities (see Section 2.3)
– assigning refinement parameters (see Section 2.3)
– visualizing of deformed geometries (see Figure 62)
– ...
Depending on the CAD program, some of these functions are already available. To give an example, the integrated
CAD/CAE/CAM software Siemens NX 8.5 R
[41] can assign analysis related data to B-Rep entities (see Figure 2),

R
whereas in the CAD program Rhinoceros [40], this needs to be realized by the interface software.
Note that in this publication, it is assumed that a CAD program provides access to geometric data and functions
such as the visualization of NURBS surfaces and space curve mapping (see Section 3.3). From our experience, this
assumption holds for almost all major CAD programs. At the authors’ chair, the AiCAD workflow was realized using
two different programs: Rhinoceros 4,5 R
and Siemens NX 8.5 R
.
6.2. Interface
Interface software is used for communicating between the CAD program and the IBRA code. The interface
software needs to provide among others, the following functions:
– supporting pre- and postprocessing for IBRA in the CAD program
– reading data from the CAD program
– providing data set for the IBRA code
– reading results from the IBRA code
For this publication, the interface is realized with the software TEDA (Tool to Enhance Design by Analysis),
developed at the authors’ chair. TEDA is a C++ library, which has user interfaces for Rhinoceros 4, 5
R
and the

R
integrated CAD/CAE/CAM software Siemens NX 8.5 . Latter interface for TEDA was implemented by Long Chen,
that we gratefully acknowledge.
6.3. Analysis - IBRA Code

The IBRA code is used for solving the IBRA problem, which is given as a set of data for example in a textfile.
This textfile is then provided by the interface software.
31
Start
Read
input data
Preprocessing
Refine geometries
Find intersection points
Loop through elements
Parametrize trimmed domains
Determine quadrature points
Precompute function values
Create isogeometric elements Ke = 0 and Fe = 0
Create B-Rep elements
Loop through
Build connectivities and quadrature points
allocate global arrays
use precomputed basis functions

and derivatives
K=0 and F=0
Add contributions to
Ke and Fe
Solve Kd = F
Write Assemble Ke K
output data and Fe F
Modification for IGA
Additional modification for IBRA Stop
Figure 22: Modifications necessary to extend an existing IGA code by IBRA (overall flowchart according to [9], p.95). The code structure illustrates
the different steps in linear finite element analysis, where K is the stiffness matrix, F is the load vector, d is the displacement vector, and index (·)e
indicates the corresponding term is on element level.
32
Figure 23: Intersection points of trimming curves and (refined) knot lines
The structure of the IBRA code is explained in Figure 22, which shows the modifications necessary for extending
an IGA code by using IBRA. The flowchart also shows the modifications necessary to extend a standard finite element
code to an IGA code [9].
In Figure 22 the changes necessary to extend an FEA code by IGA are highlighted in light grey and the parts
requiring additional modifications for IBRA are highlighted in dark grey. In detail, the following slight modifications
are necessary to the following parts of the code:
– read input data
– preprocessing
– add contributions from the B-Rep elements
– output data
These modifications are explained briefly in the following sections.
6.3.1. Read Input Data

The block ”read input data” has to be extended in order to read the following information:
– B-Rep description of the trimmed domains for all NURBS surfaces (see parameter space in Figure 9 )
– analysis-related data assigned to the B-Rep entities (see Section 2.3)
Note that for IBRA, only the trimmed domain needs to be read, but the corresponding entities in the geometry space,
e.g., space curves are not required (see also Figure 9).
The trimmed domain of a trimmed NURBS surface is simply described by one outer loop and an arbitrary number
of inner loops. The loops, in turn, consist of properly ordered lists of correctly oriented trimming curves (see also
Section 3.3).
The analysis-related data such as boundary conditions or material properties are assigned to the B-Rep entities (e.g.
trimmed domains and trimming curves).
For coupling boundary conditions, a pair of trimming curves (master and slave) needs to be specified. In addition,
the start and end parameters as well as their relative orientation to each other (see also Section 3.5.2) is required. For
setting up the master-slave relation (see minimal distance approach in Section 5.5), the space curves need to be within
a certain tolerance. This can be checked by the CAD program before exporting the IBRA input.
After these datasets are available in the IBRA code, preprocessing is required for creating finite elements. In
principle, this preprocessing could be done in the interface software as well.
6.3.2. Preprocessing
The preprocessing step is needed to create finite elements for IBRA. It is recommended to solve the following
tasks with the help of a basic CAD kernel/library or to shift these tasks to the interface software, which already has
access to a CAD kernel. The following sequence is recommended:
– apply the geometry refinement (see Section 3.2.5) that is assigned to the surfaces
– find all intersection points of the trimming curves and the (refined) knot lines (see Figure 23 and Figure 24a)
33
ξ̃end C0
ξ̃ start C0
ξ̃
(a) Find intersection points of trimming curve with the knot lines (b) Insert C 0 continuities in trimming curve at intersection points
ξ̂
η̂
η̂
η̂
η̂
ξ̂ ξ̂
(c) Create Surface (d) NURBS Patch
Figure 24: Algorithm for parametrizing a trimmed knot span with NURBS surfaces
34
η̂
ξ̂
η̂ η̂ η̂
ξ̂ ξ̂ ξ̂
(a) Parametrization with a negative Jacobian mapping (b) Improved parametrization
Figure 25: Parametrization of special cases
– introduce a C 0 continuity in the trimming curves at each intersection point (see Section 3.2.5 and Figure 24b).
– identify active elements (see Section 5.3) and eventually assign proper flags
– extract the parts of the trimming curves between the intersection points (in the meanwhile they have C 0 conti-
nuity) and assign them to the corresponding elements
– take the assigned pieces of the trimming curves and use them to create NURBS surfaces for parametrizing the
trimmed domain (see Figure 24c and Figure 24d)
– use parts of the trimming curves between the intersection points to create B-Rep elements (see Section 5.4)
– determine the locations of the quadrature points, their weighting factors (see Eq. (42)), and required basis
functions; assign them to the corresponding isogeometric elements
– determine quadrature points locations, weighting factors, and required basis functions of the trimming curves
and the underlying geometries (master and slave); assign them to the corresponding B-Rep elements
Note that for better integration of B-Rep elements with a coupling formulation, clipping of the master curve might
be necessary to ensure that it conforms to the slave curve (see also Section 3.2.5). For this publication, for each
quadrature point of the master curve, a corresponding point on the slave side is assigned, and the minimal distance
approach is employed (see also Figure 55).
Once the information pertaining to all elements is gathered, the remainder of the procedure is exactly the same as
that in classical finite element analysis.
6.3.3. Special Cases

In the parametrization of the trimmed elements, there might occur cases where the standard procedure (see Figure
24c) fails, as is shown in Figure 25a. The figure shows a case where parametrization leads to a negative Jacobian. In
this case, the knot span is divided into four sub cells, which are used as the new basis for parameterizing the trimmed
knot span. Such subdivisions have to be repeated until a non-negative Jacobian is obtained and the trimmed domain
is completely parameterized.
6.3.4. Output Data

The output of IBRA basically consists of the displacement field of the active nodes, meaning the ones that influence
the trimmed domain. The displacement field can be used for determining the deformed geometry and stresses. In
principle, the output of IBRA is the same as of IGA, with the exception that the description of the trimmed domain
needs to be considered. The description is exactly the same for the output as it is for the input. For visualizing
deformed geometries in a CAD program, the following data are sufficient:
– description of NURBS geometry (e.g. NURBS surface)
– B-Rep description of trimmed domain (see parameter space in Figure 9)
35
– displacements values of active control points resp. nodes
The space curves are then mapped in the CAD program (see Section 3.3). Note that the inactive control points do not
have displacement because they have no influence on the trimmed domain.
IBRA has been realized in CARAT++, a finite element software, developed at authors’ chair.
6.3.5. Postprocessing
As mentioned above, for the visualization of trimmed surfaces, description of NURBS surfaces and B-Rep de-
scription of the trimmed domain is required (as in the input file) in addition to the displacements of the active control
points. The stresses can be evaluated in exactly the same manner as in IGA with the additional specification that they
are defined only within the trimmed domain.
6.3.6. Flying nodes

In IBRA, control points with very small influence on the trimmed domain are termed flying nodes, because it is
possible that the displacements of these nodes change during nonlinear analysis. Nevertheless, the solution converges
owing to the negligible contribution of the flying nodes on the trimmed domain resp. solution.
These flying nodes are closely related to the control points in the fictitious FCM domain, which also make a minor
contribution to the solution. According to [69], these points should move as freely as possible to permit a smooth
extension of the physical solution into the fictitious domain. This rationale can be applied to the flying nodes as well.
For structural problems, the small contribution of the flying nodes has no influence on the solution. This is illustrated
by the benchmark example in Section 7.5, which is a complex nonlinear shell problem. For the sake of visualization,
it is advantageous to eventually identify these flying nodes for speeding up the process of visualization and sorting
out meaningless information.
7. Benchmarking and Investigations
To highlight the accuracy, flexibility and robustness of IBRA, some well selected benchmark examples with in-
creasing complexity are discussed. To demonstrate the correctness and to quantify the errors of the results, the results
are compared with solutions that are either computed analytically, with IGA, or with the commercial finite element
software ANSYS R
[74]. The benchmark examples range from geometrically linear plates to nonlinear shells, includ-
ing weak boundary conditions. For all examples, a linear elastic, isotropic material is used.
This section can be separated into two main parts. The first part deals with problems that are discretized with one
trimmed NURBS surface for investigation of NEJA with NURBS surfaces (see Section 5.3.1). The second part deals
with isogeometric B-Rep elements (according to Section 5.1), specifically with enforcing weak boundary conditions.
At the beginning, the example of a simple plate in membrane action is solved for showing the accuracy of NEJA
with NURBS surfaces. The second benchmark deals with a geometrically linear (i.e. linear kinematics) plate in
bending, which is used for demonstrating the robustness of NEJA against distortions related to geometric description.
In the last two nonlinear shell examples within the first block, flexibility in terms of trimming tolerances and the
accuracy of the proposed NEJA technique is demonstrated.
The first example of the second block deals with a geometrically nonlinear cantilever subjected to an end moment,
a so called mainspring. This benchmark is used for discussing weak coupling boundary conditions for different
discretizations, influence of ”dirty” geometries on the result, and choice of appropriate penalty factors (according
to Section 5.5). The second example deals with a nonlinear shell discretized with two trimmed surfaces. The final
benchmark example is the analysis of an industrial oil sump. Moreover, the examples highlight the simplicity and
straightforwardness of analysis model preparation using IBRA.
The geometries for the benchmark examples are created with the CAD program Rhinoceros 5 R
. The trimming
tolerance (according to Section 3.4.1) used for the examples is denoted separately for each benchmark example. For
some few examples, it would be possible to represent the benchmark geometries exactly, i.e. without any trimming
error, but this is not the goal of the proposed method. The goal of IBRA is to use the CAD model ”as it is”. Thus,
for all presented benchmark examples, the geometries are created in a CAD system, and the geometries are analyzed
directly. The analysis itself is performed with the in-house finite element program CARAT++. The results of the
classical finite element analysis are computed in ANSYS R
[74] by using quadrilateral-dominant meshes with bi-linear
36
y
Exact traction
Parameters:
Geometry: L = 4.0 R = 1.0
Symmetry BC
Exact traction
Material: E = 105 ν = 0.3
Far-field traction: T x = 10.0
L
Tx Analytical solution:
Tx R2 Tx R4 R2
σrr (r, θ) = 2 (1 − r2 ) + 2 (1 + 3 r4 − 4 r2 )cos(2θ)
Tx R2 Tx R4
σθθ (r, θ) = 2 (1 + r2 ) − 2 (1 + 3 r4 )cos(2θ)
R r 4 2
σ·n=0 σrθ (r, θ) = − T2x (1 − 3 Rr4 + 2 Rr2 )sin(2θ)
θ x
Symmetry BC
Figure 26: Benchmark example of infinite plate with hole - problem description and exact solution
quadrilaterals and triangular Reissner-Mindlin (RM) shell elements. For IBRA, the NURBS-based Kirchhoff-Love
(KL) shell element [61] is used. Note that KL theory is not always mentioned explicitly. The postprocessing is realized
with the interface software TEDA (see Section 6.2) plugged into the CAD program Rhinoceros 5 R
[40].
7.1. Definition of Error Measures

In order to verify the proposed technique, a few quantitative error measures need to be defined for the benchmark
examples. The error over a domain can be expressed by the sum of the L2 -norms of each component of the stress
tensor as follows:
X
er = k σi exact − σi num kL2 (Ω) , (62)
i
where i represents the components xx, yy and xy. In cases where a computed result (e.g. displacement) at a specific
point is compared to a reference solution, the following relative error measure is used:
|uref − uIBRA |
er = (63)
|uref |
For a better illustration of the convergence, the relative difference between different geometry refinements is used.
This difference should decrease steadily for a correctly formulated finite element analysis. The measure of this relative
difference is defined as follows:
|ufine − ucoarse |
dr = (64)
|ufine |
For a general comparison of the relative difference between two results (ref 1 and ref 2) the following formula is
used:
|uref1 − uref2 |
dr = (65)
|uref1 |
7.2. Infinite Plate with Circular Hole

The first benchmark example deals with the two-dimensional linear elastic problem shown in Figure 26. This
setup is used to solve an infinite plate with a circular hole under in-plane tension at infinity. The exact solution to this
problem is given in [75], and Hughes et al. solved the problem with IGA [8]. This problem had already been solved
by Kim et al. in [38] using a trimmed NURBS surface. Thus, with IBRA using NEJA with NURBS surfaces, one
should obtain the same results as published in [38].
37
(a) Untrimmed NURBS surface description used for IGA (b) Trimmed NURBS surface description used for IBRA
Figure 27: Different geometric descriptions (trimmed and untrimmed NURBS surfaces) of the geometry in the benchmark example, shown in
Figure 26. The dashed lines show the control point polygon, and the black points represent the control points.
The presented results (see Figure 29 and following) clearly show that the approximated solutions with trimmed
NURBS surfaces converge to the exact solution with increasing geometric refinement. It can be noticed that the
convergence properties of this example (using the geometry definitions shown in Figure 27b) are not as good as those
of an untrimmed surface (see Figure 27a). The reason for this difference will explained later in detail.
To compare the IGA and IBRA results, as well as the analytical solution, the stress distributions σθθ along the hole
are plotted (see Figure 29), and the error as the sum of the L2 -norm of each stress tensor component (see Eq. (62))
over the domain is investigated (see Figure 31).
The different discretizations (trimmed and untrimmed surfaces) for IGA and IBRA are shown in Figure 27. The
untrimmed surface used for IGA employs the same definition as proposed in [8], and the trimmed surface used for
IBRA is created using the CAD program Rhinoceros 5 R
with a trimming tolerance of 10−8 units (see Section 3.4.1).
For this specific example, one could define the trimmed NURBS surface without any trimming error by using an exact
NURBS description of the trimming curve, but this is not common when using a CAD system, and doing so would
contradict the proposed philosophy of integrating design and analysis (see Section 2.3). Note that due to the reduced
modeling effort, the standard approach in practice would be to create a trimmed surface even if it does not represent
the geometry in Figure 27 exactly. However, classical IGA requires complicated and costly geometric models suitable
for analysis (see Figure 27a). In contrast, IBRA allows the use of the trimmed surface directly for analyis without any
modifications. One just needs to refine the surface according to the desired solution accuracy, which is a requirement
of mechanics, not of the workflow. Different levels of refinement for the initial discretizations (see Figure 27) are
shown in Figure 28. In both cases, the existing knot spans in the parameter space are subdivided sequentially. In the
step from mesh Figure 28a to Figure 28b, only one direction is divided in order to ensure that the number of elements
in both directions is the same. Figure 28 graphically summarizes the refinement strategy for untrimmed and trimmed
surfaces. These discretizations are the basis for comparative evaluation of the analysis results of both approaches.
Figure 29 shows a comparison of the stress distribution σθθ at r = 1 along the boundary for different geometry
refinements shown in Figure 27. Each plot shows the analytical, IGA, and IBRA solutions using a specific number
of elements. The polynomial degree is p=q=3 for each plot. In Figure 29a, the untrimmed discretization (see Figure
28c) fits better than the one from IBRA (see Figure 28g). This is because of the fact that for IGA, the NURBS surface
parametrization coincides with the boundary curve. Therefore, approximation power at the boundary is very good,
whereas the trimmed discretization has just one element to approximate the solution at the hole boundary. Owing to
the global definition of the basis functions (see Figure 30b) along the boundary, the solution to the trimmed version
oscillates around the analytical solution. In contrast, the basis functions of the untrimmed discretization have local
support (see Figure 30a), which allows for better approximation.
The same effect can be observed for geometries with 8x8 elements (see Figure 28d and Figure 28h). Again, for
38
(a) IGA geometry 1 (b) IGA geometry 2 (c) IGA geometry 3 (d) IGA geometry 4
(e) IBRA geometry 1 (f) IBRA geometry 2 (g) IBRA geometry 3 (h) IBRA geometry 4
Figure 28: Different (refined) geometry descriptions of the example shown in Figure 26
this level of refinement, the trimmed discretization has only three elements for approximating the stresses at the hole
boundary, as opposed to eight elements for the untrimmed case. In the trimmed case result, one can again observe
small oscillations around the analytical solution. With a greater number of elements, the approximations become
better and finally converge to the exact solution. The conclusion is that the untrimmed discretization uses better-suited
parametrization for this specific problem. Thus, it is not surprising that the trimmed discretization does not perform
as well as the untrimmed one.
To confirm overall convergence to the exact solution, the sum of the L2 -error norm of each stress tensor component
over the domain (see Eq. (62)) is plotted for different refinements in Figure 31. This figure shows the error norm
versus the largest diameter found in the mesh and the number of degrees of freedom (DOF). It can be seen that the
plots of the untrimmed and trimmed discretization cases are similar. As expected, the untrimmed discretization shows
faster convergence. This is because the largest deviations for the stresses are in the region around the hole, where the
untrimmed version is optimally refined, as explained before.
The plot in Figure 31a shows the same behavior as the one presented in [8], but the geometry refinement is
different. In the present work, the parameter space is refined uniformly, whereas in [8], the geometry was refined
uniformly.
The results of the plots in Figure 31c and Figure 31d are similar to those published in [38]. The only distinction
is the different integration technique.
Similar plots obtained by using the B-Spline version of FCM have been published by Schillinger et al. [3]. In the
case of FCM, the parameter space coincides with the geometry space, which means J2 (see Figure 14) is constant.
Thus, for this example, the only difference between the B-Spline version of FCM and IBRA lies in how the integration
of the visible domain is performed. Hence, the IBRA results are similar to the ones published by Schillinger et al. [3].
With local refinement in the relevant regions, the convergence rate of the trimmed discretization could be increased
(see also [3]).
From this example, it can be concluded that the untrimmed discretizations better match the resolution requirement
of the problems mechanics than do the trimmed discretizations. The price to be paid for this ideal discretization is
greater modeling effort. Nevertheless, very accurate results can be achieved even with the trimmed NURBS surface.
7.3. Plate in Bending with Circular Hole

This benchmark example deals with a linear (i.e. linear kinematics) plate in bending. The setup is described in
Figure 32. The goal of this example is to investigate the robustness of IBRA to distortions in the geometric description.
Note that the trimmed geometries of the previous example have no such distortions because the parameter space and
the geometry space coincide, i.e., the mapping J1 (see Section 5.3.1) is constant over the entire domain.
39
Analytical solution Analytical solution
Tangential stress σθθ at hole boundary

Trimmed Trimmed
30 30
Untrimmed Untrimmed
20 20
10 10
0 0
-10 -10
π π 3π π π π 3π π
0 8 4 8 2 0 8 4 8 2
Angle θ Angle θ
(a) 4x4 elements (b) 8x8 elements
Analytical solution Analytical solution

Trimmed Trimmed
30 30
Untrimmed Untrimmed
20 20
10 10
0 0
-10 -10
π π 3π π π π 3π π
0 8 4 8 2 0 8 4 8 2
Angle θ Angle θ
(c) 16x16 elements (d) 32x32 elements
Figure 29: Stress distribution of σθθ at hole boundary r = 1 for different refinements. Compared are results of untrimmed and trimmed discretiza-
tion.
40
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
(a) Basis functions of hole boundary for untrimmed discretization (b) Basis functions of hole boundary for trimmed discretization
shown in Figure 28c shown in Figure 28g
Figure 30: Comparison of the basis functions at the hole boundary for different discretizations
Independent of the initial geometric description, the results of the refined geometries should converge to the same
solution. Thus, the distortion of the geometric description must not influence the (converged) results, but only lead to
a difference in the convergence rate, which is shown in the following.
The example described in Figure 32 is solved for different geometric descriptions, where one set of descriptions
has a geometric distortion. This geometric distortion is shown in Figure 33, where the red curve (CP 2 - CP 8)
represents a parameter curve defined in the middle of the parameter space. It is used to illustrate the distortion of
the geometric description. By intention, the distortion is quite substantial, which allows for a demonstration of the
robustness of NEJA using NURBS surfaces to numerical errors in the J1 mapping.
The different geometric descriptions for this benchmark are shown in Figure 34, where the trimming was per-
formed with a tolerance of 10−8 units. The geometries in Figure 34 (a-d) are described by four untrimmed NURBS
surfaces, those in (e-h) are described by simple trimmed NURBS surfaces, and those in (i-l) are described by trimmed
NURBS surfaces with a distorted geometric description (see Figure 33).
The problem (see Figure 32) is solved for the geometric descriptions shown in Figure 34 and for some further
refinements following the same logic. For the NURBS-based KL shell element, a polynomial degree of p=q=3 is
used. In addition the problem is solved with classical finite elements by using quadrilateral-dominant meshes with
bi-linear RM elements. Note that it is not possible to solve this problem using only one untrimmed geometry because
the KL shell element requires C 1 continuity between elements. Therefore, bending strips have to be used [76].
Figure 35 shows the convergence of displacement in the z-direction at point A (see Figure 32) for the four different
types of discretizations. The displacement at point A converges for the KL theory to uz = −6.3499 and for the RM
theory to uz = −6.3748. This results in a difference of 0.0249, which corresponds to a 0.3% discrepancy; this can be
explained by the different mechanical models.
Moreover, the graph shows that the trimmed discretizations perform the best. Moreover, the IBRA solution with
high distortion of the geometric description converges very fast. In addition it can be seen that the untrimmed descrip-
tions with the bending strips converge from above, but to the correct solution. For the bending strip factor, a value of
5 is used, which corresponds to a penalty factor of 105 · E (Young’s modulus). According to Kiendl et al. [76], this is
in the optimal range.
For highlighting the good convergence behavior of IBRA, the convergence of the displacements at points A and E
(see Figure 32) of the simple trimmed IBRA geometries are investigated. The convergence graphs of the displacement
in z-direction for these two points are plotted in Figure 36 by using different polynomial degrees . The convergence
behavior of point E is shown for illustrating the very good point-wise convergence achieved at points on the boundaries
as well.
41
1 1
p=2
p=3
0.01 0.01 p=4
p=5
Error in L2 norm
Error in L2 norm
0.0001 0.0001
1e-006 1e-006
p=2
1e-008 p=3 1e-008
p=4
p=5
1e-010 1e-010
0.1 1 10 10 100 1000 10000
hmax Number of DOFs
(a) Error measured in L2 -Norm
versus largest element dimension (b) Error measured in L2 -Norm versus number of degrees of free-
found in untrimmed discretization (IGA) dom for IGA (untrimmed)
1 100
p=2
0.1 10 p=3
1 p=4
0.01 p=5
Error in L2 norm
0.1
Error in L2 norm
0.001 0.01
0.0001 0.001
0.0001
1e-005 p=2
p=3 1e-005
1e-006 p=4
1e-006
p=5
1e-007 1e-007
0.1 1 10 100 1000 10000
hmax Number of DOFs
(c) Error measured in L2 -Norm versus largest element dimension (d) Error measured in L2 -Norm versus number of degrees of free-
found in trimmed discretization (IBRA) dom for IBRA (trimmed)
Figure 31: Convergence of sum of L2 -Norm of each stress tensor component in Cartesian coordinates x,y and xy (see Eq. (62))
undeformed geometry
clamped f
E = 107 mN2
C B ν = 0.0
t = 0.1m
R E
f = 10.0 m
N
L
L = 10.0m
y E
z R = 3.0m
y A
D x A
deformed geometry
L x
Figure 32: Benchmark example of a plate in bending - problem description
42
CP7 (0, 10) CP8 (5, 10) CP9 (10, 10)
C B
CP4 CP6 CP5

(0, 5) (10, 5) (15, 5)
y
D x A
CP1 (0, 0) CP2 (5, 0) CP3 (10, 0)
Figure 33: Distortion of geometric description of benchmark example shown in Figure 32
(a) Untrimmed geometry 1 (b) Untrimmed geometry 2 (c) Untrimmed geometry 3 (d) Untrimmed geometry 4
(e) Trimmed geometry 1 (f) Trimmed geometry 2 (g) Trimmed geometry 3 (h) Trimmed geometry 4
(i) Trimmed geometry with a dis- (j) Trimmed geometry with a dis- (k) Trimmed geometry with a (l) Trimmed geometry with a dis-
torted geometric description 1 torted geometric description 2 distorted geometric description 3 torted geometric description 4
Figure 34: Different geometry descriptions of geometry shown in Figure 32
43
-6.60
-6.50
-6.40
-6.30
z-displacement
-6.20
-6.10
-6.00
untrimmed
-5.90 trimmed
trimmed distorted
-5.80
bi-linear RM
-5.70
101 102 103 104 105 106
Number of degrees of freedom
Figure 35: Displacement of point A in z-direction (see Figure 32) using different discretizations. For the NURBS-based discretizations a polyno-
mial degree of p=q=3 is used.
-6.40 -1.66
-6.20 -1.64
-6.00 -1.62
z-displacement
z-displacement
-5.80 -1.60
-5.60 -1.58
-5.40 p=2 -1.56 p=2

p=3 p=3
-5.20 p=4 -1.54 p=4
p=5 p=5
-5.00 -1.52
101 102 103 104 105 101 102 103 104 105
Number of degrees of freedom Number of degrees of freedom
(a) Displacement in z-direction at point A (b) Displacement in z-direction at point E
Figure 36: Convergence plots of displacement in z-direction for points A and E shown in Figure 32
44
100 100
p=2 p=2
10−1 p=3 10−1 p=3
relative error in z-displacement

relative error in z-displacement
10−2 p=4 10−2 p=4

p=5 p=5
10−3 10−3
10−4 10−4
10−5 10−5
10−6 10−6
10−7 10−7
10−8 10−8
10−9 10−9
101 102 103 104 105 101 102 103 104 105
(a) Relative error in z-direction at point A (b) Relative error in z-direction at point E
Figure 37: Relative error of displacement in z-direction (see Eq. (63)) for points A and E shown in Figure 32
Figure 36 shows that the displacement in z-direction converges very quickly for both points. For better evaluation
of the convergence behavior, the two corresponding error plots are shown in Figure 37. As reference, an IBRA
solution with a polynomial degree of p=q=3 and 64x64 elements is used. The plots show the relative error, as defined
in Eq. (63). The polynomial degrees p=q=2 and p=q=3 yield optimal slope. However, for the polynomial degrees
p=q=4 and p=q=5, the convergence rate is not as good as it should be. This could be ascribed to the fact that the
solution is not compared to the exact one but only to a very accurate solution. Note that for this case, the relative error
is smaller than 10−8 .
A typical deformed geometry within IBRA is shown in Figure 38, where the result of the plate in bending exam-
ples using simple trimmed geometry is given. The points in the middle of the hole have no influence on the solution
and they are thus called inactive.
The conclusion of this example is that in contrast to the previous benchmark, the trimmed geometries perform
better and yet can be modeled with greater ease. It is shown that IBRA is very robust against distortions in the
geometric description owing to the use of a NEJA technique with NURBS surfaces. Moreover, the very good point-
wise convergence achieve at points on the boundaries as well, is shown.
7.4. Untrimmed Double Curved Shell

This example is used for demonstrating the accuracy of numerical integration for the proposed NEJA technique
with NURBS surfaces (see Section 5.3) in the context of IBRA. The problem description is given in Figure 39, where
the geometry is described by an untrimmed NURBS surface. For this example fully nonlinear kinematics are used.
For investigating the accuracy of numerical integration, the problem is solved with the following three settings, and
the results are compared:
– IGA with one untrimmed surface and standard integration (see Section 4.2.2)
– IGA with one untrimmed surface and p+1 quadrature points per direction
– IBRA with one surface and two complementary trimmed domains (see Figure 39); to this end, p+1 quadrature
points per direction (see Section 5.3) are used
The dashed curve in Figure 39 shows the interface of the two trimmed domains. Different settings are used for
both the computation of the stiffness matrix and the load vector. For this example, the displacements in the z-direction
at point C are compared.
Figure 40a shows the relative difference (see Eq. (65)) in displacement along the z-direction at point C using IGA
with p+1 resp. p+2 quadrature points per direction. It can be seen that the relative difference with fewer degrees of
freedom (DOFs) is around 10−3 , whereas the relative difference decreases very quickly to 10−11 with a greater number
of DOFs.
45
inactive nodes
active nodes
Figure 38: Deformed geometry of example shown in Figure 32, including control point net and examples of active and inactive nodes resp. control
points
D B dead load fz on whole structure D

E = 107 mN2
L
2 ν = 0.2
fixed support
z t = 0.1m
A C fz = −100.0 mN2
x L = 10.0m
A
D C R = 3.0m
B
R undeformed geometry
L
B D
C
y
A B deformed geometry
A C
x
L
Figure 39: Description of benchmark example with untrimmed double curved shell. The radius R is used to split the integration domain into two
domains
46
1.00e-002 1.00e-002
1.00e-004 1.00e-004
1.00e-006 1.00e-006
relative difference
relative difference
1.00e-008 1.00e-008
1.00e-010 1.00e-010
1.00e-012 1.00e-012
p=2 p=2
1.00e-014 p=3 1.00e-014 p=3
p=5 p=5
1.00e-016 1.00e-016
101 102 103 104 105 101 102 103 104 105
(a) Relative difference of displacements in z-direction (see Eq. (b) Relative difference of displacement in z-direction (see Eq.
(65)) at point C using IGA with p+1 resp. p+2 integration rule (65)) of point C in Figure 39 using IGA with p+1 integration rule
resp. IBRA with two complementary trimmed domains
Figure 40: Relative difference in displacement along z-direction (see Eq. (65)) of point C (see Figure 39) by using different integration schemes
For the sake of comparison, Figure 40b shows the relative difference (see Eq. (65)) in displacement along the z-
direction at point C, as obtained using IGA and IBRA with two integration domains. For IGA, a standard integration
rule is used, whereas for IBRA, the NEJA technique with NURBS surfaces and the standard integration rule is used
(see Section 5.3). For the trimming curve, a polynomial degree of p=3 is used. In Figure 40b, it can be observed that
the difference with fewer DOFs is in the same range as the one shown in Figure 40a. For a high number of DOFs, the
relative difference is in the range of 10−15 (i.e., around machine precision).
The conclusion of this example is that the NEJA technique with NURBS surfaces, as it is used for IBRA, is very
accurate and, thus, well suited for integrating trimmed surfaces. In order to increase the efficiency of NEJA with
NURBS surfaces the number of quadrature points should be decreased, which requires further research (see also [77]
and [66]).
7.5. Double Curved Shell with Circular Hole

The last example of the first part deals with a geometrically nonlinear double-curved shell problem. This example
illustrates that IBRA works for complex geometries with challenging mechanics as well. In addition, the load (see
Section 5.7) is applied to the trimmed domain. The benchmark example is described in Figure 41. Easy modeling of
such geometry with untrimmed NURBS surfaces is almost impossible, which means that a broader concept such as
IBRA should be used.
For this example, the influence of trimming tolerances (see Section 3.4.1) on the analysis result is discussed.
Here, trimming is performed with a high tolerance of 10−2 units, which is higher than the standard CAD tolerance for
such an example, as well as with a very low tolerance of 10−7 units. The results for the different trimming tolerances
are compared for judging from the analysis perspective, whether trimming with a tolerance lower than the standard
tolerance used in CAD is effective or whether doing so would simply lead to unnecessary computational cost. To
verify the IBRA results, they are compared with the solutions of classical finite element analysis.
Note that for IBRA, geometry refinement can be done very easily and rapidly (see Section 3.2.5), whereas for
classical FEA, non-trivial meshing is necessary. Figure 42 shows the difference between the geometries used for
IBRA and typical classical finite element discretization. The latter has been created using ANSYS R
, a commercial
software package.
For better understanding the benchmark example and illustrating the geometrical non-linearity of the problem, the
load-displacement curve for the z-displacement of point C is shown in Figure 43. It can be seen that the structure
becomes much softer after reaching the load factor of approximately 0.55. The deformed geometry for this problem,
obtained using IBRA, is shown in Figure 44a. The corresponding undeformed geometry is shown in Figure 42a.
47
D B dead load fz on whole D
structure E = 107 mN2
L
2 ν = 0.2
fixed support
z t = 0.1m
A C fz = 100.0 mN2
x L = 10.0m
A
D C E R = 3.0m
B
R E undeformed geometry
L
B D
C
y
A B deformed geometry
A C
x
L
Figure 41: Benchmark example of a double curved trimmed shell with a circular hole - problem description
(a) Geometry used for IBRA (b) Typical finite element mesh used for the classical FEA
Figure 42: Comparison between an IBRA geometry and a finite element mesh used for solving the problem in Figure 41
1.00
0.80
Load factor
0.60
0.40
0.20
0.00
0.00 0.50 1.00 1.50 2.00 2.50 3.00

z-displacement
Figure 43: Load-displacement curve for point C, shown in Figure 41
48
(a) Deformed IBRA geometry using one trimmed NURBS surface (b) Deformed IBRA geometry using two coupled trimmed
NURBS surfaces with different non-matching refinements
Figure 44: Examples of nonlinear IBRA solutions of problem described in Figure 41 resp. Figure 59
-3.00 -0.25
-2.50 -0.20
-2.00 -0.15
z-displacement
z-displacement
-1.50 -0.10
-1.00 -0.05
p=2
p=3 p=2
-0.50 0.00
p=5 p=3
Bilinear RM p=5
0.00 0.05
101 102 103 104 105 106 101 102 103 104 105
(a) Displacement in z-direction of point C (b) Displacement in z-direction of point E
Figure 45: Displacement in z-direction of points C and E (see Figure 41) using trimming tolerance of 10−2 units
49
100 101
p=2 p=2
p=3 p=3
10−1 100
p=5 p=5
10−2 10−1
relative difference
relative difference
10−3 10−2
10−4 10−3
10−5 10−4
10−6 10−5
101 102 103 104 105 101 102 103 104 105
(a) Relative difference in displacement of point C along z-direction (b) Relative difference in displacement of point E along z-direction
for different polynomial degrees for different polynomial degrees
Figure 46: Relative difference (see Eq. (64)) in displacement of points C and E along z-direction (see Figure 41) using trimming tolerance of 10−2
units
1.00e-003
p=2 p=2
p=3 p=3
p=5 p=5
1.00e-001
relative difference
relative difference
1.00e-004
1.00e-002
1.00e-005 1.00e-003
101 102 103 104 105 101 102 103 104 105
(a) Relative difference for point C (b) Relative difference for point E
Figure 47: Relative difference (see Eq. (65)) in displacement of points C and E along z-direction in Figure 41 using trimming tolerance of 10−2
and 10−7 units
50
z, w
M E = 1.2 · 106
x, u h ν = 0, L = 12
b = 1, h = 0.1
L b Mmax = 2πEI/L = 50π/3
M = 0.35Mmax
0.8
End moment (x π50/3)
wtip
−utip
0.6 M = 0.7Mmax
0.4 M
0.2
esh
exact ial m
Init
0
0 5 10 15 20 M = Mmax
Tip deflections
Figure 48: Benchmark example of mainspring - problem description
To show the point-wise convergence of this benchmark example, the two points C and E (see Figure 41) are
investigated. Figure 45 shows the displacement in the z-direction at points C and E for different polynomial degrees.
For these plots, a trimming tolerance of 10−2 units is used. A comparison of the results obtained using different
trimming tolerances will be shown later.
Figure 45a additionally shows the convergence of classical FEA using quadrilateral-dominant meshes with bi-
linear and triangular RM elements. It can be seen that classical FEA needs approximately 2000 DOFs for yielding
an accurate solution, whereas with IBRA employing p=q=3 good results can be achieved with only 400 DOFs.
Nevertheless all solutions converge to almost the same result (see Figure 45). For the displacement along the z-
direction at point C, the IBRA results converge to −2.66029, whereas the FEA results converge to a slightly higher
value of −2.6973. Again, the discrepancy can be explained by use of different shell theories.
For demonstrating the good convergence behavior of IBRA, Figure 46 shows the relative difference between the
results of sequentially refined geometries (see Eq. (64)). For both plots the relative differences decrease steadily, which
means that the IBRA results for both points converge well. The relative difference for point C using a polynomial
degree of p=q=5 and 10000 DOFs is less than 10−5 , which underlines the good convergence properties of IBRA
with NEJA using NURBS surfaces. The relative difference at point E is almost the same. This shows that point-wise
convergence over the whole domain is achieved, even at points on the boundary curve.
In order to show the influence of trimming tolerance on the IBRA results, the example is also solved for various
geometries with two different trimming tolerances, namely 10−2 and 10−7 . Figure 47 shows the relative difference
(see Eq. (65)) between the results at points C and E by using two different trimming tolerances. For point C the
relative difference is less than 10−4 and for point E it is less than 2 · 10−2 . This clearly indicates that the influence of
the trimming tolerance, i.e., using a low resp. high trimming tolerance is negligible w.r.t. the IBRA result.
The conclusion of this example is that the analysis results are not sensitive to the trimming tolerance. It can be
stated that the trimming tolerance in CAD is definitively accurate enough for structural analysis. Note that the mesh
for classical FEA approximates the trimming curves in cases where the trimming error is already introduced. The
error coming from the mesh approximation is an additional source of inaccuracies. This latter source of error is not
introduced in an IBRA geometry. In addition, this example demonstrates the good point-wise convergence behavior
of IBRA, even for single points on the boundary.
51
Strong Dirichlet boundary condition Weak Dirichlet boundary condition
Single patch
1 1
12 12
Matching patches
1 1
2 10 2 10
Non-matching patches
1 1
2 10 2 10
Straight-trimmed patches
1 1
2 10 2 10
Curve-trimmed patches
1 1
2 10 2 10
Patches with a gap
1 1
2 10 2 10
Arbitrarily-trimmed
patches
- Not applicable - 1
12
Figure 49: Different geometry descriptions of problem in Figure 48. The points represent characteristic control points. For all multi-patch
geometries, a weak coupling boundary condition is used. Both types of Dirichlet boundary conditions, strong and weak, can be used.
52
Single patch Matching Non-matching Straight-trimmed Curve-trimmed Patches Arbitrarily-trimmed
patches patches patches patches with gap patches
Figure 50: IBRA solutions of problem in Figure 48 with M = Mmax using different geometry descriptions from Figure 49
10 10
p=2 p=2
1 p=3 p=3
1
p=5 p=5
0.1
0.1
z-displacement
z-displacement
0.01
0.001 0.01
0.0001
0.001
1e-005
0.0001
1e-006
1e-007 1e-005
10 100 1000 10 100 1000
Number of elements Number of elements
(a) Error of displacement in z-direction (middle of patch) for the (b) Error of displacement in z-direction (middle of patch) for the
geometry ”non-matching patches” in Figure 49 geometry ”straight-trimmed patches” in Figure 49
Figure 51: Error of displacement in z-direction (middle of tip) for the problem in Figure 48 by using different geometry descriptions with M = Mmax
7.6. Cantilever Subject to End Moment

In the first benchmark example of the second part several aspects of weak Dirichlet and coupling boundary con-
ditions are discussed. IBRA coupling and Dirichlet boundary conditions are enforced by using the new proposed
isogeometric B-Rep edge elements (see Section 5.1) following the element formulation from Section 5.5 resp. Sec-
tion 5.6. The proposed formulation is purely geometry based and uses a penalty approach. Thus, it requires the
estimation of penalty factors. For this publication, only one globally defined factor is used, which is applied for both
constraining the displacements and rotations at the common boundaries. Similar to case of bending-strips [76], the
proposed formulation achieves good results in this case for wide range of penalty factors. However, in contrast to the
case of the bending-strips, the proposed formulation can deal with non-matching patches along trimming curves. For
the following examples a relative penalty factor of αrel = 100 (see Section 5.5) is used. This choice of the relative
penalty factor will be explained later.
The benchmark example itself treats the mainspring problem. The problem description is given in Figure 48 (see
also [78]). This example is used for demonstrating the flexibility of IBRA from the perspective of the problem’s
geometry description. Thus, the problem (see Figure 48) is solved with IBRA for the different geometry descriptions
53
0.8
0.8 elements=24 elements=24
0.7
elements=48 elements=48
elements=96 0.6 elements=96
0.6 elements=192 elements=192
0.5
z-displacement
z-displacement
0.4
0.4
0.3
0.2
0.2
0.1
0 0
-0.1
10−4 10−2 100 102 104 106 10−4 10−2 100 102 104 106
Relative penalty factor αrel
rot Relative penalty factor αrel
rot
(a) Displacement in z-direction (middle of tip) using polynomial (b) Displacement in z-direction (middle of tip) using polynomial
degree of p=q=2 for different elements degree of p=q=3 for different elements
Figure 52: Displacement in z-direction for solution of problem in Figure 48 using ”matching patches” geometry shown Figure 49 with M = Mmax
versus relative penalty factor for rotations. The displacement is strongly coupled for both investigations.
shown in Figure 49. Here, for all geometries, a polynomial degree of p=q=4 is used. The multi-patch geometries
use weakly enforced coupling boundary conditions (see Section 5.5) with a relative penalty factor of αrel = 100 at
common boundaries. The Dirichlet boundary conditions can be enforced strongly or weakly, where strong boundary
conditions are restricted to the patch boundaries (untrimmed surface). The geometry ”arbitrarily-trimmed patches”
demonstrates the necessity of weakly enforced boundary conditions along trimming curves. The acting moment is
also applied with isogeometric B-Rep elements referring to the rotation at the boundary.
The IBRA solutions corresponding to the different geometries are shown in Figure 50. All of them represent the
expected solution for a closed ring with M = Mmax .
For detailed investigation of the solution quality, two convergence plots are discussed in Figure 51. The two
plots show the displacement of the tip (point in the middle) in the z-direction for the ”non-matching” and ”straight-
trimmed” geometries (see Figure 49). It can be seen that both geometries have an absolute error of less than 10−4 in
displacement along the z-direction for equilibrium accuracy of 10−6 . Note that a small discrepancy is introduced by
using the penalty approach.
For discussing the influence of the penalty factor on the solution Figure 52 shows the z-displacement versus the
relative rotation penalty factor. Here, the displacements at the common boundary are coupled strongly to isolate the
effect of rotation coupling. It can be observed that a relative penalty factor larger than 10 gives good results regardless
of the polynomial degree.
For quantifying the error of rotation at the boundary, the relative difference between the angles at the common
edge versus the relative penalty factor is plotted in Figure 53. Here, again, the displacement coupling is enforced
strongly. It can be observed that by increasing the relative penalty factor, the error decreases uniformly, independent
of the polynomial degree. Thus, it can be stated that the higher the penalty factor, the smaller will the error be. Note
that a too high penalty factor could lead to numerical problems, see also [65].
For the considered problem, Figure 53 shows that for a relative penalty factor of 100, the difference in the angles
is less than 10−6 .
The same behavior can be observed for weak displacement coupling. Thus, a globally defined penalty factor of
100 is used in this publication because it can be considered sufficiently high for achieving good results and not too
high for avoiding numerical problems.
For practical applications, IBRA needs to deal with non-watertight models because almost all CAD models are
non-watertight as explained in Section 3.4.1. Figure 54 shows the solution of an extremely non-watertight model (see
”patches with gap” geometry in Figure 49) for the problem shown in Figure 48. For the proposed B-Rep element
54
1
p=2
p=3
0.01
p=4
Difference of the angles
0.0001
1e-006
1e-008
1e-010
1e-012
10−4 10−2 100 102 104 106
Relative penalty factor αrel
rot
Figure 53: Difference of angles (between common boundary) for solution of the problem in Figure 48 using ”matching patches” geometry (see
Figure 49). In this study, the displacement is strongly coupled.
overlap
gap
Figure 54: IBRA solution of problem in Figure 48 using extremely non-watertight model (see ”patches with gap” geometry in Figure 49)
55
L = 1.37 + 10.63 + f
1.37 10.63
f
overlap gap Closest points
gap situation Master

Slave
f
Closest points
overlap situation
Figure 55: Problem description for investigating non-watertight (gaps and overlaps) models used in the problem shown in Figure 48 with the
”curved-trimmed patch” geometry with additional gap description shown in Figure 49. The variable f expresses non-watertightness in the hori-
zontal direction. For the coupling formulation, the closest points are used. The master-slave relationship for the gap and overlap situation is shown
in detail.
formulation (see Section 5.5), the gaps are simply preserved, as can be seen in Figure 54. The displacements and
rotations at the master and slave boundaries are the same. Thus, the overlap in Figure 54 corresponds to the gap size
which at its turn is kept constant throughout the analysis. Note that such a non-watertight model is an extreme case
and that all analysis techniques (FEA, IGA, FCM, ...) need to interpret such (”dirty”) geometries (adding material,
preserving the gap,...).
Figure 55 investigates the influence of non-watertightness (see Section 3.4.1) on the analysis result. It shows
the ”curve-trimmed patches” geometry shown in Figure 49 with the gap-overlap variable f . The variable is used to
describe the size of the gap resp. overlap. A positive value expresses a gap, whereas a negative value indicates an
overlap.
Figure 56 shows the influence of the gap-overlap function f on the difference in the displacement along the
z-direction, which is given by
d1 = |uTip ( f ) − uref |, (66)

where uref is the z-displacement when f = 0. For the curve d1 , the total length L = 1.37 + 10.63 + f is used as position
for applying the moment.
In Figure 56, it can be seen that the displacement difference in the final result for a small | f |, especially within
the CAD tolerance of 0.001 units (see ”small geometries” Section 3.4.1), is negligible. The reason for the different
influence of gaps resp. overlaps is shown in Figure 55. In this example, given that the coupling technique is based
on the minimal distance approach (master-slave), the behavior for gaps is better because there is a unique assignment
of the quadrature points for the master and slave curves. For the overlaps, the assignment of quadrature points is not
unique on the interface.
Figure 57 shows the moments m1 in the longitudinal direction, as evaluated at the boundaries of the IBRA solution
for the ”straight-trimmed patch” geometry shown in Figure 49. They are represented by straight lines perpendicular
to the deformed geometry, whereas the length of these lines indicates the magnitude. The moment is constant over the
entire domain, as expected.
56
0.04
d1
0.035
Difference for z-displacement
CAD tolerance (0.001)

0.03
0.025
0.02
0.015
0.01
0.005
0
-6 -4 -2 0 2 4 6
f · 103
Figure 56: Influence of non-watertightness on analysis result of problem in Figure 48 with M = Mmax using ”curved-trimmed patch” geometry
with additional gap shown in Figure 49 with a gap/overlap of f compared with a gap/overlap of f = 0. The default CAD tolerance of 0.001 units
is represented by the light grey colored bar.
Figure 57: Evaluation of stress resultants at boundaries of two non-matching trimmed NURBS surfaces, here for moment m1 in longitudinal
direction. The stress resultants are represented by straight lines perpendicular to the deformed geometry. As expected, m1 is constant along the
ring.
57
F
fixed
Figure 58: IBRA for shell structures with kinks. The kink is preserved by the proposed B-Rep element formulation (see Section 5.5)
D B
dead load fz on whole E = 107 mN2
L structure D
2 ν = 0.2
master
z t = 0.1m
A C fixed support fz = −100.0 mN2
x L = 10.0m
D C R1 = 3.0m
R2 A
α R2 = 5.0m
B R3 = 7.0m
R1 M α = 45◦
L
R3
undeformed geometry
y slave C
A B
x deformed geometry
L
Figure 59: Benchmark example of double-curved trimmed shell including weak coupling boundary conditions - problem description
The conclusion of this benchmark example is that IBRA is demonstrably very flexible from the viewpoint of the
geometric description of a problem, regardless of how ”stupid” the geometric description may look like. Of course,
a simple geometry description is generally better suited for analysis. Nevertheless, IBRA provides the framework
to handle complex geometries as well. In addition, it can be shown that IBRA can treat non-watertight models in a
robust manner.
7.7. L-Shaped Cantilever

The second benchmark example of the second part is used for demonstrating that the proposed coupling formu-
lation (see Section 5.5) can preserve kinks as well. Because kinks are treated in exactly the same way as explained
above this example is not evaluated quantitatively. Deformation of the L-shaped cantilever is shown in Figure 58.
7.8. Double Curved Shell with Circular Hole including Weak Coupling
The third benchmark example of the second part is used for testing the B-Rep element formulation (see Section
5.5) for a mechanically complex problem. The problem is exactly the same as that presented in Section 7.5, albeit
with a different geometry representation. Instead of using a single-patch geometry, two weakly coupled non-matching
patches are used (see Figure 59).
58
-2.70
-2.68
z-displacement
-2.66
-2.64
single p=5
-2.62 multi p=3
multi p=4
multi p=5
-2.60
102 103 104 105
Number of degrees of freedom
Figure 60: Z-direction displacement for problem shown in Figure 59 using different polynomial degrees. For the coupled multi-patches, the
single-patch solution can be used as reference. Therefore, the figure contains the single-patch solution with p=q=5 from Figure 45a as well.
A A
E = 2.1 · 105
load fz ν = 0.2
t=5
W fz = 10
L = 440
W = 240
H = 50
z
H
x
y
L fixed
Figure 61: Problem description of oil sump example. The dark grey region is loaded. The Dirichlet boundary conditions is applied to the edge with
the dashed curve.
For this example, the slave patch uses 2 more elements per direction than the master patch. The displacement
of point C in the z-direction is shown in Figure 60. It can clearly be seen that the multi-patch results converge to
the same solution as those for the single-patch (see Figure 45a). A comparison between the deformed single- and
multi-patch solutions is shown in Figure 44.
The conclusion of this example is that the proposed B-Rep element formulation works for complex nonlinear
(double curved) shell problems as well producing accurate results regardless of the discretization.
7.9. Industrial Application: Oil Sump

The very last example, described in Figure 61, deals with the structural analysis of an oil sump. This illustrative
example is used for demonstrating that IBRA can be used for industrial problems. The oil sump examples consist
of more than 50 trimmed surfaces and about 120 common edges. Figure 62 shows the comparison of the deformed
discretizations for the classical FEA (realized with Nastran NX in Siemens NX 8.5 R
) and IBRA. Both geometries
represent the problem solution, but the corresponding geometry representations are different. IBRA uses a CAD
geometry description, whereas classical FEA uses a finite element mesh. The colors indicate the magnitude of the
displacement.
The conclusion of this example is that IBRA provides a general framework for analyzing complex CAD models
(industrial applications) in a consistent manner by decomposing the model such that a finite element procedure can be
59
(a) FEA solution in Siemens NX 8.5
15.1
0.0
Section A-A
undeformed fixed
fixed
deformed
(b) IBRA solution in Rhinoceros 5 and section A-A with undeformed and scaled deformed geometries (see Figure 61)
Figure 62: Comparison of FEA and IBRA solutions of oil sump example described in Figure 61
60
applied.
8. Conclusion and Outlook
The new Analysis in Computer Aided Design design-through-analysis workflow is presented. This scheme allows
for unification of the design and the analysis models. In addition, the new concept of isogeometric B-Rep analysis
(IBRA) is introduced as an extension of isogeometric analysis (IGA). It is shown that IBRA is a necessary and a
consequent generalization of IGA for efficient integration of design and analysis models. The IBRA concept provides
a finite-element-like framework for analyzing complex CAD models in a consistent manner (see Section 6).
IBRA is successfully applied to CAD models (NURBS-based B-Rep models), with the Nested Jacobian Approach
(NEJA) with NURBS surfaces being used for numerical integration of the trimmed NURBS surfaces. The new type of
elements called isogeometric B-Rep elements are used to enforce weak Dirichlet, Neumann, and coupling boundary
conditions. A corresponding element formulation is introduced. Appropriate solutions for local refinement and the
application of arbitrarily located loads within IBRA are presented as well.
The different benchmark examples clearly demonstrate the power of IBRA. IBRA is flexible for trimming toler-
ances, robust to parametric distortion, and can handle non-watertight models. These features have been demonstrated
by the benchmark examples, which cover a gamut of problems ranging from linear to nonlinear shell problems, dis-
cretized with trimmed and non-matching patches. The main features of the AiCAD workflow resp. IBRA are as
follows:
– usage of good analysis properties of NURBS basis functions, as demonstrated already by IGA
– unification of the design and analysis models because IBRA uses an uniform geometry representation for design
and analysis
- geometry transformations with associated possible drawbacks, e.g., loss of information are not necessary
(no meshing)
- accelerated design cycle
– the IBRA solution is provided as a CAD representation, which allows for high quality postprocessing (high
continuity, evaluation of stress at arbitrary points,...)
– current shortcomings of IGA with respect to a full design-through-analysis workflow, i.e., the creation of a
direct and complete analysis model from CAD, are now overcome in a general and consistent manner
A prototype of the AiCAD workflow was developed within the CAD programs Rhinoceros 4, 5
R
and Siemens NX

R
8.5 in combination with the in-house developed finite-element code CARAT++ and the interface software TEDA.
From the experience gained by developing the IBRA concept, we are convinced that it provides an analysis frame-
work that is general and robust enough for application to very large and complex industrial problems (see also Figure
61). Although the general concept has been established successfully, further research is necessary. Some interesting
fields are as follows:
– more efficient variants of the nested Jacobian approach (NEJA) for integrating the trimmed domain
– development of alternative B-Rep element formulations for shell coupling for example without a penalty factor
that must be chosen manually
– development of B-Rep element formulations for other applications (contact, optimization, FSI,...)
– development of efficient integration and displacement mapping techniques for IBRA solids
With the present concept and future developments in the mentioned fields of research, consistent integration of
design and analysis in industry would be closer to realization.
9. Acknowledgment
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) as part of the SPP project ”Leicht
Bauen mit Beton”. The support is gratefully acknowledged.
61
References
[1] J. Parvizian, A. Düster, E. Rank, Finite cell method (2007). doi:10.1007/s00466-007-0173-y.
[2] A. Düster, J. Parvizian, Z. Yang, E. Rank, The finite cell method for three-dimensional problems of solid mechanics, Computer Methods in
Applied Mechanics and Engineering 197 (45-48) (2008) 3768–3782. doi:10.1016/j.cma.2008.02.036.
[3] D. Schillinger, L. Dedè, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, T. J. 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 249-252 (2012) 116–150. doi:10.1016/j.cma.2012.03.017.
[4] D. Schillinger, M. Ruess, N. Zander, Y. Bazilevs, A. Düster, E. Rank, Small and large deformation analysis with the p- and B-spline versions
of the Finite Cell Method (2012). doi:10.1007/s00466-012-0684-z.
[5] E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, A. Düster, Geometric modeling, isogeometric analysis and the finite cell method,
Computer Methods in Applied Mechanics and Engineering 249-252 (2012) 104–115. doi:10.1016/j.cma.2012.05.022.
[6] L. Kantorovich, V. Krylov, Approximate methods of higher analysis, Interscience Publishers, 1958.
[7] Scan&Solve, http://www.intact-solutions.com.
[8] T. Hughes, J. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer
Methods in Applied Mechanics and Engineering 194 (39-41) (2005) 4135–4195. doi:10.1016/j.cma.2004.10.008.
[9] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs, Isogeometric analysis: Toward integration of CAD and FEA, Wiley, Chichester and West Sussex
and U.K and Hoboken and NJ, 2009.
[10] L. A. Piegl, W. Tiller, The NURBS book, 2nd Edition, Springer, Berlin and New York, op. 1997.
[11] D. F. Rogers, An introduction to NURBS: With historical perspective, Morgan Kaufmann Publishers, San Francisco, 2001.
[12] J. Cottrell, A. Reali, Y. Bazilevs, T. Hughes, Isogeometric analysis of structural vibrations, Computer Methods in Applied Mechanics and
Engineering 195 (41-43) (2006) 5257–5296. doi:10.1016/j.cma.2005.09.027.
[13] F. Auricchio, L. Beirão da Veiga, T. Hughes, A. Reali, G. Sangalli, Isogeometric collocation for elastostatics and explicit dynamics, Computer
Methods in Applied Mechanics and Engineering 249-252 (2012) 2–14. doi:10.1016/j.cma.2012.03.026.
[14] R. Schmidt, R. Wüchner, K.-U. Bletzinger, Isogeometric analysis of trimmed NURBS geometries, Computer Methods in Applied Mechanics
and Engineering 241-244 (2012) 93–111. doi:10.1016/j.cma.2012.05.021.
[15] T. Belytschko, R. W. Lewis, D. Perić, R. L. Taylor, Isogeometric analysis of nearly incompressible solids, International Journal for Numerical
Methods in Engineering 87 (1-5) (2011) 273–288. doi:10.1002/nme.3048.
[16] N. Nguyen-Thanh, J. Kiendl, H. Nguyen-Xuan, R. Wüchner, K. Bletzinger, Y. Bazilevs, T. Rabczuk, Rotation free isogeometric
thin shell analysis using PHT-splines, Computer Methods in Applied Mechanics and Engineering 200 (47-48) (2011) 3410–3424.
doi:10.1016/j.cma.2011.08.014.
[17] D. Benson, S. Hartmann, Y. Bazilevs, M.-C. Hsu, T. Hughes, Blended isogeometric shells, Computer Methods in Applied Mechanics and
Engineering 255 (2013) 133–146. doi:10.1016/j.cma.2012.11.020.
[18] J. Cottrell, T. Hughes, A. Reali, Studies of refinement and continuity in isogeometric structural analysis, Computer Methods in Applied
Mechanics and Engineering 196 (41-44) (2007) 4160–4183. doi:10.1016/j.cma.2007.04.007.
[19] Y. Bazilevs, V. Calo, J. Cottrell, T. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for
large eddy simulation of incompressible flows, Computer Methods in Applied Mechanics and Engineering 197 (1-4) (2007) 173–201.
doi:10.1016/j.cma.2007.07.016.
[20] Y. Zhang, Y. Bazilevs, S. Goswami, C. L. Bajaj, T. J. Hughes, Patient-specific vascular NURBS modeling for isogeometric analysis of blood
flow (2007). doi:10.1016/j.cma.2007.02.009.
[21] Y. Bazilevs, V. M. Calo, Y. Zhang, T. J. R. Hughes, Isogeometric Fluid–structure Interaction Analysis with Applications to Arterial Blood
Flow (2006). doi:10.1007/s00466-006-0084-3.
[22] M.-C. Hsu, I. Akkerman, Y. Bazilevs, High-performance computing of wind turbine aerodynamics using isogeometric analysis, Computers
& Fluids 49 (1) (2011) 93–100. doi:10.1016/j.compfluid.2011.05.002.
[23] Y. Bazilevs, M.-C. Hsu, M. Scott, Isogeometric fluid–structure interaction analysis with emphasis on non-matching discretiza-
tions, and with application to wind turbines, Computer Methods in Applied Mechanics and Engineering 249-252 (2012) 28–41.
doi:10.1016/j.cma.2012.03.028.
[24] İ. Temizer, P. Wriggers, T. Hughes, Three-dimensional mortar-based frictional contact treatment in isogeometric analysis with NURBS
(2012). doi:10.1016/j.cma.2011.10.014.
[25] R. Dimitri, L. d. Lorenzis, M. Scott, P. Wriggers, R. Taylor, G. Zavarise, Isogeometric large deformation frictionless contact using T-splines,
Computer Methods in Applied Mechanics and Engineering 269 (2014) 394–414. doi:10.1016/j.cma.2013.11.002.
[26] L. d. Lorenzis, İ. Temizer, P. Wriggers, G. Zavarise, A large deformation frictional contact formulation using NURBS-based isogeometric
analysis, International Journal for Numerical Methods in Engineering (2011) n/adoi:10.1002/nme.3159.
[27] W. A. Wall, M. A. Frenzel, C. Cyron, Isogeometric structural shape optimization (2008). doi:10.1016/j.cma.2008.01.025.
[28] J. Kiendl, R. Schmidt, R. Wüchner, K.-U. Bletzinger, Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and
sensitivity weighting, Computer Methods in Applied Mechanics and Engineering 274 (2014) 148–167. doi:10.1016/j.cma.2014.02.001.
[29] X. Qian, O. Sigmund, Isogeometric shape optimization of photonic crystals via Coons patches, Computer Methods in Applied Mechanics
and Engineering 200 (25-28) (2011) 2237–2255. doi:10.1016/j.cma.2011.03.007.
[30] T. Sederberg, J. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCs, ACM Trans. Graphics 22 22 (477-484).
[31] T. Sederberg, D. L. Cardon, J. Zheng, Lyche T., T-spline Simplication and Local Refinement, ACM Trans. Graphics 22 23 (276-283).
[32] Y. Bazilevs, V. Calo, J. Cottrell, J. Evans, T. Hughes, S. Lipton, M. Scott, T. Sederberg, Isogeometric analysis using T-splines (2010).
doi:10.1016/j.cma.2009.02.036.
[33] D. Burkhart, Hamann B., Umlauf G., Iso-geometric Finite Element Analysis Based on Catmull-Clark Subdivision Solids, Comput. Graph.
Forum 29 (5) (2010) 1575–1784.
62
[34] F. Cirak, M. Ortiz, P. Schröder, Subdivision surfaces: a new paradigm for thin-shell finite-element analysis, International Journal for Numer-
ical Methods in Engineering 47 (2000) 2039–2072.
[35] F. Cirak, Q. Long, Subdivision shells with exact boundary control and non-manifold geometry (2011). doi:10.1002/nme.3206.
[36] M. Mäntylä, An introduction to solid modeling, Vol. 13 of Principles of computer science series, Computer Science Press, Rockville and Md,
1988.
[37] M. E. Mortenson, Geometric modeling, Wiley, New York, 1985.
[38] H.-J. Kim, Y.-D. Seo, S.-K. Youn, Isogeometric analysis for trimmed CAD surfaces, Computer Methods in Applied Mechanics and Engi-
neering 198 (37-40) (2009) 2982–2995. doi:10.1016/j.cma.2009.05.004.
[39] Y.-W. Wang, Z.-D. Huang, Y. Zheng, S.-G. Zhang, Isogeometric analysis for compound B-spline surfaces, Computer Methods in Applied
Mechanics and Engineering 261-262 (2013) 1–15. doi:10.1016/j.cma.2013.04.001.
[40] Rhinoceros, http://www.rhino3d.com.
[41] Siemens NX, http://www.plm.automation.siemens.com.
[42] J. Rooney, P. Steadman, Computer-aided Design, Open University Press, Milton Keynes, 1987.
[43] I. Stroud, Boundary representation modelling techniques, Springer, London, 2006.
[44] K. Reed, D. Harrod, W. Conroy, The Initial Graphics Exchange Specification (IGES) version 5.0, Vol. 4412 of NISTIR, National Institute of
Standards and Technology, Gaithersburg and MD, 1990.
[45] P. M. Knupp, S. Steinberg, The fundamentals of grid generation, CRC Press, Boca Raton [u.a.], 1994.
[46] B. H. V. Topping, Finite element mesh generation, Saxe-Coburg publications on computational engineering, Saxe-Coburg Publications,
Stirling and Scotland, 2004.
[47] G. A. Hansen, R. W. Douglass, A. Zardecki, Mesh enhancement: Selected elliptic methods, foundations and applications, Imperial College
Press, London, 2005.
[48] T. W. Sederberg, G. T. Finnigan, X. Li, H. Lin, H. Ipson, Watertight trimmed NURBS, ACM Transactions on Graphics 27 (3) (2008) 1.
doi:10.1145/1360612.1360678.
[49] Abel Coll Sans, Robust volume mesh generation for non-watertight geometries, Ph.D. thesis, CIMNE, Barcelona (May 2014).
[50] B. K. Choi, Surface modeling for CAD/CAM, Vol. 11 of Advances in industrial engineering, Elsevier, Amsterdam and New York, 1991.
[51] E. Cohen, R. F. Riesenfeld, G. Elber, Geometric modeling with splines: An introduction, AK Peters, Natick and Mass, 2001.
[52] C. d. Boor, A practical guide to splines: With 32 figures, rev. ed. Edition, Vol. 27 of Applied mathematical sciences, Springer, New York,
2001.
[53] L. Piegl, W. Tiller, Geometry-based triangulation of trimmed NURBS surfaces, Computer Aided Design 1 (30) (1998) 11–18.
[54] G. Renner, V. Weiß, Exact and approximate computation of B-spline curves on surfaces, Computer-Aided Design 36 (4) (2004) 351–362.
doi:10.1016/S0010-4485(03)00100-3.
[55] N. M. Patrikalakis, Surface-to-surface intersections, IEEE Computer Graphics and Applications 13 (1) (1993) 89–95.
[56] T. W. Sederberg, R. J. Meyers, Loop detection in surface patch intersections, Computer Aided Geometric Design 5 (2) (1988) 161–171.
[57] S. Krishnan, D. Manocha, An efficient surface intersection algorithm based on lower-dimensional formulation, ACM Transactions on Graph-
ics (TOG) 16 (1) (1997) 74–106.
[58] Y. L. Ma, W. Hewitt, Point inversion and projection for NURBS curve and surface: Control polygon approach, Computer Aided Geometric
Design 20 (2) (2003) 79–99. doi:10.1016/S0167-8396(03)00021-9.
[59] P. Wriggers, Nichtlineare Finite-Element-Methoden, Springer, Berlin [etc.], 2001.
[60] W. Wunderlich, W. D. Pilkey, Mechanics of structures: Variational and computational methods, 2nd Edition, CRC Press, Boca Raton and FL,
2003.
[61] J. Kiendl, K.-U. Bletzinger, J. Linhard, R. Wüchner, Isogeometric shell analysis with Kirchhoff–Love elements, Computer Methods in
Applied Mechanics and Engineering 198 (49-52) (2009) 3902–3914. doi:10.1016/j.cma.2009.08.013.
[62] D. Benson, Y. Bazilevs, M. Hsu, T. Hughes, Isogeometric shell analysis: The Reissner–Mindlin shell, Computer Methods in Applied Me-
chanics and Engineering 199 (5-8) (2010) 276–289. doi:10.1016/j.cma.2009.05.011.
[63] W. Dornisch, S. Klinkel, B. Simeon, Isogeometric Reissner–Mindlin shell analysis with exactly calculated director vectors, Computer Meth-
ods in Applied Mechanics and Engineering 253 (2013) 491–504. doi:10.1016/j.cma.2012.09.010.
[64] R. Echter, B. Oesterle, M. Bischoff, A hierarchic family of isogeometric shell finite elements (2013). doi:10.1016/j.cma.2012.10.018.
[65] J. M. Kiendl, Isogeometric analysis and shape optimal design of shell structures, Shaker, Aachen, 2011.
[66] A. P. Nagy, D. J. Benson, On the numerical integration of trimmed isogeometric elements, Computer Methods in Applied Mechanics and
Engineeringdoi:10.1016/j.cma.2014.08.002.
[67] Y. Basar, W. Krätzig, Mechanik der Flächentragwerke. Theorie, Berechnungsmethoden, Anwendungsbeispiele, Vieweg, Braunschweig, 1985.
[68] A.-V. Vuong, C. Giannelli, B. Jüttler, B. Simeon, A hierarchical approach to adaptive local refinement in isogeometric analysis, Computer
Methods in Applied Mechanics and Engineering 200 (49-52) (2011) 3554–3567. doi:10.1016/j.cma.2011.09.004.
[69] D.Schillinger, The p- and B-spline versions of the geometrically nonlinear finite cell method and hierarchical refinement strategies for adaptive
isogeometric and embedded domain analysis, Ph.D. thesis, Technischen Universität München, München (22.May 2012).
[70] M. Scott, X. Li, T. Sederberg, T. Hughes, Local refinement of analysis-suitable T-splines, Computer Methods in Applied Mechanics and
Engineering 213-216 (2012) 206–222. doi:10.1016/j.cma.2011.11.022.
[71] G. Kuru, C. Verhoosel, K. van der Zee, E. van Brummelen, Goal-adaptive Isogeometric Analysis with hierarchical splines, Computer Methods
in Applied Mechanics and Engineering 270 (2014) 270–292. doi:10.1016/j.cma.2013.11.026.
[72] M. Ruess, D. Schillinger, A. I. Özcan, E. Rank, Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geome-
tries, Computer Methods in Applied Mechanics and Engineering 269 (2014) 46–71. doi:10.1016/j.cma.2013.10.009.
[73] M. Ruess, D. Schillinger, Y. Bazilevs, V. Varduhn, E. Rank, Weakly enforced essential boundary conditions for NURBS-embedded and
trimmed NURBS geometries on the basis of the finite cell method, International Journal for Numerical Methods in Engineering 95 (10)
(2013) 811–846. doi:10.1002/nme.4522.
[74] ANSYS, http://www.ansys.com.
63
[75] P. L. Gould, Introduction to linear elasticity, 2nd Edition, Springer-Verlag, New York, 1994.
[76] J. Kiendl, Y. Bazilevs, M.-C. Hsu, R. Wüchner, 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) (2010) 2403–2416.
doi:10.1016/j.cma.2010.03.029.
[77] T. Hughes, A. Reali, G. Sangalli, Efficient quadrature for NURBS-based isogeometric analysis, Computer Methods in Applied Mechanics
and Engineering 199 (5-8) (2010) 301–313. doi:10.1016/j.cma.2008.12.004.
[78] K. Sze, X. Liu, S. Lo, Popular benchmark problems for geometric nonlinear analysis of shells, Finite Elements in Analysis and Design 40 (11)
(2004) 1551–1569. doi:10.1016/j.finel.2003.11.001.
64

Breitenberger IBRA

Uploaded by

Copyright:

Available Formats

Breitenberger IBRA

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Breitenberger IBRA

Uploaded by

Copyright:

Available Formats

Analysis in Computer Aided Design:

Nonlinear Isogeometric B-Rep Analysis of Shell Structures

M. Breitenbergera,∗, A. Apostolatosa , B. Philippa , R. Wüchnera , K.-U. Bletzingera

2. Design and Analysis Modeling

2.1. Design Modeling

2.2. Analysis Modeling – Standard Analysis (CAD-CAE) Workflow

2.3. Analysis Modeling – Analysis in Computer Aided Design (AiCAD) Workflow

3. Design – NURBS based B-Rep Model

CAD geometry kernel Definition of analysis- CAD geometry kernel

CAD geometry kernel

analysis modeling (preprocessing)

finite element mesh refined design model

CAD geometry kernel

standard analysis model AiCAD analysis model

classical finite element analysis isogeometric B-Rep analysis

approximated solution represented approximated solution represented

b) CSG model: composition of simpler solids

c) NURBS based B-Rep model: composition of

Figure 3: Different geometry representations of same object (after [50])

3.2. Non-Uniform Rational B-Splines (NURBS)

3.2.1. Knot Vector

3.2.2. Basis Functions

3.2.3. NURBS Curves

3.2.4. NURBS Surfaces

and NURBS surfaces are defined by

3.2.5. Geometry Refinement

3.3. Trimmed NURBS Surfaces

Ck (ξ̃) = S(uk (ξ̃), vk (ξ̃)), k = 1, 2, ..., M (8)

3.4.1. Trimming Tolerances

a) Surface B-Rep model

b) Trimmed B-Spline surface 1 c) Trimmed NURBS surface 2 d) Trimmed NURBS surface 3

Figure 9: NURBS-based surface B-Rep model

a) Original B-Spline surface 1 b) Original NURBS surfaces 2 and 3

c) Intersecting surfaces before trimming operations

d) Intersecting surfaces after trimming operations

Figure 10: Boolean operations including trimming of intersecting surfaces

3.5.2. Topology of Surface B-Rep Models

4.1. Analysis of Shells

4.1.1. Geometry Description

X(θ1 , θ2 , θ3 ) = Xsurf (θ1 , θ2 ) + θ3 A3 , x(θ1 , θ2 , θ3 ) = xsurf (θ1 , θ2 ) + θ3 a3 , (9)

A1 and A2 resp. a1 and a2 are defined by

u = usurf = xsurf − Xsurf . (12)

δW = δWint + δWext = 0. (13)

4.2. Isogeometric Analysis

4.2.1. NURBS-based Kirchhoff-Love (KL) Shell Element

4.2.2. Numerical Integration of NURBS Surfaces

w̃l = J3 (ξl , ηl ) wl (25)

5. Isogeometric B-Rep Analysis

5.1. Isogeometric B-Rep Elements

Ω(1) Ω(2) Γ(2)

– Coupling boundary conditions

∂Wint ∂Wext X δWB-Rep

5.2. NURBS-based Isogeometric B-Rep Analysis of Shells

∂Xsurf ∂θ1 ∂θ2

∂θ1 ∂θ2 ∂θ1 ∂θ2

5.3. Numerical Integration of Trimmed Surface Domains

The knot spans (elements) can be categorized as follows:

5.3.1. Nested Jacobian Approach

Z (1) (1) (1)