# CAD-compatible structural shape optimization with a movable Bézier
tetrahedral mesh

Jorge López  Cosmin Anitescu  Timon Rabczuk<EMAIL_ADDRESS>Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh
City, Viet Nam  Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi
Minh City, Viet Nam  Institute of Structural Mechanics, Bauhaus-Universität
Weimar, Marienstr. 15, D-99423 Weimar, Germany

###### Abstract

This paper presents the development of a complete CAD-compatible framework for
structural shape optimization in 3D. The boundaries of the domain are
described using NURBS while the interior is discretized with Bézier
tetrahedra. The tetrahedral mesh is obtained from the mesh generator software
Gmsh. A methodology to reconstruct the NURBS surfaces from the triangular
faces of the boundary mesh is presented. The description of the boundary is
used for the computation of the analytical sensitivities with respect to the
control points employed in surface design. Further, the mesh is updated at
each iteration of the structural optimization process by a pseudo-elastic
moving mesh method. In this procedure, the existing mesh is deformed to match
the updated surface and therefore reduces the need for remeshing. Numerical
examples are presented to test the performance of the proposed method. The use
of the movable mesh technique results in a considerable decrease in the
computational effort for the numerical examples.

###### keywords:

Shape optimization , Bézier tetrahedra , Isogeometric analysis , Moving mesh ,
IGES

##  1 Introduction

Modeling and analysis have typically been distinguished as two different yet
complementary phases for the engineering practice. However, for the analysis
using the conventional numerical methods, a discrete mesh is essential for the
geometries that have been modeled using Computer-Aided Design (CAD). Various
methods exist in the literature [44, 5] to discretize the geometry and the
solution fields. In this regard, Isogeometric Analysis (IGA) [10] emerged as
an alternative to Finite Element Method (FEM), proposing the merging of the
analysis and design domains. In this framework, the same basis functions which
are employed in CAD to describe the geometry are also used for the
approximation of the solution fields at the analysis stage.

In the context of shape optimization, IGA provides the opportunity to
represent complex geometries with few design variables [37, 27, 21, 24, 25,
41]. This is in contrast to the node-based optimization approach using FEM,
where element vertices are used as design variables. Moreover, B-splines and
NURBS basis functions, commonly used in CAD software, offer a variety of
refinement algorithms that keep the geometry intact.

Even though IGA has several advantages, some difficulties have been identified
particularly in case of three-dimensional domains. First, the CAD models might
parameterize only the boundaries and not the interior of the domain. A CAD
model frequently uses a collection of connected surfaces to represent the
geometry, omitting the information required to parametrize the volume. Also,
complex geometries cannot be represented with tensor product basis functions.
Often, a designer will model geometries that are not analysis-suitable.
Furthermore, compatibility between the different CAD software packages and the
analysis platforms is crucial. Considering this, in the context of shape
optimization, a CAD-compatible model must meet the following requirements:

  * •

It should employ commonly used CAD basis functions to parameterize the
boundary of the geometry.

  * •

It must be capable of modelling complex geometries.

  * •

Exchange of information between the geometry and the analysis software should
be possible.

  * •

It should be able to maintain a link between the design variables and the
analysis mesh.

To overcome the issues discussed above, we broadly present the features of the
proposed approach in the context of CAD-compatible shape optimization. First,
we propose to use the IGES exchange format to communicate between the design
and the analysis mesh. Second, our approach uses Bézier elements for the
discretizetion of the domain, where the geometry can be exactly represented in
a CAD-compatible format. This might appear an extension of the NURBS Enhanced
Finite Element Method (NEFEM) [29, 30, 31] with the use of IGA. However, in
NEFEM the parent element is modified for the curved elements at the
boundaries, and therefore the numerical quadrature is altered for the
integration of those elements, which is in contrast to the method proposed in
this paper. In this work, we provide a linkage between CAD-based design and
optimization analysis in a more flexible way than what is provided by standard
IGA. The NURBS-based IGA still has some difficulties with volumetric mesh
generation since the problem is essentially equivalent to that of hexahedral
meshing. Moreover, multiple patches or parameterization singularities are
still required for classic IGA. In this regard, the Bézier tetrahedra can
function as a bridge between surface CAD models and volumetric
discretizations.

To make the proposed method efficient for 3D models, we develop a method based
on the point inversion algorithm [26] to reconstruct the NURBS boundaries.
This approach is particularly important in three-dimensional domains, where
converting a NURBS surface into triangular elements requires that the degree
of the triangle elements is at least twice of the degree of the tensor product
description of the initial surface [9], thereby creating a large number of
control points. By contrast, the link between the boundary representation and
the analysis mesh is more easily achieved for two-dimensional problems by
using $C^{0}$ continuous curve segments as one side of individual triangular
elements [7, 27, 19].

Lastly, the proposed approach exploits the advantages of linking the NURBS
boundary with the Bézier tetrahedral mesh. The sensitivities can be computed
analytically in the analysis mesh with respect to the coordinates of the NURBS
control points. Moreover, it provides the possibility to update the mesh
without remeshing at each iteration, avoiding changes in the connectivity of
the elements and resulting in faster convergence rates without oscillations.

The use of simplex elements in an isogeometric framework has been of interest
in recent years. Similar work that has already been carried out include the
implementation of Bézier triangles in the framework of IGA for applications in
vibrations [16], thin plates and shells [20, 42], damage mechanics[17] and
phase separation [43]. The increase of continuity to $C^{1}$ Bézier triangles
has been studied in [11, 40]. Further extensions of 3D implementations using
Bézier tetrahedral meshes have been developed in [39, 13], and for more
general Bézier elements in [8]. In the context of shape optimization, Bézier
triangles have been implemented for the design and optimization of 2D auxetic
materials in [38], and for stress constraints in [19]. Also, the
implementation of the projection procedures to approximate the Bézier elements
faces to the NURBS surface has been shown in [8]. Lastly, some work on moving
mesh techniques have been commonly presented for fluid-structure interaction
[1, 2, 32, 33, 35]. In the context of optimization, the moving mesh approach
has been used in [4, 18, 38], where the update of the mesh by an elastic body
model has been preferred.

The remainder of the paper is organized as follows. In Section 2, we introduce
the basics of Computer Aided Graphic Design (CAGD). We discuss the basics of
elements and basis functions typically used in this modeling approach. In
Section 3, an overview of the construction of volumetric objects using the
IGES exchange format is given. In Section 4, we describe the approach to
reconstruct the NURBS boundaries and the method for smoothing the weights in
the interior of the domain. In Section 5, the implementation of the proposed
approach for shape optimization is introduced and further explained. The
analytical sensitivity analysis and the method to update the mesh based on a
pseudo-elastic mesh is also presented. Numerical examples illustrating the
performance of the proposed approach are presented in Section 6, followed by
concluding remarks and possibilities for future work.

##  2 Fundamentals of CAGD

###  2.1 Tetrahedral elements

Unstructured meshes, typically tetrahedral elements for three-dimensional
modeling, do not present a regular pattern on the topology or connectivity of
their elements. For this reason, the storage of the connectivity between
elements is required, which implies higher computational cost in comparison to
structured meshes. However, their flexibility to model complex geometries
gives them an advantage over their structured counterparts. In this section,
we introduce the basic concepts of tetrahedral elements typically used in CAD.

####  2.1.1 Barycentric coordinate system

First, we describe the barycentric coordinates for a tetrahedron, which later
will be used to define basis functions over the element. This system, first
introduced by Möbius [23], determines the location of a point in
$\mathbb{R}^{3}$ as the barycenter of a tetrahedron. Namely, the value of
masses located at the 4 vertices of the tetrahedron are given, such that the
center of mass of the element can be located. Keeping this interpretation in
mind, the barycentric coordinates
$\bm{\lambda}=[\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}]$ can be
written in terms of the volume of a tetrahedron as:

$\lambda_{i}=\frac{V_{i}}{V},$ (1)

where $V$ denotes the volume of the element and $V_{i}$ is the volume of the
tetrahedron formed by the point of interest and the three opposite vertices to
the vertex $i$, as shown in Fig. 1. Let us consider a tetrahedron $\bm{T}$ in
the Euclidean space $\mathbb{R}^{3}$ with vertices at $\bm{P}_{1}$,
$\bm{P}_{2}$, $\bm{P}_{3}$ and $\bm{P}_{4}$. Using the barycentric coordinate
system, any point $\bm{P}\in\mathbb{R}^{3}$ can be defined as a linear
combination of the vertices of the tetrahedron with the barycentric
coordinates $\bm{\lambda}=[\lambda_{1},$ $\lambda_{2},$ $\lambda_{3},$
$\lambda_{4}]$:

$\bm{P}(\bm{\lambda})=\lambda_{1}\bm{P}_{1}+\lambda_{2}\bm{P}_{2}+\lambda_{3}\bm{P}_{3}+\lambda_{4}\bm{P}_{4},$
(2)

satisfying:

$\lambda_{1}+\lambda_{2}+\lambda_{3}+\lambda_{4}=1.$ (3)

Figure 1: Barycentric coordinate system in a tetrahedron.

####  2.1.2 Bézier tetrahedra

Bézier curves are parametric curves defined as a linear combination of control
points and Bernstein polynomials. Bernstein polynomials can be defined over
tetrahedral domains by making use of the previously introduced barycentric
coordinates $\bm{\lambda}$. We consider the parent element in the space
$(\xi,\eta,\zeta)$ with vertices as shown in Fig. 2. Some of the properties
exhibited by the Bézier basis are partition of unity, pointwise nonnegativity
and endpoint interpolation [14, 8]. The Bézier-Bernstein basis functions of
degree $p$ over the parent tetrahedron are defined as:

Figure 2: Parent element and its mapping into the physical space.

$B_{ijkl}^{p}(\bm{\lambda})=\frac{p!}{i!j!k!l!}\lambda_{1}^{i}\lambda_{2}^{j}\lambda_{3}^{k}\lambda_{4}^{l},$
(4)

and its derivatives as:

$\begin{array}[]{r@{}l}\frac{d}{d\lambda_{1}}B_{ijkl}^{p}(\bm{\lambda})&=\frac{p!}{(i-1)!j!k!l!}\lambda_{1}^{i-1}\lambda_{2}^{j}\lambda_{3}^{k}\lambda_{4}^{l},\vspace{5pt}\\\
\frac{d}{d\lambda_{2}}B_{ijkl}^{p}(\bm{\lambda})&=\frac{p!}{i!(j-1)!k!l!}\lambda_{1}^{i}\lambda_{2}^{j-1}\lambda_{3}^{k}\lambda_{4}^{l},\vspace{5pt}\\\
\frac{d}{d\lambda_{3}}B_{ijkl}^{p}(\bm{\lambda})&=\frac{p!}{i!j!(k-1)!l!}\lambda_{1}^{i}\lambda_{2}^{j}\lambda_{3}^{k-1}\lambda_{4}^{l},\vspace{5pt}\\\
\frac{d}{d\lambda_{4}}B_{ijkl}^{p}(\bm{\lambda})&=\frac{p!}{i!j!k!(l-1)!}\lambda_{1}^{i}\lambda_{2}^{j}\lambda_{3}^{k}\lambda_{4}^{l-1},\end{array}$
(5)

where $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$, $\lambda_{4}$ are the
barycentric coordinates on the reference element; and the indexes $i+j+k+l=p$.
Usually, the derivatives of the basis functions are computed with respect to
the parametric coordinates as:

$\begin{array}[]{r@{}l}\frac{d}{d\xi}B_{ijkl}^{p}(\bm{\lambda})&=\frac{\partial\lambda_{1}}{\partial\xi}\frac{d}{d\lambda_{1}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{2}}{\partial\xi}\frac{d}{d\lambda_{2}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{3}}{\partial\xi}\frac{d}{d\lambda_{3}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{4}}{\partial\xi}\frac{d}{d\lambda_{4}}B_{ijkl}^{p}(\bm{\lambda}),\vspace{8pt}\\\
\frac{d}{d\eta}B_{ijkl}^{p}(\bm{\lambda})&=\frac{\partial\lambda_{1}}{\partial\eta}\frac{d}{d\lambda_{1}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{2}}{\partial\eta}\frac{d}{d\lambda_{2}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{3}}{\partial\eta}\frac{d}{d\lambda_{3}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{4}}{\partial\eta}\frac{d}{d\lambda_{4}}B_{ijkl}^{p}(\bm{\lambda}),\vspace{8pt}\\\
\frac{d}{d\zeta}B_{ijkl}^{p}(\bm{\lambda})&=\frac{\partial\lambda_{1}}{\partial\zeta}\frac{d}{d\lambda_{1}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{2}}{\partial\zeta}\frac{d}{d\lambda_{2}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{3}}{\partial\zeta}\frac{d}{d\lambda_{3}}B_{ijkl}^{p}(\bm{\lambda})+\frac{\partial\lambda_{4}}{\partial\zeta}\frac{d}{d\lambda_{4}}B_{ijkl}^{p}(\bm{\lambda}),\\\
\end{array}$ (6)

For the parent tetrahedra from Fig. 2 we use the following relations:

$\begin{split}\lambda_{1}&=-\xi-\eta-\zeta+1,\\\
\lambda_{2}&=\xi,\quad\lambda_{3}=\eta,\quad\lambda_{4}=\zeta.\end{split}$ (7)

Using the vertex coordinates from the parent element, Eq. (6) can be written
as:

$\begin{array}[]{r@{}l}\frac{d}{d\xi}B_{ijkl}^{p}(\bm{\lambda})&=-\frac{d}{d\lambda_{1}}B_{ijkl}^{p}(\bm{\lambda})+\frac{d}{d\lambda_{2}}B_{ijkl}^{p}(\bm{\lambda}),\vspace{8pt}\\\
\frac{d}{d\eta}B_{ijkl}^{p}(\bm{\lambda})&=-\frac{d}{d\lambda_{1}}B_{ijkl}^{p}(\bm{\lambda})+\frac{d}{d\lambda_{3}}B_{ijkl}^{p}(\bm{\lambda}),\vspace{8pt}\\\
\frac{d}{d\zeta}B_{ijkl}^{p}(\bm{\lambda})&=-\frac{d}{d\lambda_{1}}B_{ijkl}^{p}(\bm{\lambda})+\frac{d}{d\lambda_{4}}B_{ijkl}^{p}(\bm{\lambda}).\\\
\end{array}$ (8)

Fig. 3 shows the convention used to order the basis functions for a cubic
tetrahedral parent element. Following this convention, a tetrahedron in the
physical space (see Fig. 2) is described as the linear combination of such
shape functions with the control points in the physical space $P_{ijkl}$:

$T(\bm{\lambda})=\sum_{i+j+k+l=p}B_{ijkl}^{p}(\bm{\lambda})P_{ijkl}.$ (9)

In the case that rational basis functions are used, the tetrahedron in the
physical space is defined by:

$T(\bm{\lambda})=\sum_{i+j+k+l=p}R_{ijkl}^{p}(\bm{\lambda})w_{ijkl}P_{ijkl},$
(10)

where $w_{ijkl}$ are the weights and the rational basis
$R_{ijkl}^{p}(\bm{\lambda})$ are computed by:

$R_{ijkl}^{p}(\bm{\lambda})=\frac{B_{ijkl}^{p}(\bm{\lambda})w_{ijkl}}{W(\bm{\lambda})},$
(11)

where $W(\bm{\lambda})$ is the weight function:

$W(\bm{\lambda})=\sum_{i+j+k+l=p}B_{ijkl}^{p}(\bm{\lambda})w_{ijkl}.$ (12)

####  2.1.3 Lagrange Tetrahedra

Though Bézier tetrahedra are commonly used in the CAD community, mesh
generation programs typically generate Lagrange elements, mainly because of
the interpolatory property at all nodes [28, 44]. Thus, the control points of
these elements can be located exactly over the boundaries of CAD models, which
usually use NURBS for the geometry description.

(a)

(b)

(c)

Figure 3: Node numbering of a a) quadratic, b) cubic and c) quartic Bézier
tetrahedron. Figure 4: Lagrange and Bézier control points required to model
exactly the same tetrahedron.

The Lagrange basis of degree $p$ over the same parent element as in the case
of Bézier tetrahedra, are defined as [44]:

$L_{ijkl}^{p}(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4})=\ell_{i}(\lambda_{1})\ell_{j}(\lambda_{2})\ell_{k}(\lambda_{3})\ell_{l}(\lambda_{4}),$
(13)

with $i+j+k+l=p$. The functions $\ell_{t}(\lambda_{s})$ are one-dimensional
Lagrange polynomials of degree $t$, defined by:

$\ell_{t}(\lambda_{s})=\prod_{\begin{subarray}{c}r=0\end{subarray}}^{t-1}\left(\frac{p\lambda_{s}-r}{t-r}\right).$
(14)

The same physical tetrahedron as in the last subsection $\bm{T}(\bm{\lambda})$
is in this case described as a linear combination of the Lagrange basis and
the corresponding control points:

$T(\bm{\lambda})=\sum_{i+j+k+l=p}L_{ijkl}^{p}(\bm{\lambda})P^{L}_{ijkl},$ (15)

where $P^{L}_{ijkl}$ are the control points of the Lagrange tetrahedron. Even
though the basis and control points are different, Bézier and Lagrange basis
can be used to model geometrically the same physical tetrahedron, provided
that all the weights are equal to one, as shown in Fig. 4.

####  2.1.4 Mesh size control

Figure 5: Refinement by splitting of a tetrahedron. One tetrahedron is split
into eight tetrahedra.

In order to control the size of the tetrahedral mesh, two different approaches
can be used together or separately. First the mesh size can be controlled by
setting a characteristic length $h_{l}$ at diferent points in the domain.In
this case, $h_{l}$ indicates the radius of a circumscribed sphere for a
tetrahedron at some location of the mesh. That is, the size of the mesh is
obtained by linearly interpolating the values of all given $h_{l}$. In this
work, a homogeoneous $h_{l}$ is given for all vertices of the boundary.

Additionally, the mesh can be further refined by splitting each element into
smaller tetrahedra. Fig. 5 shows the split of an element into eight elements
by cutting planes through the middle points of each edge of the tetrahedron.

###  2.2 NURBS

(a)

(b)

Figure 6: a) Cubic B-spline basis functions over the parameter space with knot
vector $\bm{\xi}=[0,0,0,0,\frac{1}{3},\frac{2}{3},1,1,1,1]$. b) B-Spline curve
with its control points, the red crosses show the knot span interfaces in the
physical space.

Let
$\bm{\Xi}=[\hat{\xi}_{1},\hat{\xi}_{2},...,\hat{\xi}_{n},...\hat{\xi}_{n+p+1}]$
be a vector of non-decreasing parameter values known as the knot vector, where
$n$ is the number of basis functions associated to $\bm{\hat{\xi}}$, and
$\hat{p}$ is the degree of those functions. A NURBS curve is a parametric
curve defined by:

$\bm{C}(\hat{\xi})=\frac{\sum_{i=1}^{n}N_{i,p}\hat{w}_{i}\hat{P}_{i}}{\sum_{j=1}^{n}N_{j,p}\hat{w}_{j}},$
(16)

where $\hat{w}_{i}$ is the $ith$ weight of the NURBS curve, $\hat{P}_{i}$ the
corresponding NURBS control point and $N_{i,p}$ are B-Spline basis functions
of degree $\hat{p}$. B-Splines are piecewise polynomials obtained recursively
by the Cox-de-Boor formula:

$N_{i,0}(\hat{\xi})=\begin{cases}1\hskip 5.0pt\mbox{if}\hskip
5.0pt\hat{\xi}_{i}\leq\hat{\xi}<\hat{\xi}_{i+1},\\\ 0\hskip
5.0pt\mbox{otherwise},\end{cases}$ (17)

$N_{i,p}(\hat{\xi})=\frac{\hat{\xi}-\hat{\xi}_{i}}{\hat{\xi}_{i+\hat{p}}-\hat{\xi}_{i}}N_{i,\hat{p}-1}(\hat{\xi})+\frac{\hat{\xi}_{i+\hat{p}+1}-\hat{\xi}}{\hat{\xi}_{i+\hat{p}+1}-\hat{\xi}_{i+1}}N_{i+1,\hat{p}-1}(\hat{\xi}).$
(18)

Fig. 6 shows the cubic B-Splines basis functions over the parameter space
$\hat{\xi}\in[0,1]$. A NURBS surface can be obtained by applying tensor
product of two NURBS curves with knot vectors
$\bm{\Xi}=[\hat{\xi}_{1},\hat{\xi}_{2},...,\hat{\xi}_{n},...\hat{\xi}_{n+p+1}]$
and
$\bm{H}=[\hat{\eta}_{1},\hat{\eta}_{2},...,\hat{\eta}_{n},...\hat{\eta}_{m+q+1}]$as:

$\bm{S}(\hat{\xi},\hat{\eta})=\sum_{i=1}^{n}\sum_{j=1}^{m}\hat{R}_{i,j}^{\hat{p},\hat{q}}(\hat{\xi},\hat{\eta})\hat{w}_{i,j}\hat{P}_{i,j},$
(19)

The basis $\hat{R}_{i,j}^{\hat{p},\hat{q}}$ are rational bivariate B-Splines
expressed by:

$\hat{R}_{i,j}^{\hat{p},\hat{q}}(\hat{\xi},\hat{\eta})=\frac{N_{i,\hat{p}}(\hat{\xi})M_{j,\hat{q}}\hat{\eta})\hat{w}_{i,j}}{\sum_{r=1}^{n}\sum_{s=1}^{m}N_{r,\hat{p}}(\hat{\xi})M_{s,\hat{q}}(\hat{\eta})\hat{w}_{r,s}},$
(20)

here, $N_{i,\hat{p}}$, $M_{j,\hat{q}}$ are the B-Splines basis functions of
degree $\hat{p}$, $\hat{q}$ respectively. Fig. 7 shows 2D B-Spline basis
function as a result of applying the tensor product of two 1D B-Spline basis
functions, and the physical surface.

(a)

(b)

Figure 7: a) Tensor product basis functions in 2D, cubic B-splines with knot
vector $\bm{\xi}=[0,0,0,0,\frac{1}{3},\frac{2}{3},1,1,1,1]$ are shown in the
$u$ direction, whereas quadratic B-splines over
$\bm{\eta}=[0,0,0,\frac{1}{2},1,1,1]$ are used in the $v$ direction. b)
Mapping into the physical space.

##  3 Volume representations with the IGES exchange format

Let us consider a 3D boundary representation of an object, from which a
tetrahedral mesh is generated to model the interior of the domain. If the
control points of the tetrahedra facing the boundary of the domain can
accurately represent the original boundary representation, we can use the
interior elements to solve the required PDEs for analysis while preserving the
initial geometry. In order to make the analysis CAD-compatible, the boundary
representation is generated by NURBS basis functions, which are mainly used
for CAD modeling and allow the creation of a broad variety of curves, i.e.
conic shapes. Bézier tetrahedra are used to mesh the interior of the domain.
They also can model conic surfaces, by means of rational basis, and in many
cases can give a reasonably good approximation of NURBS objects. An exchange
format, however, is required so the same model can be used by different
modeling, analysis and mesh generator software packages from a variety of
developers.

In the 1970s, efforts to solve this problems derived in the first version of
the IGES format. Later, using IGES as a foundation, the STEP file format was
published by the International Organization for Standarization (ISO) as a more
extensive standard. Together, IGES and STEP are the most used neutral file
formats used in the CAD industry. However, other formats developed by CAD
software companies are also used as exchange formats, for instance, STL is
nowadays broadly used in the 3D printing community.

In this work, we have used the IGES format, although other standards (such as
STEP) can be used analogously. IGES, Initial (initially Interim) Graphics
Exchange Specification, was the first standard created for the exchange of
models between different systems. IGES supports the translation of geometric,
topological and non-geometric information. Each basic design component is
included into the exchange format as an entity. An entity can contain
geometric or non-geometric information. Geometry entities contain information
to describe points, curves, surfaces, etc. and include also the relation or
connectivity with other geometric entities. Additionally, non-geometric
entities, mainly used for 2D drawings, include annotations or attributes such
as color.

The construction of NURBS volume domains is supported in IGES through the use
of the Manifold Boundary Solid Object (MBSO) entity. This entity defines a
volume by specifying the boundary of the 3D object. The boundary is first
defined by a set of shells, a connected boundary of a region in
$\mathbb{R}^{3}$. Shells in turn, are decomposed into faces. Each face must
point to a surface and a curve loop. The latter is provided in an
anticlockwise direction with respect to the outer face of the underlying
surface. Finally, the curves and vertices that bound the faces are listed.
Fig. 8 shows the construction of a volume domain using the MBSO entity with
the IGES exchange format.

Figure 8: Entity types used to define the MSBO.

##  4 Surface reconstruction

In the previous section, NURBS surfaces were used for a boundary
representation of a solid object. The next step is to model the interior of
the domain with a Bézier tetrahedral mesh. In this work, we used Gmsh to
obtain the finite element mesh, where the tetrahedral elements interpolate the
boundaries using Lagrange polynomial based elements. We assume that the
boundaries of the Lagrange elements approximate the original NURBS surfaces
with control points located exactly on the initial CAD description due to the
nodal interpolation property of the Lagrange polynomials. In order to use
Bézier elements instead, their control point coordinates are computed based on
the control net of the boundary NURBS surfaces.

A NURBS surface created by tensor product of curves of degree $p$ can be
exactly redefined with Bézier triangles of degree $2p$ [9], which might
greatly increase the number of degrees of freedom for the analysis phase.
Alternatively, a NURBS surface can be related to Bézier elements by using
Bézier extraction [6] and Bézier projection [36], as shown in [8]. However,
these are local operators and require each Bézier element to be located in
only one knot span of the NURBS surface. Here, we present a method to link the
elements provided by the mesh generator with the original boundary
description, regardless of whether the Bézier element faces are located in one
NURBS knot span or not.

First, the control points located at the boundaries of the Lagrange
tetrahedral elements are selected. Then, for each tetrahedron, the NURBS
parameter space coordinates $\hat{\xi}^{e}$, $\hat{\eta}^{e}$ are computed.
The mapping of these points into the physical space is thus defined as a
linear combination of the shape functions at the computed parameter points and
the NURBS control points. For each tetrahedra at the boundary, these physical
points should match with the mapping of the corresponding boundary face of the
parent Bézier element into the physical space. For one tetrahedron, this
relation is described as:

$\hat{\bm{x}}\mid_{\Omega^{e}}=\sum_{i=1}^{n}N_{i,p}(\hat{\xi}^{e},\hat{\eta}^{e})\hat{P}_{i}=\sum_{i=1}^{n_{f}}B_{i,p}(\hat{\lambda})P_{i},$
(21)

where $n$ is the number of basis functions of the corresponding NURBS surface
and $n_{f}$ is the number of control points on one face of a Bézier element
with corresponding barycentric coordinates $\bm{\hat{\lambda}}$ in the parent
element. In matrix form, Eq. (21) can be written as:

$\hat{\bm{x}}^{e}=\bm{N}^{e}\bm{\hat{P}}=\bm{B}(\hat{\lambda})\bm{P},$ (22)

with the matrices of basis functions defined as:

$\bm{N}^{e}=\begin{bmatrix}N_{1,p}(\hat{\xi}^{e}_{1},\hat{\eta}^{e}_{1})&N_{2,p}(\hat{\xi}^{e}_{1},\hat{\eta}^{e}_{1})&\cdots&N_{n,p}(\hat{\xi}^{e}_{1},\hat{\eta}^{e}_{1})\\\
N_{1,p}(\hat{\xi}^{e}_{2},\hat{\eta}^{e}_{2})&N_{2,p}(\hat{\xi}^{e}_{2},\hat{\eta}^{e}_{2})&\cdots&N_{n,p}(\hat{\xi}^{e}_{2},\hat{\eta}^{e}_{2})\\\
\vdots&\vdots&&\vdots\\\
N_{1,p}(\hat{\xi}^{e}_{n_{f}},\hat{\eta}^{e}_{n_{f}})&N_{2,p}(\hat{\xi}^{e}_{n_{f}},\hat{\eta}^{e}_{n_{f}})&\cdots&N_{n,p}(\hat{\xi}^{e}_{n_{f}},\hat{\eta}^{e}_{n_{f}})\end{bmatrix},$
(23)

$\bm{B}(\bm{\hat{\lambda}})=\begin{bmatrix}B_{1,p}(\hat{\lambda}_{1})&B_{2,p}(\hat{\lambda}_{1})&\cdots&B_{n_{f},p}(\hat{\lambda}_{1})\\\
B_{1,p}(\hat{\lambda}_{2})&B_{2,p}(\hat{\lambda}_{2})&\cdots&B_{n_{f},p}(\hat{\lambda}_{2})\\\
\vdots&\vdots&&\vdots\\\
B_{1,p}(\hat{\lambda}_{n_{f}})&B_{2,p}(\hat{\lambda}_{n_{f}})&\cdots&B_{n_{f},p}(\hat{\lambda}_{n_{f}})\end{bmatrix}.$
(24)

Thus, the location of the Bézier control points of the face at the boundary
are obtained by inverting the matrix containing the basis functions values at
the corresponding parameter points:

$\bm{P}=\bm{B}^{-1}(\bm{\hat{\lambda}})\bm{N}^{e}\bm{\hat{P}}.$ (25)

With this approach, it is also possible to recover the geometry in cases where
the weights are different than one. Defining the coordinates in the projected
space $\mathbb{R}^{4}$ by $\bm{\tilde{P}}=[xw,yw,zw,w]$, we obtain:

$\bm{\tilde{P}}=\bm{B}^{-1}(\bm{\hat{\lambda}})\bm{N}^{e}\bm{\hat{\tilde{P}}}.$
(26)

Figure 9: Reconstruction of a NURBS surface by computing the boundary control
points of Bézier tetrahedra.

The entire procedure to reconstruct the NURBS surface from the given Lagrange
triangular mesh is summarized graphically in Fig. 9.

###  4.1 Weight smoothing

The implementation of Eq. (26) results in control points with variable weights
over the surface of the domain. However, it has been proven that the variation
of the location of the control points should be minimized in order to achieve
optimal convergence [3, 22]. In the context of rational polynomials, it means
that also the variation of the weights has to be minimized. This is achieved
by solving the Laplace equation over the domain for the weights considering
the weights at the boundary as essential boundary conditions. The weak
formulation of the Laplace problem reads: Find
$w\in\mathcal{W}=\\{w|w\in(H^{1}(\Omega))^{3},w=\bar{w}\;\text{on }\Gamma\\}$
such that for all
$\mu\in\mathcal{W}_{0}=\\{\mu|\mu\in(H^{1}(\Omega))^{3},\mu=0\;\text{on
}\Gamma\\}$:

$a_{w}(w,\mu)=0,$ (27)

with the bilinear $a(\cdot)$ operator:

$a_{w}(w,\mu)=\int_{\Omega}\nabla w\cdot\nabla\mu\,d\Omega.$ (28)

Discretizing the continuous functions $w$ and $\mu$ with the same shape
functions for each tetrahedron, the discrete weak form reads:

$a_{w}(w^{h},\mu^{h})=0,$ (29)

where $w^{h}$ and $\mu^{h}$ are the discrete trial and test funstions.
Introducing the discretization into Eq. (28), we get the linear system of
equations:

$\bm{K}_{w}\bm{W}=\bm{F}_{w},$ (30)

where $\bm{F}_{w}$ is a fictitious force vector, $\bm{W}$ is the vector of
weights. The element stiffness matrix $\bm{K}^{e}_{w}$ for the Laplace problem
is defined as:

$\bm{K}^{e}_{w}=\int_{\Omega^{e}}\bm{B}^{T}_{w}\bm{I}_{3}\bm{B}_{w}\left|\bm{J}\right|d\Omega^{e},$
(31)

with $\bm{I}_{3}$ being the 3$\times$3 identity matrix and $\bm{J}$ the
Jacobian describing the mapping from the parent element to the physical space:

$\bm{J}=\begin{bmatrix}\frac{\partial x}{\partial\xi}&\frac{\partial
y}{\partial\xi}&\frac{\partial z}{\partial\xi}\vspace{5pt}\\\ \frac{\partial
x}{\partial\eta}&\frac{\partial y}{\partial\eta}&\frac{\partial
z}{\partial\eta}\vspace{5pt}\\\ \frac{\partial
x}{\partial\zeta}&\frac{\partial y}{\partial\zeta}&\frac{\partial
z}{\partial\zeta}\vspace{5pt}\end{bmatrix}.$ (32)

The matrix $\bm{B}_{w}$ is given by:

$\bm{B}_{w}=\begin{bmatrix}\frac{\partial R^{p}_{p,0,0,0}}{\partial
x}&\frac{\partial R^{p}_{0,p,0,0}}{\partial x}&\cdots&\frac{\partial
R^{p}_{i,j,k,l}}{\partial x}\\\ \frac{\partial R^{p}_{p,0,0,0}}{\partial
y}&\frac{\partial R^{p}_{0,p,0,0}}{\partial y}&\cdots&\frac{\partial
R^{p}_{i,j,k,l}}{\partial y}\\\ \frac{\partial R^{p}_{p,0,0,0}}{\partial
z}&\frac{\partial R^{p}_{0,p,0,0}}{\partial z}&\cdots&\frac{\partial
R^{p}_{i,j,k,l}}{\partial z}\end{bmatrix},\hskip 30.0pti+j+k+l=p.$ (33)

##  5 Shape optimization

Structural shape optimization is a process that seeks for improving the design
process by coupling geometrical modeling and analysis of a structure. Together
with the use of a proper optimization algorithm, the location of the control
points of the CAD model can be found, such that the structural performance is
improved. In order to evaluate such characteristics, a typical approach is to
minimize the compliance of the system under constraints such as displacements,
volume, etc. The statement of the structural shape optimization reads:

$\begin{array}[]{r@{}l}&\mbox{minimize}\hskip 11.0ptf_{0}(\bm{s}),\\\
&\mbox{such that}\hskip 10.0ptf_{i}(\bm{s})\leq 0,\hskip 50.0pti=1,2,...,m,\\\
&\hskip 58.0pts_{j}^{min}\leq s_{j}\leq s_{j}^{max},\hskip
9.0ptj=1,2,...,n,\end{array}$ (34)

where $f_{0}(\bm{s})$, $f_{i}(\bm{s})$ are the objective and constraints
functions respectively. In this work, we consider the NURBS control points
coordinates of the boundary description as design variables $\bm{s}$. A
standard Galerkin discretization is used in order to solve the linear
elasticity problem and evaluate the objective and constraints functions. The
weak form of the linear elasticity equation is: Find
$\bm{u}\in\mathcal{U}=\\{\bm{u}|\bm{u}\in(H^{1}(\Omega))^{3},\bm{u}=\bm{\bar{u}}\;\text{on
}\Gamma_{u}\\}$ such that for all
$\bm{\psi}\in\mathbf{\Psi}=\\{\bm{\psi}|\bm{\psi}\in(H^{1}(\Omega))^{3},\bm{\psi}=\bm{0}\;\text{on
}\Gamma_{u}\\}$:

$a_{u}(\bm{u},\bm{\psi})=l_{u}(\bm{\psi}),$ (35)

with:

$a_{u}(\bm{u},\bm{\psi})=\int_{\Omega}\bm{\epsilon}(\bm{u}):\mathbb{C}:\bm{\epsilon}(\bm{\psi})d\Omega,$
(36)

and

$l_{u}(\bm{\psi})=\int_{\Omega}\bm{b}\cdot\bm{\psi}d\Omega+\int_{\Gamma^{N}}\bm{t}\cdot\bm{\psi}d\Gamma^{N},$
(37)

where $\bm{\epsilon}$ denotes the strain, $\mathbb{C}$ is the fourth-order
elasticity tensor, $\bm{b}$ is the body force over the domain and $\bm{t}$ is
the traction over the boundary $\Gamma^{N}$. Using the isoparametric concept
the displacement vector $\bm{u}$, and the test function $\bm{\psi}$ are
discretized using the same basis functions on an element $\Omega^{e}$. Then,
the weak formulation can be written in discrete form as:

$a_{u}(\bm{u}^{h},\bm{\psi}^{h})=l_{u}(\bm{\psi}^{h}),$ (38)

where $\bm{u}^{h}$, $\bm{\psi}^{h}$ are the semidiscrete displacement and test
functions respectively. Using these relations and (38) for each tetrahedron,
we can obtain the discrete form of the balance equation:

$\bm{K}_{u}\bm{U}=\bm{t},$ (39)

where $\bm{U}$ is the vector containing the discrete values $u_{ijkl}$, and
$\bm{K}_{u}$ is the stiffness matrix for the elasticity problem, defined for
each tetrahedron as:

$\displaystyle\bm{K}^{e}_{u}$
$\displaystyle=\int_{\hat{\Omega}_{e}}\bm{B}_{u}^{T}\bm{C}\bm{B}_{u}\left|\bm{J}\right|d\hat{\Omega}_{e},$
(40)

here $\bm{C}$ is the constitutive matrix, and the matrix $\bm{B}_{u}$ is
described by:

$\bm{B}_{u}=\begin{bmatrix}\frac{\partial R^{p}_{p,0,0,0}}{\partial
x}&0&0&\cdots&\frac{\partial R^{p}_{i,j,k,l}}{\partial x}&0&0\\\
0&\frac{\partial R^{p}_{p,0,0,0}}{\partial y}&0&\cdots&0&\frac{\partial
R^{p}_{i,j,k,l}}{\partial y}&0\\\ 0&0&\frac{\partial R^{p}_{p,0,0,0}}{\partial
z}&\cdots&0&0&\frac{\partial R^{p}_{i,j,k,l}}{\partial z}\\\ \frac{\partial
R^{p}_{p,0,0,0}}{\partial y}&\frac{\partial R^{p}_{p,0,0,0}}{\partial
x}&0&\cdots&\frac{\partial R^{p}_{i,j,k,l}}{\partial y}&\frac{\partial
R^{p}_{i,j,k,l}}{\partial x}&0\\\ 0&\frac{\partial R^{p}_{p,0,0,0}}{\partial
z}&\frac{\partial R^{p}_{p,0,0,0}}{\partial y}&\cdots&0&\frac{\partial
R^{p}_{i,j,k,l}}{\partial z}&\frac{\partial R^{p}_{i,j,k,l}}{\partial y}\\\
\frac{\partial R^{p}_{p,0,0,0}}{\partial z}&0&\frac{\partial
R^{p}_{p,0,0,0}}{\partial x}&\cdots&\frac{\partial R^{p}_{i,j,k,l}}{\partial
z}&0&\frac{\partial R^{p}_{i,j,k,l}}{\partial x}\end{bmatrix},\hskip
10.0pti+j+k+l=p.$ (41)

###  5.1 Sensitivity analysis

We use a gradient based optimization algorithm (MMA), therefore the
derivatives of the objective and constraints functions with respect to the
design variables are required. In this work, we consider displacement and
compliance minimization problems. For both cases, the gradient $d\bm{u}/ds$
needs to be computed. In the following, we consider $s$ to be one entry of the
vector of design variables $\bm{s}$. The derivative of the displacement vector
with respect to the design variables can be obtained from the equilibrium
equation (39):

$\displaystyle\bm{U}_{,s}$
$\displaystyle=\bm{K}^{-1}\left(\bm{f}_{,s}-\bm{K}_{,s}\right),$ (42)

where $\bm{U}_{,s}$, $\bm{f}_{,s}$ are the gradients of the displacement and
force vectors with respect to the design variables respectively. The stiffness
matrix gradient $\bm{K},s$ is defined for each element as:

$\displaystyle\bm{K}_{e,s}$
$\displaystyle=\int_{\hat{\Omega}_{e}}\left(\bm{B}^{T}_{u,s}\bm{C}\bm{B}_{u}\left|\bm{J}\right|+\bm{B}_{u}^{T}\bm{C}\bm{B}_{u,s}\left|\bm{J}\right|+\bm{B}_{u}^{T}\bm{C}\bm{B}_{u}\left|\bm{J}\right|_{,s}\right)d\hat{\Omega}_{e},$
(43)

with the gradient $\left|\bm{J}\right|_{,s}$ obtained by:

$\left|\bm{J}\right|_{,s}=\left|\bm{J}\right|\mbox{tr}\left(\bm{J}^{-1}\bm{J}_{,s}\right).$
(44)

The sensitivity of the Jacobian is computed with:

$\bm{J}_{,s}=\begin{bmatrix}\frac{\partial x_{,s}}{\partial\xi}&\frac{\partial
y_{,s}}{\partial\xi}&\frac{\partial z_{,s}}{\partial\xi}\vspace{5pt}\\\
\frac{\partial x_{,s}}{\partial\eta}&\frac{\partial
y_{,s}}{\partial\eta}&\frac{\partial z_{,s}}{\partial\eta}\vspace{5pt}\\\
\frac{\partial x_{,s}}{\partial\zeta}&\frac{\partial
y_{,s}}{\partial\zeta}&\frac{\partial
z_{,s}}{\partial\zeta}\vspace{5pt}\end{bmatrix}.$ (45)

The gradient $\bm{B}_{u,s}$ can be computed as follows, and rearranging the
entries as in (41):

$\begin{bmatrix}\frac{\partial R_{ijkl}}{\partial x}\\\\[6.0pt] \frac{\partial
R_{ijkl}}{\partial y}\\\\[6.0pt] \frac{\partial R_{ijkl}}{\partial
z}\end{bmatrix}=\bm{J}^{-1}_{,s}\begin{bmatrix}\frac{\partial
R_{ijkl}}{\partial\xi}\\\\[6.0pt] \frac{\partial
R_{ijkl}}{\partial\eta}\\\\[6.0pt] \frac{\partial
R_{ijkl}}{\partial\zeta}\end{bmatrix},$ (46)

with

$\bm{J}^{-1}_{,s}=-\bm{J}^{-1}\bm{J}_{,s}\bm{J}^{-1},$ (47)

Recalling (25) and (26), the sensitivity $\bm{J}_{,s}$ can be written in
matrix form as:

$\bm{J}_{,s}=\bm{B}_{,\bm{\xi}}\hat{\bm{B}}^{-1}(\hat{\lambda})\hat{\bm{N}}^{e}\bm{\hat{P}_{,s}},$
(48)

with

$\bm{B}_{,\bm{\xi}}=\begin{bmatrix}\frac{\partial
R^{p}_{p,0,0,0}}{\partial\xi}&\frac{\partial
R^{p}_{0,p,0,0}}{\partial\xi}&\cdots&\frac{\partial
R^{p}_{i,j,k,l}}{\partial\xi}\vspace{5pt}\\\ \frac{\partial
R^{p}_{p,0,0,0}}{\partial\eta}&\frac{\partial
R^{p}_{0,p,0,0}}{\partial\eta}&\cdots&\frac{\partial
R^{p}_{i,j,k,l}}{\partial\eta}\vspace{5pt}\\\ \frac{\partial
R^{p}_{p,0,0,0}}{\partial\zeta}&\frac{\partial
R^{p}_{0,p,0,0}}{\partial\zeta}&\cdots&\frac{\partial
R^{p}_{i,j,k,l}}{\partial\zeta}\end{bmatrix},\hskip 30.0pti+j+k+l=p,$ (49)

and

$\bm{\hat{P}}_{,s}=\begin{bmatrix}\hat{P}^{x}_{1,s}&\hat{P}^{y}_{1,s}&\hat{P}^{z}_{1,s}\vspace{5pt}\\\
\hat{P}^{x}_{2,s}&\hat{P}^{y}_{2,s}&\hat{P}^{z}_{2,s}\vspace{5pt}\\\
\vdots&\vdots&\vdots\vspace{5pt}\\\
\hat{P}^{x}_{n,s}&\hat{P}^{y}_{n,s}&\hat{P}^{z}_{n,s}\vspace{5pt}\\\
\end{bmatrix},$ (50)

As mentioned before, we consider the design variables to be one of the
coordinates of the NURBS control points $s=\hat{P}^{i}_{j}\cdot e_{i}$. Hence,
we can compute the sensitivities $\hat{P}^{\alpha}_{\beta,s}$ by:

$\frac{\partial\hat{P}^{\alpha}_{\beta}}{\partial\hat{P}^{i}_{j}}=\begin{cases}1\hskip
5.0pt\mbox{if}\hskip 5.0pt\alpha=i\hskip 5.0pt\mbox{and}\hskip
5.0pt\beta=j,\\\ 0\hskip 5.0pt\mbox{otherwise},\end{cases}$ (51)

with $\alpha=$ 1, 2, 3 corresponding to the $x$, $y$, $z$ directions.

Figure 10: Boundary conditions for the pseudo-elastic problem applied for the
mesh moving algorithm on a 2D cantilever beam.

###  5.2 Method of moving asymptotes

The method of moving asymptotes (MMA) [34] is a non-linear optimization
algorithm within the Conservative Convex Separable Approximations (CCSA)
methods. These methods are characterized by solving a sequence of convex
subproblems at each iteration, which are approximations of the original
problem. In general, the MMA solves the problem statement from (34) by
employing the following formulation:

$\begin{array}[]{r@{}l}&\mbox{minimize}\hskip
11.0ptg_{0}^{(k)}(\bm{s})+a_{0}z+\sum_{i=1}^{m}(c_{i}y_{i}+\frac{1}{2}d_{i}y_{i}^{2}),\\\
&\mbox{such that}\hskip 10.0ptg_{i}^{(k)}(\bm{s})-a_{i}z-y_{i}\leq 0,\hskip
80.0pti=1,2,...,m,\\\ &\hskip 58.0pt\hat{\alpha}_{j}^{(k)}\leq
s_{j}\leq\hat{\beta}_{j}^{(k)},\hskip 107.0ptj=1,2,...,n,\\\ &\hskip
58.0ptz\geq 0\hskip 7.0pt\mbox{and}\hskip 7.0pty_{i}\geq 0\hskip
102.0pti=1,2,...,m,\end{array}$ (52)

where $y_{i}$, $z$ are artificial variables included to ensure the feasibility
of the problem. To approximate the original problem from (34), the values for
the constants are set $a_{0}=1$, $a_{i}=0$, $c_{i}=1000$, $d_{i}=1$ for all
$i$. The approximation functions $g_{i}$ used in the MMA are:

$g_{i}^{(k)}=\sum_{j=1}^{n}\left(\frac{p_{ij}^{(k)}}{u_{j}^{(k)}-s_{j}^{(k)}}+\frac{q_{ij}^{(k)}}{s_{j}^{(k)}-l_{j}^{(k)}}\right),\hskip
25.0pti=0,1,...,m,$ (53)

with:

$p_{ij}^{(k)}=(u_{j}^{(k)}-s_{j}^{(k)})^{2}\Bigg{(}0.001\hskip
2.0pt\mbox{max}\bigg{\\{}\frac{\partial f_{i}}{\partial
s_{j}}(\bm{s}^{(k)},0)\bigg{\\}}+\frac{10^{-5}}{s_{j}^{max}-s_{j}^{min}}\Bigg{)},$
(54)

$q_{ij}^{(k)}=(s_{j}^{(k)}-l_{j}^{(k)})^{2}\Bigg{(}0.001\hskip
2.0pt\mbox{min}\bigg{\\{}-\frac{\partial f_{i}}{\partial
s_{j}}(\bm{s}^{(k)},0)\bigg{\\}}+\frac{10^{-5}}{s_{j}^{max}-s_{j}^{min}}\Bigg{)},$
(55)

where the superscript $(k)$ denotes the iteration number. The design variable
boundaries $\alpha_{j}^{(k)}$, $\beta_{j}^{(k)}$ used for the approximation
problem are computed by:

$\begin{array}[]{r@{}l}\alpha_{j}^{(k)}&=\mbox{max}\\{s_{j}^{\mbox{min}},l_{j}^{(k)}+\theta(s_{j}^{(k)}-l_{j}^{(k)})\\},\\\
\beta_{j}^{(k)}&=\mbox{min}\\{s_{j}^{\mbox{max}},u_{j}^{(k)}+\theta(u_{j}^{(k)}-s_{j}^{(k)})\\},\end{array}$
(56)

where $l_{j}^{(k)}$, $u_{j}^{(k)}$ are the lower and upper asymptotes of the
convex approximations. The distance between the design variables and their
corresponding $\alpha_{j}^{(k)}$, $\beta_{j}^{(k)}$ bounds is restrained with
the factor $\theta$. The lower and upper asymptotes are computed for the first
two iterations as:

$\begin{array}[]{r@{}l}l_{j}^{(k)}&=s_{j}^{(k)}-\gamma^{(0)}(s_{j}^{\mbox{max}}-s_{j}^{\mbox{min}}),\\\
u_{j}^{(k)}&=s_{j}^{(k)}+\gamma^{(0)}(s_{j}^{\mbox{max}}-s_{j}^{\mbox{min}}),\end{array}$
(57)

and for $k\geq 3$ by:

$\begin{array}[]{r@{}l}l_{j}^{(k)}&=s_{j}^{(k)}-\gamma^{(k)}(s_{j}^{(k-1)}-l_{j}^{(k-1)}),\\\
u_{j}^{(k)}&=s_{j}^{(k)}+\gamma^{(k)}(u_{j}^{(k-1)}-s_{j}^{(k-1)}),\end{array}$
(58)

with the variable $\gamma^{(k)}$ as:

$\gamma^{(k)}=\begin{cases}\gamma_{\ell}\hskip 6.0pt\mbox{if}\hskip
3.0pt(s_{j}^{(k)}-s_{j}^{(k-1)})(s_{j}^{(k-1)}-s_{j}^{(k-2)})<0,\\\
\gamma_{u}\hskip 8.0pt\mbox{if}\hskip
3.0pt(s_{j}^{(k)}-s_{j}^{(k-1)})(s_{j}^{(k-1)}-s_{j}^{(k-2)})>0,\\\ 1\hskip
12.0pt\mbox{if}\hskip
3.0pt(s_{j}^{(k)}-s_{j}^{(k-1)})(s_{j}^{(k-1)}-s_{j}^{(k-2)})=0,\end{cases}$
(59)

here, $\gamma_{\ell}$, $\gamma_{u}$ are constants that control the movement of
the asymptotes depending on the direction the design variables are taking.

###  5.3 Mesh update

During the optimization process, the control points chosen as design variables
can move according to the objective, constraints and algorithm. Therefore, an
appropriate method to update the mesh is needed. A simple approach is to
regenerate the mesh at each iteration after changing the boundary of the
geometry. However, depending on the size of the problem, this might lead to an
increase of the computational time. Moreover, the mesh topology changes can
also cause a non-monotonic optimization convergence.

Algorithm 1 Shape optimization with Bézier tetratedra and a NURBS boundary
description

| Input: |  Set of NURBS surfaces describing the boundary domain
---|---|---
|  |  Initial values of the design variable vector($\bm{s}_{0}$)
|  |  Minimum and maximum limits for the design variable vector
|  |  Maximum number of iterations ($it_{max}$)
|  |  Termination tolerance $\epsilon$
|  |  Element distortion tolerance $\delta$
| Output: |  Set of NURBS surfaces describing the optimal boundary shape
|  |  Optimal design variable vector ($\bm{s}_{opt}$)
|  |  Objective function value ($f_{0}$)
|  |  Constraints values ($f_{i}(\bm{s})$)

Method:

1:Initial Boundary representation of a CAD solid

2:Represent the volume with a CAD Exchange format (IGES)

3:Generate mesh (Gmsh)

4:Read mesh file: extract nodes coordinates and topology connectivity

5:Compute Bézier control points at the boundary based on the NURBS surfaces

6:Smooth weights

7:while $\left|\frac{f_{0}^{it}-f_{0}^{{it}-1}}{f_{0}^{0}}\right|>\epsilon$
and $it\leq it_{max}$ do

8: Solve PDE: perform analysis and sensitivity analysis

9: Optimization algorithm (MMA): compute new location of design variables,
update the current boundary design

10: Solve pseudo-elastic problem: compute new location of the mesh nodes

11: Update geometry to a CAD Exchange format (IGES)

Body-fitted approaches of moving mesh problems consider the mesh as a
deformable body. The problem consists on finding the location of the node
locations of the interior of the domain given the displacement of the
boundaries. A common method used in fluid-structure interaction and
optimization problems is to solve a linear elasticity problem (also called
pseudo-elastic problem) with pure Dirichlet boundary conditions. Consider the
boundary of the domain $\Gamma=\Gamma_{d}\cup\Gamma_{s}$. Here, $\Gamma_{d}$
denotes the design boundary where the design control points are located and a
fictitious displacement $\bar{\bm{d}}$ is imposed, while on $\Gamma_{s}$, the
control points are allowed to slip along the boundary. Fig. 10 shows the
boundary conditions imposed over a 2D mesh in order to illustrate the mesh
moving problem. The pseudo-elastic problem with no body force nor traction is
as follows: Find
$\bm{d}\in\mathcal{D}=\\{\bm{d}|\bm{d}\in(H^{1}(\Omega))^{3},\bm{d}=\bm{\bar{d}}\;\text{on
}\Gamma_{d}\\}$ such that for all
$\bm{\chi}\in\mathbf{X}=\\{\bm{\chi}|\bm{\chi}\in(H^{1}(\Omega))^{3},\bm{\chi}=\bm{0}\;\text{on
}\Gamma_{d}\\}$:

$a_{d}(\bm{d},\bm{\chi})=0,$ (60)

with:

$a_{d}(\bm{d},\bm{\chi})=\int_{\Omega}\bm{\epsilon}(\bm{d}):\mathbb{C}_{mesh}:\bm{\epsilon}(\bm{\chi})d\Omega,$
(61)

where $\bm{d}$ is the displacement of the control points of the mesh, $\chi$
is the test function and with $\mathbb{C}_{mesh}$ as the fourth-order elastic
tensor for the mesh, which is defined for each element in order to control
their stiffness and retain mesh quality as much as possible. Similarly as in
the linear elasticity problem, the trial and test functions are discretized
into $\bm{d}^{h}$, $\bm{\chi}^{h}$ using the same basis functions for each
tetrahedron to obtain the discrete equation:

$a_{d}(\bm{d}^{h},\bm{\chi}^{h})=0,$ (62)

Again, introducing the discretization equations into Eq. (60), we can obtained
the discrete balance equation for the pseudo-elastic problem:

$\bm{K}_{m}\bm{D}=\bm{F}_{m}$ (63)

$\bm{K}_{m}$ and $\bm{F}_{m}$ are the fictitious stiffness matrix and force
vector respectively, and $\bm{D}$ is the vector of displacements $d_{ijkl}$.
As in [12], the same elastic matrix used in the solution of the linear
elasticity problem is used as elastic matrix for the mesh,
$\bm{C}_{mesh}=\bm{C}$. Also, the Jacobian is dropped from the definition of
$\bm{K}_{m}$, so bigger elements have lower stiffness values as:

$\displaystyle\bm{K}^{e}_{m}$
$\displaystyle=\int_{\hat{\Omega}_{e}}\bm{B}_{m}^{T}\bm{C}\bm{B}_{m}d\hat{\Omega}_{e},$
(64)

Finally, the control point coordinates of the tetrahedral mesh at the $t^{th}$
iteration is obtained by adding the displacement vector to the control points
location at the previous iteration:

$\bm{P}_{t}=\bm{P}_{t-1}+\bm{D},\hskip 10.0ptt=2,3,...,it_{max}.$ (65)

The complete procedure of the structural optimization approach is summarized
in Algorithm 1.

(a)

(b)

Figure 11: a) Dimensions and boundary conditions for the hollow sphere
validation example. b) Boundary representation of one eight of the sphere.

##  6 Numerical examples

In this section, we present some numerical examples that demonstrate the
effectiveness of the method. The first example validates the convergence of
the Bézier tetrahedron element developed. Next, numerical examples for the
solution of 3D structural shape optimization problems are presented. The
results presented were obtained using a dual Intel(R) Xeon(R) CPU E5-2650 v2 @
2.60GHz processor with 64 Gb RAM. The optimization algorithm was implemented
in Matalab 2018b, and the mesh was generated using the Gmsh version 4.2.2. For
all the examples, the termination and distortion tolerances are set
$\epsilon=10^{-4}$, $\delta=1.8$.

###  6.1 Validation test: Pressurized sphere

To validate the Bézier tetrahedron, we consider the problem of a hollow sphere
under internal pressure. The sphere has an external radius $R_{a}=4$ m and
internal radius $R_{i}=1$ m, the internal pressure is $P=1$ Pa. The material
properties used are the elasticity modulus $E=1000$ Pa, and Possin’s ratio
$\nu=0.3$. Due to symmetry, only one eighth of the domain is modelled. Fig. 11
shows the dimensions and the imposed boundary conditions. The exact stresses
using the spherical coordinate system $[r,\phi,\theta]$ are as follows:

$\begin{split}\sigma_{rr}&=\frac{PR_{i}^{3}(R_{a}^{3}-r^{3})}{r^{3}(R_{a}^{3}-R_{i}^{3})},\\\
\sigma_{\phi\phi}&=\frac{PR_{i}^{3}(R_{a}^{3}+2r^{3})}{2r^{3}(R_{a}^{3}-R_{i}^{3})},\\\
\sigma_{\theta\theta}&=\sigma_{\phi\phi}.\end{split}$ (66)

Fig. 12 shows the energy norm convergence for Lagrange tetrahedra of degree
$p=$1, 2, 3 and 4, as well as for quadratic, cubic and quartic rational Bézier
tetrahedra. The results were obtained by generating a mesh with element sizes
of $h_{l}=$ 4, 2, 1, 0.5, 0.3. With this example, we show that the accuracy of
the $C^{0}$ continuous Bézier tetrahedra is comparable to the Lagrange
elements, which is expected since both bases span the same space of piecewise
polynomials of a given degree. However, it can be seen that the $H_{1}$ norm
error is slightly lower for rational Bézier elements, especially for higher
polynomial degrees.

Figure 12: Energy norm error convergence for the hollow sphere example.

(a)

(b)

(c)

Figure 13: a) Initial set of surfaces used to describe the boundary of the
cantilever beam. b) Initial location of the design and linked control points.
c) Dimensions and boundary conditions for the problem.

(a)

(b)

(c)

(d)

Figure 14: Optimal geometries for the cantilever beam problem using cubic
Bézier tetrahedra with size mesh a) $h_{l}=2$ and c) $h_{l}=1.5$, and their
corresponding convergence curves b), d).

###  6.2 Cantilever beam

The first example dealing with shape optimization is a 3D cantilever beam. The
initial shape is constructed with 6 tensor product B-spline surfaces, one for
each face. The beam has length $L=30$ m, width $b=5$ m and an initial high
$h=6$ m. The beam is fixed on one of its faces and a vertical load is applied
on the lower edge of the free end, as shown in Fig. 13. The problem is
formulated as follows:

Table 1: Optimization results for the cantilever beam problem. Numerical experiment | Mesh update strategy
---|---
Remeshing |  | Pseudo-elastic mesh
---
iter | $f^{opt}_{0}$ |  | running
---
time (s)
iter | $f^{opt}_{0}$ |  | running
---
time (s)
$p$=2 | $h_{l}$ | 4 | 16 | 0.4922 | 269 | 23 | 0.4851 | 151
dof’s | 1203
$p$=2 | $h_{l}$ | 2 | 19 | 0.4945 | 1051 | 23 | 0.4918 | 596
dof’s | 4650
$p$=2 | $h_{l}$ | 1.5 | 19 | 0.4963 | 1595 | 26 | 0.4946 | 1175
dof’s | 7968
$p$=3 | $h_{l}$ | 4 | 20 | 0.4968 | 639 | 25 | 0.4951 | 252
dof’s | 3384
$p$=3 | $h_{l}$ | 2 | 17 | 0.4988 | 2041 | 24 | 0.4972 | 1050
dof’s | 13737
$p$=3 | $h_{l}$ | 1.5 | 17 | 0.5003 | 3143 | 24 | 0.4990 | 1931
dof’s | 23907
$p$=4 | $h_{l}$ | 4 | 18 | 0.4996 | 1092 | 25 | 0.4980 | 513
dof’s | 7275

(a)

(b)

(c)

(d)

Figure 15: a) Dimensions for the cube with a hole example for optimization. b)
Initial boundary representation. c) Initial shape for the cube with a hole
problem showing the design control points. d) Boundary conditions for an eight
of the domain.

$\begin{array}[]{r@{}l}&\mbox{minimize}\hskip
45.0ptf_{0}(\bm{s})=\frac{1}{n_{edge}}\sum_{i=1}^{n_{edge}}(u_{i}^{*})^{2},\\\
&\mbox{such that}\hskip 28.0ptV(\bm{s})\leq\bar{V},\end{array}$ (67)

where $u_{i}^{*}$ is the displacement at the $ith$ point over the edge where
the load is imposed, $n_{edge}$ refers to the number of control points on the
edge and $V(\bm{s})$ is the volume of the domain, which is constrained to a
maximum value of $\bar{V}=1050$ m3. The design variables are the $z$-component
of the control points coordinates of the top surface. The design variables are
bounded by $1\leq s_{i}\leq 10$. The MMA parameters used in this example are
$\theta=0.1$, $\gamma^{(0)}=0.5$, $\gamma_{\ell}=1.05$, $\gamma_{u}=0.65$.

Table 1 summarizes the results obtained using quadratic, cubic and quartic
Bézier tetrahedra and various mesh sizes. It can be seen that even though the
number of iterations required to achieve convergence is lower for the
remeshing approach, the use of the mesh update based on the pseudo-elastic
mesh reduces the running time considerably, in contrast with remeshing the
domain after every iteration has finished. In order to avoid degenerated
elements, especially on the free end of the beam, the mesh is updated by
calling Gmsh after 8 iterations when the pseudo-elastic approach is used. The
number of degrees of freedom shown in Table 1 refers to the initial mesh.

(a)

(b)

(c)

(d)

Figure 16: Optimal geometries of the cube with a hole example using cubic
Bézier tetrahedra with b) no refinement and c) one refinement, and their
corresponding convergence curves b), d).

###  6.3 Cube with a hole

The shape of a hole inside a cube is optimized in order to minimize the strain
energy. This example is a 3D version of the plate with a hole commonly used as
a benchmark for 2D structural shape optimization. The initial domain consists
on a cube with side length $L=4$ with an octahedral inclusion inscribed in a
circle of radius $a=1.5$ as shown in Fig. 15(a). Symmetry is considered and
therefore one eight of the cube is modelled. A total of 10 surfaces are used
for the boundary representation. This example shows that an analysis suitable
discretization can be generated even though some of the surfaces include
singularities in their parameterization. Fig. 15 show the boundary
representation, design variables and boundary conditions of the problem. A
traction $t=100$ N is imposed at the outer faces of the cube. The problem
consists on minimizing the strain energy with an imposed maximum volume
constraint $\bar{V}=62.23$ m3. The statement is formulated mathematically as:

$\begin{array}[]{r@{}l}&\mbox{minimize}\hskip
30.0ptf_{0}(\bm{s})=\bm{u}^{T}\bm{f},\\\ &\mbox{such that}\hskip
28.0ptV(\bm{s})\leq\bar{V},\end{array}$ (68)

Table 2: Optimization results for the cube with hole example, where a characteristic length $h_{l}=2$ was used over the whole domain. Numerical experiment | Mesh update strategy
---|---
Remeshing | Pseudo-elastic mesh update
iter | $f^{opt}_{0}$ |  | running
---
time (s)
iter | $f^{opt}_{0}$ |  | running
---
time (s)
p=2 | refinements | 0 | 14 | 5.3566 | 108 | 7 | 5.3395 | 31
dof’s | 1203
p=2 | refinements | 1 | 8 | 4.5746 | 293 | 7 | 4.5681 | 191
dof’s | 4650
p=3 | refinements | 0 | 11 | 5.3857 | 152 | 7 | 5.3702 | 53
dof’s | 7968
p=3 | refinements | 1 | 8 | 4.5827 | 728 | 7 | 4.5758 | 515
dof’s | 3384
Table 3: Design variables initial values and their corresponding box constraints of the hammer problem. | Design
---
variable $s_{i}$
| Initial
---
value $s^{0}_{i}$
| Minimum
---
value $s^{min}_{i}$
| Maximum
---
value $s^{max}_{i}$
$s_{1}$ | 0 | -2 | 0.8
$s_{2}$ | 0 | -2 | 0.8
$s_{3}$ | 2 | 1.2 | 3
$s_{4}$ | 2 | 1.2 | 3
$s_{5}$ | 13 | 11 | 14
$s_{6}$ | 11 | 9 | 12
$s_{7}$ | 11 | 9 | 12

The design variables are the $x$, $y$ and $z$ coordinates of the control
points describing the NURBS surface of the hole with the restrictions imposed
by symmetry. The control points at the edges of the hole are restricted to
respect symmetric conditions. In total, 27 design variables are considered,
with boundary constraints set as $0\leq s_{i}\leq 1$ for all $s_{i}$. Cubic
elements and one refinement by partitions are used to solve the governing
equations. For this example, the MMA parameters were selected as $\theta=0.1$,
$\gamma^{(0)}=0.2$, $\gamma_{\ell}=0.2$, $\gamma_{u}=0.2$.

The optimal shape approximates to a spherical inclusion as expected, in
correspondence with the 2D problem where the optimal shape resembles a
circular hole. The results of the optimization are shown in Table 2, using
quadratic and cubic Bézier tetrahedra. A constant characteristic length
$h_{l}=2$ was used to control the mesh size. Afterwards, refinement by uniform
splitting of the elements was performed. Again, the time required for
convergence with a pseudo-elastic mesh is considerably reduced in comparison
to remeshing the domain after every iteration. Moreover, it is shown that the
number of iterations is reduced and the objective function is minimized
slightly more with the pseudo-elastic approach.

The final shapes using a movable mesh with no refinements and with 1
refinement by splitting are shown in Fig. 16(a) and 16(c) respectively. A
comparison of the corresponding convergence curves is shown in Fig. 16(b) and
Fig. 16(d).

Table 4: Optimization results for the hammer example. Numerical experiment | Mesh update strategy
---|---
Remeshing | Pseudo-elastic mesh update
it | $f^{opt}_{0}$ |  | running
---
time (s)
it | $f^{opt}_{0}$ |  | running
---
time (s)
p=2 | refinements | 0 | 8 | 235.1560 | 534 | 8 | 233.9033 | 156
dof’s | 2571
p=2 | refinements | 1 | 8 | 237.9417 | 2076 | 8 | 237.4860 | 950
dof’s | 16107
p=3 | refinements | 0 | 8 | 237.8953 | 1162 | 8 | 237.4585 | 348
dof’s | 7389
p=3 | refinements | 1 | 8 | 233.7232 | 9091 | 8 | 238.5702 | 5757
dof’s | 49578

###  6.4 Hammer

Now, we consider the problem to optimize a hammer, first introduced in [15].
The problem consists in optimizing a T-shaped geometry to obtain a hammer-like
structure by minimizing the strain energy with a volume constraint
$V(\bm{s})\leq 80$ cm3. The design variables, shown on Fig. 17 and Table 3,
are the $x$ coordinates of the control points marked as $s_{1}$, $s_{2}$,
$s_{3}$ and $s_{4}$, and the $z$ coordinates for the control points $s_{5}$,
$s_{6}$ and $s_{7}$. A traction $t=100$ N is applied on one of the ends of the
T-shape, and it is fixed in the bottom face, see Fig. 17(b).

The initial geometry is constructed with 18 surface patches to describe the
boundary. Quadratic NURBS surfaces are used for the boundary representation,
and Bézier tetrahedra for the interior of the domain. The mesh is generated
using a homogeneous characteristic length $h_{l}=1$. After the mesh is
generated, one more refinement by splitting can be performed. Table 4 shows
the optimization results using quadratic and cubic elements with and without
refinement. In this case, the same MMA parameters used in the cube with a hole
example were also used here.

(a)

(b)

(c)

Figure 17: Initial shape of the hammer problem showing a) the design
variables, b) boundary conditions and c) set of surfaces used for the boundary
representation.

(a)

(b)

(c)

(d)

Figure 18: Optimal shape and convergence results for the hammer example using
cubic Bézier tetrahedra with a), b) No refinement and c), d) 1 refinement.

The results in Fig. 18 show that the optimal shape minimizes the strain energy
by approximately 29% of the initial condition. Also the amount of material
required is around 92% of the initial volume. This is in contrast with the
results shown in [15], where it was reported an increase of the compliance by
about 18% higher than the initial value, constraining the volume constraint to
2% more than the initial volume. Furthermore, the use of a movable mesh shows
again a considerable reduction in the running time to achieve convergence of
the optimization problem.

##  7 Conclusions

In this work, we presented a framework for structural shape optimization using
CAD-compatible 3D models. NURBS surfaces were used to describe the boundaries
of the geometry and Bézier tetrahedra for the discretization of the domain.
The mesh was generated using the software Gmsh, and later the elements were
converted into Bézier elements by computing the coordinates of the surface
control points with respect to the initial NURBS description. This framework
allows the modeling of complicated geometries using one patch.

An advantage of getting the control points of the Bézier tetrahedral mesh as a
function of the NURBS control points is that the sensitivity analysis for the
shape optimization problem, can be performed by implementing analytical
sensitivities in the analysis mesh with design variables located in the set of
design boundary surfaces.

Additionally, the use of a moving mesh was proposed in order to reuse the mesh
originally generated and avoid the use of discretizations with different
number of elements along the optimization procedure. This approach results in
faster convergence in terms of computational time, though in some cases with a
larger number of iterations. Moreover, for the cube with a hole example, the
use of the pseudo-elastic mesh update approach also achieves convergence with
a reduced number of iterations and minimizing slightly more the objective
function.

Possible extensions of this work include the use of $C^{1}$ continuous
elements and its implementation for other type of physical problems. The use
of stress constraints can also be considered. Moreover, improvements in the
boundary triangulation approach as well as for the mesh update can be done in
order to avoid remeshing for longer periods.

##  8 Acknowledgements

Jorge López appretiates the support of the Country-related cooperation
programme with Mexico: CONACYT-DAAD funding programme number 57177537. The
authors acknowledge Prof. Svanberg from KTH Royal Institute of Technology,
Stockholm, for providing the MMA toolbox. The research leading to these
results has received funding from the Deutsche Forschungsgemeinschaft (DFG) -
Project number 392023639.

## References

  * Alauzet [2014] F. Alauzet.  A changing-topology moving mesh technique for large displacements.  _Engineering with Computers_ , 30(2):175–200, 2014.
  * Alauzet and Loseille [2016] F. Alauzet and A. Loseille.  A decade of progress on anisotropic mesh adaptation for computational fluid dynamics.  _Computer-Aided Design_ , 72:13–39, 2016.
  * Bazilevs et al. [2006] Y. Bazilevs, L. Beirao da Veiga, J. A. Cottrell, T. J. Hughes, and G. Sangalli.  Isogeometric analysis: approximation, stability and error estimates for h-refined meshes.  _Mathematical Models and Methods in Applied Sciences_ , 16(07):1031–1090, 2006.
  * Belegundu and Rajan [1988] A. D. Belegundu and S. Rajan.  A shape optimization approach based on natural design variables and shape functions.  _Computer Methods in Applied Mechanics and Engineering_ , 66(1):87–106, 1988.
  * Belytschko et al. [2001] T. Belytschko, N. Moës, S. Usui, and C. Parimi.  Arbitrary discontinuities in finite elements.  _International Journal for Numerical Methods in Engineering_ , 50(4):993–1013, 2001.
  * Borden et al. [2011] M. J. Borden, M. A. Scott, J. A. Evans, and T. J. Hughes.  Isogeometric finite element data structures based on Bézier extraction of NURBS.  _International Journal for Numerical Methods in Engineering_ , 87(1-5):15–47, 2011.
  * Engvall and Evans [2016] L. Engvall and J. A. Evans.  Isogeometric triangular Bernstein–Bézier discretizations: Automatic mesh generation and geometrically exact finite element analysis.  _Computer Methods in Applied Mechanics and Engineering_ , 304:378–407, 2016.
  * Engvall and Evans [2017] L. Engvall and J. A. Evans.  Isogeometric unstructured tetrahedral and mixed-element Bernstein–Bézier discretizations.  _Computer Methods in Applied Mechanics and Engineering_ , 319:83–123, 2017.
  * Goldman and Filip [1987] R. N. Goldman and D. J. Filip.  Conversion from Bézier rectangles to Bézier triangles.  _Computer-Aided Design_ , 19(1):25–27, 1987.
  * Hughes et al. [2005] T. J. Hughes, J. A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, nurbs, exact geometry and mesh refinement.  _Computer methods in applied mechanics and engineering_ , 194(39-41):4135–4195, 2005.
  * Jaxon and Qian [2014] N. Jaxon and X. Qian.  Isogeometric analysis on triangulations.  _Computer-Aided Design_ , 46:45–57, 2014.
  * Johnson and Tezduyar [1994] A. A. Johnson and T. E. Tezduyar.  Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces.  _Computer methods in applied mechanics and engineering_ , 119(1-2):73–94, 1994.
  * Kadapa [2019] C. Kadapa.  Novel quadratic Bézier triangular and tetrahedral elements using existing mesh generators: Applications to linear nearly incompressible elastostatics and implicit and explicit elastodynamics.  _International Journal for Numerical Methods in Engineering_ , 117(5):543–573, 2019.
  * Lai and Schumaker [2007] M.-J. Lai and L. L. Schumaker.  _Spline functions on triangulations_.  Number 110. Cambridge University Press, 2007.
  * Lian et al. [2017] H. Lian, P. Kerfriden, and S. Bordas.  Shape optimization directly from CAD: An isogeometric boundary element approach using T-splines.  _Computer Methods in Applied Mechanics and Engineering_ , 317:1–41, 2017.
  * Liu and Jeffers [2018a] N. Liu and A. E. Jeffers.  A geometrically exact isogeometric Kirchhoff plate: Feature-preserving automatic meshing and C 1 rational triangular bézier spline discretizations.  _International Journal for Numerical Methods in Engineering_ , 115(3):395–409, 2018a.
  * Liu and Jeffers [2018b] N. Liu and A. E. Jeffers.  Rational Bézier triangles for the analysis of isogeometric higher-order gradient damage models.  In _13th World Congress on Computational Mechanics (WCCM XIII) 2nd Pan American Congress on Computational Mechanics (PANACM II)_ , pages 22–27, 2018b.
  * Liu and Korvink [2008] Z. Liu and J. G. Korvink.  Adaptive moving mesh level set method for structure topology optimization.  _Engineering Optimization_ , 40(6):529–558, 2008\.
  * López et al. [2019] J. López, C. Anitescu, N. Valizadeh, T. Rabczuk, and N. Alajlan.  Structural shape optimization using Bézier triangles and a CAD-compatible boundary representation.  _Engineering with Computers_ , pages 1–16, 2019.
  * Ludwig et al. [2019] T. Ludwig, C. Hühne, and L. De Lorenzis.  Rotation-free Bernstein–Bézier elements for thin plates and shells—development and validation.  _Computer Methods in Applied Mechanics and Engineering_ , 348:500–534, 2019.
  * Manh et al. [2011] N. D. Manh, A. Evgrafov, A. R. Gersborg, and J. Gravesen.  Isogeometric shape optimization of vibrating membranes.  _Computer Methods in Applied Mechanics and Engineering_ , 200(13-16):1343–1353, 2011.
  * Michoski et al. [2016] C. Michoski, J. Chan, L. Engvall, and J. A. Evans.  Foundations of the blended isogeometric discontinuous Galerkin (BIDG) method.  _Computer Methods in Applied Mechanics and Engineering_ , 305:658–681, 2016.
  * Möbius [1827] A. F. Möbius.  _Der barycentrische Calcul, ein Hülfsmittel zur analytischen Behandlung der Geometrie (etc.)_.  Barth, 1827.
  * Nørtoft and Gravesen [2013] P. Nørtoft and J. Gravesen.  Isogeometric shape optimization in fluid mechanics.  _Structural and Multidisciplinary Optimization_ , 48(5):909–925, 2013.
  * Park et al. [2013] B.-U. Park, Y.-D. Seo, O. Sigmund, and S.-K. Youn.  Shape optimization of the Stokes flow problem based on isogeometric analysis.  _Structural and Multidisciplinary Optimization_ , 48(5):965–977, 2013.
  * Piegl and Tiller [2012] L. Piegl and W. Tiller.  _The NURBS book_.  Springer Science & Business Media, 2012.
  * Qian [2010] X. Qian.  Full analytical sensitivities in NURBS based isogeometric shape optimization.  _Computer Methods in Applied Mechanics and Engineering_ , 199(29-32):2059–2071, 2010.
  * Schillinger et al. [2016] D. Schillinger, P. K. Ruthala, and L. H. Nguyen.  Lagrange extraction and projection for NURBS basis functions: A direct link between isogeometric and standard nodal finite element formulations.  _International Journal for Numerical Methods in Engineering_ , 108(6):515–534, 2016.
  * Sevilla et al. [2008] R. Sevilla, S. Fernández-Méndez, and A. Huerta.  NURBS-enhanced finite element method (NEFEM).  _International Journal for Numerical Methods in Engineering_ , 76(1):56–83, 2008.
  * Sevilla et al. [2011a] R. Sevilla, S. Fernández-Méndez, and A. Huerta.  NURBS-enhanced finite element method (NEFEM).  _Archives of Computational Methods in Engineering_ , 18(4):441, 2011a.
  * Sevilla et al. [2011b] R. Sevilla, S. Fernández-Méndez, and A. Huerta.  3d NURBS-enhanced finite element method (NEFEM).  _International Journal for Numerical Methods in Engineering_ , 88(2):103–125, 2011b.
  * Stein et al. [2003] K. Stein, T. Tezduyar, and R. Benney.  Mesh moving techniques for fluid-structure interactions with large displacements.  _Journal of applied mechanics_ , 70(1):58–63, 2003.
  * Stein et al. [2004] K. Stein, T. E. Tezduyar, and R. Benney.  Automatic mesh update with the solid-extension mesh moving technique.  _Computer Methods in Applied Mechanics and Engineering_ , 193(21-22):2019–2032, 2004.
  * Svanberg [1987] K. Svanberg.  The method of moving asymptotes—a new method for structural optimization.  _International journal for numerical methods in engineering_ , 24(2):359–373, 1987.
  * Tezduyar and Sathe [2007] T. E. Tezduyar and S. Sathe.  Modelling of fluid–structure interactions with the space–time finite elements: solution techniques.  _International Journal for Numerical Methods in Fluids_ , 54(6-8):855–900, 2007.
  * Thomas et al. [2015] D. C. Thomas, M. A. Scott, J. A. Evans, K. Tew, and E. J. Evans.  Bézier projection: a unified approach for local projection and quadrature-free refinement and coarsening of [NURBS and T-splines with particular application to isogeometric design and analysis.  _Computer Methods in Applied Mechanics and Engineering_ , 284:55–105, 2015.
  * Wall et al. [2008] W. A. Wall, M. A. Frenzel, and C. Cyron.  Isogeometric structural shape optimization.  _Computer methods in applied mechanics and engineering_ , 197(33-40):2976–2988, 2008.
  * Wang et al. [2018] C. Wang, S. Xia, X. Wang, and X. Qian.  Isogeometric shape optimization on triangulations.  _Computer Methods in Applied Mechanics and Engineering_ , 331:585–622, 2018.
  * Xia and Qian [2017] S. Xia and X. Qian.  Isogeometric analysis with Bézier tetrahedra.  _Computer Methods in Applied Mechanics and Engineering_ , 316:782–816, 2017.
  * Xia et al. [2015] S. Xia, X. Wang, and X. Qian.  Continuity and convergence in rational triangular Bézier spline based isogeometric analysis.  _Computer Methods in Applied Mechanics and Engineering_ , 297:292–324, 2015.
  * Yoon et al. [2013] M. Yoon, S.-H. Ha, and S. Cho.  Isogeometric shape design optimization of heat conduction problems.  _International Journal of Heat and Mass Transfer_ , 62:272–285, 2013.
  * Zareh and Qian [2019] M. Zareh and X. Qian.  Kirchhoff–Love shell formulation based on triangular isogeometric analysis.  _Computer Methods in Applied Mechanics and Engineering_ , 347:853–873, 2019.
  * Zhang and Qian [2019] R. Zhang and X. Qian.  Triangulation-based isogeometric analysis of the Cahn–Hilliard phase-field model.  _Computer Methods in Applied Mechanics and Engineering_ , 357:112569, 2019.
  * Zienkiewicz et al. [2005] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu.  _The finite element method: its basis and fundamentals_.  Elsevier, 2005.