Shape optimisation with multiresolution subdivision surfaces and immersed finite elements
Kosala Bandara, Thomas Rüberg, Fehmi Cirak∗
Department of Engineering, University of Cambridge, Trumpington Street, Cambridge CB2 1PZ, U.K.

arXiv:1510.02719v1 [math.NA] 9 Oct 2015

Abstract
We develop a new optimisation technique that combines multiresolution subdivision surfaces for boundary description with immersed finite elements for the discretisation of the primal and adjoint problems of optimisation. Similar to wavelets multiresolution
surfaces represent the domain boundary using a coarse control mesh and a sequence of detail vectors. Based on the multiresolution decomposition efficient and fast algorithms are available for reconstructing control meshes of varying fineness. During shape
optimisation the vertex coordinates of control meshes are updated using the computed shape gradient information. By virtue of
the multiresolution editing semantics, updating the coarse control mesh vertex coordinates leads to large-scale geometry changes
and, conversely, updating the fine control mesh coordinates leads to small-scale geometry changes. In our computations we start
by optimising the coarsest control mesh and refine it each time the cost function reaches a minimum. This approach effectively
prevents the appearance of non-physical boundary geometry oscillations and control mesh pathologies, like inverted elements. Independent of the fineness of the control mesh used for optimisation, on the immersed finite element grid the domain boundary is
always represented with a relatively fine control mesh of fixed resolution. With the immersed finite element method there is no
need to maintain an analysis suitable domain mesh. In some of the presented two- and three-dimensional elasticity examples the
topology derivative is used for creating new holes inside the domain.
Keywords: Shape optimisation, immersed finite elements, multiresolution surfaces, subdivision schemes, isogeometric analysis,
Catmull-Clark, b-splines

1. Introduction
We consider the shape optimisation of two- and three-dimensional solids by combining multiresolution subdivision surfaces
with immersed finite elements. As widely discussed in isogeometric analysis literature, the geometry representations used in
today’s computer aided design (CAD) and finite element analysis (FEA) software are inherently incompatible [1]. This is
particularly limiting in shape optimisation during which a given
CAD geometry model is to be iteratively updated based on the
results of a finite element computation. The inherent shortcomings of present geometry and analysis representations have
motivated the proliferation of various shape optimisation techniques. In the most prevalent approaches a surrogate geometry
model [2–8] or the analysis mesh [9, 10] instead of the true
CAD model is optimised, see also [11] and references therein.
Generally, it is tedious or impossible to map the optimised surrogate geometry model or analysis mesh back to the original
CAD model, which is essential for continuing with the design
process and later for manufacturing purposes. Moreover, geometric design features are usually defined with respect to the
CAD model and cannot be easily enforced on the surrogate
model. Recently, the shape optimisation of shells, solids and
other applications using isogeometric analysis has been explored;
∗ Corresponding

author
Email address: f.cirak@eng.cam.ac.uk (Fehmi Cirak)

that is, through directly optimising the CAD geometry model [12–
15].
In the present work the domain boundary is represented
with subdivision curves (in 2D) or surfaces (in 3D). Although
historically subdivision and related techniques have originated
in computer graphics, they recently became available in several
CAD software packages, including Autodesk Fusion 360, PTC
Creo and CATIA. As will be demonstrated in this paper, subdivision curves/surfaces provide an elegant isogeometric, bidirectional mapping between the geometry and analysis models.
In subdivision a geometry is described using a control mesh and
a limiting process of repeated refinement [16, 17]. The refinement rules are usually adapted from knot refinement rules for
b-splines [18–20]. The specific subdivision rules used in this
work are derived from cubic b-splines. For surfaces we use the
subdivision rules proposed by Catmull and Clark [20], which
lead to smooth surfaces even in case of unstructured meshes
with extraordinary vertices (i.e., domain vertices with number
of adjacent edges different than four). The hierarchy of control
meshes underlying a subdivision surface lends itself naturally
to multiresolution decomposition of geometry [21, 22]. To this
end, suitable operators are available for decomposing a geometry in a coarse control mesh and a sequence of detail vectors
similar to wavelets. Subsequently, it becomes possible to reconstruct on-the-fly control meshes of any fineness and to edit their
vertex positions. The size of the geometric region influenced
by each vertex depends on the resolution of the control mesh,
February 22, 2018

editing coarser levels leads to large-scale changes while editing finer levels lead to small-scale changes. Importantly, after
applying large-scale changes to the limit geometry the available
finer detail vectors can be automatically added to obtain the new
limit geometry, cf. Figure 1.
During gradient-based shape optimisation the shape gradient is computed to determine the boundary perturbations which
lead to a reduction in a given cost function. We consider the
adjoint approach as applied to the continuous shape optimisation problem for computing the shape gradient [23–25]. The
shape gradient is a function of the solution of the original (primal) and a complimentary adjoint boundary value problem. In
case of compliance minimisation the primal and adjoint solutions are up to the sign identical and, hence, only the primal
boundary value problem has to be solved. We discretise the primal and adjoint problems with immersed, or embedded, finite
elements, see, e.g., [26–29], which have clear advantages when
applied to structural shape optimisation [25, 30, 31]. The geometry of the domain boundary can be updated without needing to
generate a new, or to smooth, an existing domain mesh. On the
fixed non-boundary-conforming immersed finite element grid,
we use tensor-product b-splines as basis functions and enforce
Dirichlet boundary conditions with the Nitsche technique.
As mentioned, the developed multiresolution optimisation
approach relies on multiresolution subdivision curves/surfaces
for geometry description and the immersed finite element method
for domain discretisation. The multiresolution representation of
the domain boundaries allows us to describe the same geometry with control meshes of different resolution for analysis and
optimisation purposes. For finite element analysis a relatively
fine control mesh is used in order to faithfully describe the domain boundaries on the immersed grid. In contrast, the degrees
of freedom in optimisation (i.e., design variables) are chosen as
the vertex coordinates of a coarser control mesh. Usually, better shapes can be found by starting with a coarse control mesh
and optimising increasingly finer control meshes. During the
optimisation iterations the refinement level of the control mesh
is increased each time a minimum is reached. Informally, the
superior performance of the multiresolution approach can be
explained with the correlation between the number of design
variables and the number of local minima of the cost function.
Having initially fewer local minima reduces the possibility of
landing in a local minimum. It bears emphasis that with the employed, wavelet-like multiresolution decomposition, the coarse
control mesh for optimisation can be chosen independently of
the size of the features present in the geometry. For instance,
the control mesh for optimisation can be chosen much coarser
than any fillets or small-scale surface undulations present on the
to be optimised geometry.
There are a number of prior works, especially in aerodynamic shape optimisation, that use hierarchical and adaptive
geometry representations, see e.g. [32] and references therein.
For instance, in [7] Bezier basis functions and degree elevation
and in [8] b-splines and knot insertion are considered to create
a hierarchy of geometry representations. Most of these papers
primary aim to speed up the optimisation process by reducing
the number of optimisation variables or by employing multigrid

techniques. However, the added benefit of hierarchical representations in reducing the parameterisation dependency of the
final results is also often noted. In comparison to the mentioned
techniques, the advantages of multiresolution subdivision surfaces are: (i) the ability to represent geometries with arbitrary
topology, (ii) wavelet-like multiresolution representation of the
geometry, and (iii) the ease of integration with CAD packages
that use subdivision.
This paper is organised as follows. Section 2 introduces
the governing equations for linear elastic shape optimisation.
Specifically, the derivation of the continuous shape gradient
using the adjoint approach is illustrated. Subsequently, Section 3 discusses the discretisation of the primary and the adjoint boundary value problems with an immersed finite element method. Section 4 forms the core of the paper and introduces multiresolution optimisation. First the derivation of
the cubic b-spline and Catmull-Clark subdivision rules from the
well known b-spline refinement relations is demonstrated. After that, the multiresolution decomposition of subdivision surfaces and its use for multiresolution editing are shown. Upon
this basis we introduce our multiresolution shape optimisation
algorithm. Finally, Section 5 presents several two- and threedimensional examples of increasing complexity to demonstrate
the efficiency and robustness of multiresolution optimisation.
In some of the examples the topology derivative is considered
to introduce new holes in the domain.
2. Governing equations
We introduce in this section the shape optimisation of twoand three-dimensional linear elastic solids. In addition to the
linear elasticity boundary value problem a cost function that is
to be minimised is considered. For computing the shape derivatives, that is the derivatives of the cost function with respect to
the domain perturbations, we chose to use the analytic adjoint
approach. As it will become clear in Section 3, with the analytic
approach it is straightforward to compute the shape gradients
using the immersed finite element technique. In the following,
we briefly review few key results from shape calculus. See, for
instance, the monographs [23–25, 33] for a more detailed discussion.
2.1. Linear elasticity
The equilibrium equation for a solid body with the domain Ω
is given by
∇ · σ(u) + f = 0

in Ω ,

(1a)

u=0

on ΓD ,

(1b)

σ(u)n = t

on ΓN ,

(1c)

where σ is the stress tensor, u is the displacement vector, f
is the external load vector and t is the prescribed traction on
the Neumann boundary ΓN with the outward normal n. On the
Dirichlet boundary ΓD , for simplicity, only homogenous boundary conditions are prescribed.
2

(a) Original control mesh and corresponding geometry.

(b) Edited control mesh and corresponding geometry.

Figure 1: Multiresolution editing of a cylindrical component with small scale geometric details in form of small bumps. The
original geometry is shown on the left. The edited version on the right contains both large and small scale changes. The change of
the overall shape from a cylinder to a cone constitutes a large scale change and the additional bumps close to the bottom constitute
small scale changes. The bumps in the original geometry are automatically preserved during multiresolution editing.
where λ is a Lagrange parameter. L(Ω, u, λ) depends on the unknown domain Ω, the displacement field u and the Lagrange
parameter λ. The stationarity condition for the Lagrangian (6),
i.e. δL(Ω, u, λ) = 0, provides the complete set of equations
that describe the shape optimisation problem. In the following we consider one after the other the variation of L(Ω, u, λ)
with respect to the Lagrange parameter λ, displacements u and
domain Ω.

We assume a homogenous linear elastic material model
σ(u) = C : (u)

(2)

with the fourth order constitutive tensor C and linear elastic
strain tensor
1
(3)
(u) = (∇u + ∇T u) .
2
2.2. Shape optimisation
2.2.1. Formulation
The cost function to be minimised
J(Ω, u) → min

2.2.2. Primal problem
The variation of the Lagrangian L(Ω, u, λ) with respect to λ
reads
Z
∂L
∂L
δλ +
δ(∇λ) = −
δλ · [∇ · σ(u) + f ] d Ω
∂λ
∂(∇λ)
ZΩ
−
u · (C : ∇(δλ)) n d Γ
(7)
ΓD
Z
−
δλ · [t − σ(u)n] d Γ = 0 ,

(4)

depends on the to be optimised domain Ω and the unknown
solution u. The most common cost functions are integrals over
the domain Ω or its boundary Γ. For brevity and without loss
of generality, we consider in the following only the structural
compliance as the cost function
Z
J(Ω, u) =
σ(u) : (u) d Ω
Ω
Z
Z
(5)
=
f · u dΩ +
t · u dΓ.
Ω

ΓN

where we used the divergence theorem and δ(∇λ) = ∇(δλ). For
arbitrary δλ it is evident that (7) is equivalent to the linear elasticity equations (1).

ΓN

It is straightforward to adapt the subsequent derivations to other
common cost functions, such as integrals of stresses or displacements.
The minimisation of the cost function J(Ω, u) with the boundary value problem (1) as a constraint can be achieved by means
of the Lagrangian function
L(Ω, u, λ) = J(Ω, u)
Z
Z
+
∇λ : σ(u) d Ω −
λ · f dΩ
Ω
ZΩ
−
u · (C : ∇λ)n + λ · σ(u)n d Γ
ΓD
Z
−
λ · t dΓ,

2.2.3. Adjoint problem
Next, we consider the variation of the Lagrangian L(Ω, u, λ)
with respect to the displacements u
∂L
∂L
∂J
δu +
δ(∇u) =
δu
∂u
∂(∇u)
∂u
Z
+
∇λ : C : δ(∇u) d Ω
(8)
ZΩ
−
δu · (C : ∇λ)n + λ · (C : δ(∇u)) n d Γ = 0 .

(6)

ΓD

ΓN

3

After introducing the cost function (5) and reformulating the
domain term with the divergence theorem we obtain
Z
Z
Z
f · δu d Ω +
t · δu d Γ −
δu · (∇ · σ(λ)) d Ω
Ω
ΓN
Ω
Z
Z
−
λ · (C : ∇(δu)) n d Γ +
δu · σ(λ)n d Γ = 0 .
ΓD

with the determinant of the mapping
j(xt ) = det

which has, according to, e.g., [33, 34], the derivative
d
= ∇ · δv .
j(xt )
t=0
dt

ΓN

(9)

λ=0
σ(λ)n = −t

in Ω ,

(10a)

on ΓD ,

(10b)

on ΓN .

(10c)

Γ

where n is the unit normal to the boundary. Notice that this integral is zero when the perturbation direction δv is chosen tangential to the boundary. Perturbations tangential to the boundary do
not lead to a change in I1 .
Although not used in this work, we give for completeness
the variation of the boundary integral of the scalar function ψ(xt ),

By comparing this adjoint problem with the primal problem (1)
we deduce that λ = −u. Note that this holds only when the cost
function is the structural compliance (5).
2.2.4. Shape derivative
Finally, we consider the variation of L(Ω, u, λ) with respect
to the problem domain Ω, which is also referred to as the shape
derivative. To this end, we first define a linear mapping which
maps a given domain Ω into a perturbed domain Ωt , see Figure 2. With this mapping a material point with the coordinate
x ∈ Ω is mapped onto
xt = x + tδv ,

I2 (Γ) =

Z
Γ

ψ(x) d Γ .

where δv is a prescribed vector field and t is a scalar parameter.
In the usual continuum mechanics terminology, Ω is the reference configuration, Ωt is the current configuration, δv is the
material velocity vector and t is the (pseudo-) time. In shape
optimisation literature the mapping (11) is usually expressed as

where H(x) is the mean curvature on Γ, see [24, 33].
We can now write the variation of the Lagrangian L(Ω, u, λ)
in the direction of δv using the results (17) and (19). For practical shape optimisation problems in solid mechanics we usually
have only boundary variations of the form

(12)

δv̂ = 0

on ΓD ,

Before attempting the variation of L(Ω, u, λ), we first give
the variation of generic volume and surfaces integrals with respect to δv. The variation of a domain integral of the scalar
function ψ(xt )
Z
I1 (Ω) =
ψ(x) d Ω
(13)

δv̂ = 0

on ΓN with σn = t ,

δv̂ , 0

on ΓN with σn = 0 .

∂L
∂J(Ω, u)
δv̂ =
δv̂
∂Ω
Z ∂Ω

at the reference configuration Ω in the direction of δv is defined
as
Z
d
∂I1
d
δv = I1 (x + tδv)
=
ψ(xt ) d Ωt
. (14)
t=0
t=0
∂Ω
dt
dt Ωt

+

Ω

Ω

ψ(xt ) j(xt ) d Ω

ΓN

(∇λ : σ(u) − λ · f )( δv̂ · n) d Γ .

(21)

For structural compliance (5) as the cost function (i.e., λ = −u)
we obtain
Z
∂L
(2u · f − ∇u : σ(u)) (δv̂ · n) d Γ
δv̂ =
∂Ω
ΓN
Z
(22)
=
g(u)(δv̂ · n) d Γ ,

Transforming this integral into the reference configuration yields
Z

(20)

This means that only parts of the Neumann boundary ΓN with
no traction are free to move during shape optimisation. With
the result (17) at hand the variation of the Lagrangian (6) in the
direction δv̂ reads

Ω

∂I1
d
δv =
∂Ω
dt
Z
=

(18)

The variation of this integral at the reference configuration x in
the direction of δv reads
Z
∂I2
(∇ψ(x) · n + H(x)ψ(x)) (δv · n) d Γ ,
δv =
(19)
∂Γ
Γ

(11)

Ωt = Ω + tδv .

(16)

After introducing (16) into (15) and applying the divergence
theorem we obtain
Z
∂I1
(∇ψ(x) · δv + ψ(x)∇ · δv) d Ω
δv =
∂Ω
ZΩ
(17)
=
ψ(x)(δv · n) d Γ ,

The corresponding boundary value problem, referred to as the
adjoint problem, reads
∇ · σ(λ) − f = 0

∂xt
,
∂x

t=0

!
(15)
d
d
ψ(xt ) j(xt ) + ψ(xt ) j(xt ) d Ω
t=0
dt
dt

ΓN

4

tδv(x)

xt

Figure 2: Reference and the perturbed domains (left and right, respectively).
where g(u) is the shape kernel function. It is worth emphasising
that without restricting δv̂ as stated in (20), the shape derivate
would contain several more terms. Moreover, for cost functions
other than the structural compliance, the kernel function usually
is also dependent on the adjoint solution λ.
During the iterative shape optimisation the shape kernel function g(u) is used as gradient information. In order to achieve a
maximum decrease in the cost function the boundary perturbation is chosen proportional to
δv̂ = −g(u)n

3.1. Weak form of the equilibrium equations
For the linear elastic solid introduced in Section 2.1, the
weak form of the equilibrium equations (1) reads
Z
Z
Z
t · δu dΓ
σ(u) : (δu) dΩ =
f · δu dΩ +
Ω
Ω
ΓN
Z
(25)
+
t(u) · δu dΓ ,
ΓD

where δu are test functions, which are here assumed not to
be zero on the Dirichlet boundary. This assumption is necessary because we use non-boundary-fitting meshes. The weak
form (25) is not coercive and would lead, for instance, to a singular system of equations when discretised. In order to render (25) coercive we use the consistent penalty method proposed by Nitsche [37], which leads to
Z
Z
σ(u) : (δu) dΩ −
t(u) · δu dΓ
Ω
ΓD
Z
Z
γ
u · δu dΓ
−
u · t(δu) dΓ +
(26)
h ΓD
Γ
Z D
Z
t · δu dΓ
=
f · δu dΩ +

(23)

such that
∂L
δv̂ = −
∂Ω

Z
ΓN

g(u)2 d Γ .

(24)

3. Immersed finite element discretization
The shape derivatives introduced in the previous section depend on the solution of the primal and adjoint problems (1)
and (10), respectively. Although for compliance optimisation
the primal and adjoint solutions are identical (up to sign), this
is not generally the case for other cost functions. During the iterative shape optimisation, both boundary value problems have
to be repeatedly solved on constantly evolving domains. In a
conventional finite element setting this requires frequent mesh
smoothing or updating. Therefore, immersed, or embedded,
grid finite element approaches that do not require remeshing
have clear advantages in shape optimisation [25, 30, 35]. In the
present work, we use an immersed finite element technique that
we previously developed in the context of incompressible fluidstructure interaction [29, 36]. The key features of which are:
(i) the weak enforcement of Dirichlet boundary conditions with
the Nitsche technique, (ii) the use of isoparametric b-spline basis functions for discretisation, and (iii) the numerically robust
boundary and cut-cell treatment. In the following we provide a
brief summary of our discretisation method. Although we only
discuss the discretization of the primal problem (1), the same
derivations also apply to the adjoint problem (10).

Ω

ΓN

with the penalty parameter γ > 0 and the characteristic finite element size h. In contrast to conventional penalty methods, the
parameter γ in the Nitsche method is only required for numerical stability and typically a small value is sufficient.
3.2. Finite element discretisation
We use a logically Cartesian grid and the associated tensorproduct b-spline basis functions for discretizing the weak form (26).
The grid has to have the connectivity of a Cartesian grid but the
cell sizes need not to be uniform. Figure 3 shows a typical setup
in two space dimensions. The Cartesian grid facilitates the use
of tensor-product b-spline basis functions, which have a number of appealing properties known from isogeometric analysis.
Specifically, in shape optimisation the smoothness of higherorder b-splines leads to a shape kernel function (22), which is
continuous across element boundaries. This leads to optimisation algorithms that are more robust than the ones based on
5

changes the low-pass filtering of the geometry provided through
the switch to an implicit representation makes the finite element
computations more robust.
4. Multiresolution optimisation
On the logically Cartesian discretisation grid we represent
the domain boundaries either with subdivision curves or subdivision surfaces, depending on the dimensionality of the problem. The inherent hierarchy of subdivision schemes lends itself to multiresolution representation and editing of curves and
surfaces. The specific subdivision scheme that we use yields
in the curve case cubic b-splines and in the surface case cubic tensor-product b-splines [20]. Subdivision surfaces yield
smooth surfaces even for unstructured surface meshes with nontensor product structure. We refer to the monograph [16] as an
introduction to subdivision surfaces and multiresolution editing
in geometric modelling and animation.

Figure 3: Domain discretization with the immersed finite element method method. The top figures shows a (logically) Cartesian grid and a red spline curve describing the domain boundary. On the Cartesian grid the boundary is represented using the
signed distance function shown below. The zero isocontour of
the distance function provides an approximation to the original
spline curve.

4.1. Subdivision scheme for one-dimensional cubic b-splines
The refinability property of cubic b-splines can be utilised
to derive a corresponding subdivision scheme. To illustrate this,
we consider the coarse knot sequence ξi = 0, 1, 2, 3, . . . and the
corresponding fine knot sequence ξ̃i = 0, 0.5, 1, 1.5, 2, 2.5, 3, . . ..
We denote the b-splines on the coarse knot sequence with Bi (ξ)
and the ones on the fine knot sequence with B̃i (ξ). According to
the b-spline refinability equation, see, e.g., [16, 39], it is possible to represent the coarse b-splines as a linear combination of
the fine b-splines
X
Bi (ξ) =
S i j B̃ j (ξ) ,
(28)

C 0 -continuous shape functions and discontinuous shape gradients.
According to the isoparametric concept, we approximate
the domain geometry and the solution with tensor-product bsplines.
X
X
Bαi (ξ)ui , (27)
Bαi (ξ)xi , u ≈ uh (ξ) =
x ≈ xh (ξ) =
i

i

where ξ are the parametric coordinates in the knot space and i
is the control point index with the corresponding b-spline Bαi of
degree α. In optimisation computations, we usually use quadratic
b-splines with α = 2, which are C 1 -continuous and lead to continuous shape gradients. The control point coefficients xi and ui
are the nodal coordinates and displacements, respectively. The
discrete finite element equations are obtained by introducing
the approximation equations (27) into the weak form (26) and
subsequent element-by-element numerical integration. The elements which are only partially covered by the solid, so-called
cut-cells, are first triangulated prior to integration, see [36]. The
ill-conditioning of the system matrix due to small cut-cells is
avoided with an extension approach originally proposed in [38].
Finally, note that the shape gradient in the cut-cells is computed
by simply evaluating (22). The advantage of the analytic adjoint
formulation here is that the derivatives of the cut-cells with respect to the boundary position are not needed.
On the logically Cartesian grid we represent the domain
boundaries implicitly with a signed distance function (or, in
other terms, a level set). This is done despite the fact that we
have a parametric representation of the domain in form of a
multiresolution subdivision surface, see Section 4. The aim of
the switch from a parametric to an implicit representation is
to eliminate pathological geometries and topologies, like multiple crossings of a cell edge by the boundary. For optimisation problems with large boundary deformations and topology

j

where S i j is the subdivision matrix with the components

· · ·


· · ·



· · ·


S i j = · · ·



· · ·



· · ·

···

···

···

1
2
1
2

1
8
3
4
1
8

0

0

0

0

0

0

0

0

···

···

···

···

···

0

···

···

···

···

···

0

0

0

0

0

1
2
1
2

1
8
3
4
1
8

0

0

0

1
2
1
2

1
8
3
4
1
8
···

0
1
2
1
2
···


· · ·


· · ·



· · ·



· · · . (29)



· · ·



· · ·

···

Each row of this banded matrix has the same five non-zero components, always shifted by two columns relative two adjacent
rows. As shown in Figure 4 each row expresses how a coarse
b-spline can be obtained as the weighted sum of fine b-splines.
It is evident that the exact structure and components of the matrix (29) depends on the degree of the considered b-splines. The
components of the subdivision matrix for different polynomial
degrees can be found, e.g., in [29].
6

3
4
3
4
1
8
1
8

1
8
1
8
(a) Even vertex

1
2
1
2

1
2
1
(b)
2 Odd vertex

Figure 5: Subdivision stencils for cubic b-splines.
Hence, in terms of computer implementation, for a given coarse
polygon a corresponding fine polygon is obtained by first splitting each edge into two edges and subsequently computing control point coordinates according to (33) and (34). In computer
graphics literature the weights in (33) and (34) are usually given
in form of stencils shown in Figure 5.
In subdivision schemes the foregoing described approach
for obtaining a refined control polygon is applied successively
leading to finer and finer polygons. From the properties of the
b-splines underlying (33) and (34), it is clear that the control
points converge to a cubic b-spline.

Figure 4: Refinement relation for cubic b-splines.
Next, we consider a spline curve defined in terms of the
coarse b-splines and the corresponding control vertices with the
coordinates xi , i.e.,
X
xh (ξ) =
Bi (ξ)xi .
(30)
i

4.2. Catmull-Clark subdivision surfaces
Two-dimensional b-splines can be generated as the tensor
products of two one-dimensional b-splines. Similarly, the tensor product of two one-dimensional subdivision stencils yields
the two-dimensional subdivision stencils shown in Figure 6.
The three different stencils correspond to the three different type
of vertices which occur during the splitting of each face into
four faces. It can be easily verified that the weights given in
Figure 6 are the tensor-products of the one-dimensional weights
for cubic b-splines given in Figure 5. Hence, successively refining a control mesh and computing the control vertex coordinates with the stencils given in Figure 6 will lead in the limit to
a cubic spline surface.
It is evident that the tensor-product stencils only apply to
meshes in which each vertex within the domain is connected to
four faces. The number of the faces connected to a vertex is referred to as the valence of that vertex and is denoted with v. For
the sake of brevity, we refer to [16] for the discussion of regularity of vertices on the boundaries and corners. The domain vertices with a valence other than four are known as extraordinary
vertices or star-vertices. As originally proposed by Catmull and
Clark [20], the key idea in subdivision surfaces is to apply the
modified stencil shown in Figure 7 at the extraordinary vertices.
To summarise, in each subdivision refinement step each face
of the control mesh is split into four faces and the coordinates
of the control vertices are computed with the weights given in
Figures 6 and 7 depending on the local connectivity structure of
the vertex. There is mathematical theory which shows that the
resulting surface is C 2 continuous almost everywhere except at
the extraordinary vertices where it is only C 1 continuous [17].

The control polygon of the spline is obtained by linearly connecting the control points xi . Because of the refinability property (28), the spline curve (30) can also be represented with the
fine b-splines B̃ j (ξ). To show this, we introduce the refinement
relation (28) into (30)


XX
X
X


xh (ξ) =
S i j B̃ j (ξ)xi =
B̃ j (ξ)  S i j xi  . (31)
i

j

j

i

Hence, choosing the fine control vertex coordinates with
X
x̃ j =
S i j xi
(32)
i

ensures that both the coarse and fine b-splines represent the
same identical spline curve. Further inspection of (32) and subdivision matrix (29) reveals, that the fine control vertices are the
weighted averages of coarse control vertices, with the columns
of the subdivision matrix representing the weights. There are
two different sets of weights corresponding to the two different types of columns of the subdivision matrix. Based on the
numbering scheme implied in Figure 4 and the structure of the
subdivision matrix, it is easy to deduce that one set of weights
applies to control vertices with even indices
x̃2i =

3
1
1
xi−1 + xi + xi+2
8
4
8

(33)

and the other set of weights applies to control vertices with odd
indices
x̃2i+1 =

1
1
xi + xi+1 .
2
2

(34)
7

1
64

1
64

3
32

1
64

9
16

9
16

3
32

3
32

3
32

1
64

1
64

1
64

3
3
321 32
64

1
64

9
16
3
32
3
32

1
364 3
32 32

1
16

1
64

3
32

3
32

1
64

1
64

1
4

1
4

11
44

1
4

1
4

1
4

1
4

11
44

1
4

1
4

3
8

3
8
1
16

(a) Vertex

(b) Face

11
1616

1 1
16 16

33
88

11
1616

1
16

3
8
1 1
16 16

3
8
1
16

(c) Edge

Figure 6: Subdivision stencils for the Catmull-Clark scheme. Each of the stencils are used for computing the coordinates of vertices
of the type indicated by red dot.

Figure 7: Subdivision stencil for an irregular vertex with valence v with the weights β =

8

3
2v

and γ =

1
4v ,

see [20].

R maps the given control points at level ` + 1 to the control
points at the coarser level `,

In addition to the stencils shown in Figures 6 and 7 there are
also extended subdivision stencils for vertices on edges, creases
and corners, see, e.g., [40]. In this context, a crease is a line on
the surface across which the surface is only C 0 continuous. As
an illustrative example, Figure 8 shows the subdivision refinement of a control mesh for a T-junction geometry with extraordinary vertices and prescribed crease edges.
For later reference, we write the subdivision process as a
linear mapping that maps a coarse control mesh at level ` to a
finer control mesh at level ` + 1
x`+1 = Sx` ,

x` = Rx`+1 .

The coarsening matrix R is not unique and different choices are
possible. Essentially, the control mesh at level ` + 1 has more
control vertices than the one at level ` and may contain more
geometric information. In our implementation we obtain the
coarsening matrix R from a least squares fit of the subdivided
coarse control vertices Sx` to the fine control vertices x`+1 , i.e.,

(35)
x` = argminkx`+1 − Sy` k2 ,

where x`+1 and x` are two matrices containing the coordinates
of all the vertices at level ` and ` + 1, respectively. By definition
the initial coarse control mesh has the level ` = 0. The number of columns of x`+1 and x` is equal to the space dimension
and their number of rows is equal to the number of all vertices
in the mesh. The subdivision matrix S contains the weights
given by the subdivision stencils and its dimension depends on
the subdivision level ` considered. Lastly, according to (35)
the subdivision process can be interpreted as the chain of linear
mappings for obtaining increasingly finer control meshes, i.e.,
x0

S

x1

(37)

S

x2

S

···

S

x`−1

S

x`

(38)

y`

which leads to
x` = Rx`+1

with R = (ST S)−1 ST .

(39)

By comparing with (35) we can identify R as the pseudo-inverse
of the subdivision matrix S so that
 
x` = R Sx` = x` .
(40)

.

(36)
4.3. Multiresolution surface editing
The sequence of control meshes generated during subdivision refinement readily lends itself for multiresolution editing
of geometries [21, 22, 41, 42]. As discussed, Catmull-Clark
subdivision surfaces are based on cubic b-splines. Hence, the
support size of each subdivision basis function consists of a
two-ring of faces, cf. Figure 4. With increasing refinement level
the physical support size of basis functions becomes smaller.
Accordingly, depending on the refinement level the editing of
control vertex positions leads to changes with different spatial
extent on the limit surface. As an illustrative example the Tjunction geometry introduced earlier is considered in Figure 9.
In the middle column the coordinates of selected vertices at the
levels ` = 0, ` = 1 and ` = 2 are modified. As can be seen, in
the last column of pictures this leads to changes in the limit surface in the vicinity of the edited vertices and the spatial extent
of the changes is correlated with the refinement level.
The subdivision surfaces by itself do not provide the possibility to simultaneously edit coarse and fine control meshes. For
instance, after a fine control mesh is edited it is not possible anymore to edit a coarser level in order to apply large scale changes
to the geometry. Simultaneous editing of different levels can be
achieved with a wavelet-like multiresolution decomposition of
the control meshes, as will be discussed further below.
Before considering the multiresolution decomposition of control meshes, we introduce the coarsening of control meshes that
were obtained with subdivision. The linear coarsening matrix

9

In words, subdivision refinement of a control mesh followed by
coarsening (without editing) yields the original control mesh.
Instead of least squares fitting, the coarsening matrix R can
also be defined based on quasi-interpolation [43, 44] or smoothing [21]. On the other hand, coarsening by simply subsampling of the fine control mesh usually leads to artefacts in form
of oscillations in the coarse control mesh. The proposed least
squares fit approach is not very common in computer graphics
because of the need for interactivity and fast processing times.
Although the least squares matrix in (39) is sparse its solution
cannot be found at interactive rates. Notice also that (39) represents a system of equations with d right hand side vectors for
each of the d coordinate directions.
Similar to subdivision refinement the coarsening matrix can
be successively applied in order to obtain coarser representations of the geometry, i.e.,
x0

R

x1

R

x2

R

···

R

x`−1

R

x`

.

(41)
The dimension of the matrix R depends on the considered level `.
As mentioned during the coarsening process each of the coordinate directions are considered individually. Figure 10 shows
the coarsening of a subdivision surface with the described approach. From the shown limit surfaces it is visible that the
coarsening process leads to a smoothing of the geometry; and
the overall geometry is faithfully represented by the coarser representations. For this reason, the coarsening process is sometimes also referred to as the smoothing process.
With the introduced subdivision and coarsening operations
it is now possible to devise a wavelet-like multiresolution decomposition of a control mesh and the associated geometry.
The aim of this decomposition is to enable the simultaneous

(a) Coarse control mesh (x`=0 )

(b) Once subdivided control mesh (x`=1 )

(c) Limit surface (x`=∞ )

Figure 8: Subdivision refinement of a T-junction geometry. On the coarse control mesh the edges in red are tagged as crease edges.
The control mesh at the centre is obtained after one step of subdivision refinement. The geometry on the right is a rendering of the
limit surface. Notice that the limit surface is not smooth on the crease edges.

Figure 9: Multiresolution editing of the T-junction geometry introduced in Figure 8. The aim is to modify the shape and diameter
of the protruding tube. To achieve this the coordinates of selected vertices at either level ` = 0, 1 or 2 are modified as shown on
the middle column, where the vectors t 0 , t 1 and t 2 represent the modifications at the respective level. Notice the effect of the
modification level on the limit surface (last column).

10

(a) Fine control mesh (x`=3 )

(b) Twice coarsened control mesh (x`=1 )

(c) Coarse control mesh (x`=0 )

Figure 10: Subdivision coarsening of a fine control mesh. Successive application of the coarsening matrix R leads to increasingly
coarser control meshes.
editing of different levels of the subdivision surface. This is
achieved by storing a coarse control mesh and the differences
between successive control meshes called details. Subsequently
it is possible to first edit any control mesh level and then to add
the stored details as necessary. For a given fine control mesh at
level ` the detail vector d`−1 is computed using the subdivision
and coarsening processes as follows


(42)
d`−1 = x` − S Rx` = (I − SR)x` .

The corresponding synthesis step takes the form

In turn, when a control mesh at level ` − 1 and the detail d`−1
are given the geometry at level ` can be recovered according
to (42) and (39) with

For an efficient implementation of the analysis and synthesis
steps and the related data structures we refer to Zorin et al. [16].
A typical workflow during multiresolution editing of a cylindrical component with few small bumps is shown in Figure 11.
First the geometry is created starting from a coarse control mesh
and adding the bumps on the two times subdivided control mesh.
This is achieved by displacing few selected vertices on the fine
control mesh and yields the leftmost picture Figure 11. After
this step, in order to apply large scale changes it is necessary to
employ a multiresolution decomposition of the geometry. More
specifically, the control mesh x`=2 is decomposed into the detail vectors d0 and d1 , and the original control mesh x0 . For
this specific geometry the detail vector d0 is zero. After this decomposition it is possible to change the original control mesh
into, for instance, a cone and subsequently to subdivide and to
automatically add the stored details leading to the shown cone
geometry with bumps.

x` = Sx`−1 + d`−1 .

x0

d0

x1

x2

SR

R

d1

I−

d0

I−

SR

R

(43)

···

R

I−

···
SR

d`−2

R

I−

x`−1
SR

d`−1

R

I−

x1

d1

S

x2

···

S

···

d`−2

S

x`−1

S

x`

d`−1
(45)

The global detail vector d` is composed of local vertex detail
vectors, which can be conveniently stored at the vertices. In our
actual implementation, we express each local detail vector in a
local tangential coordinate system at its vertex. As known in
computer graphics, this is necessary so that in (43) any modifications to the geometry x`−1 should have an intuitive effect on
the details contained in d` , see, e.g., [21, 41, 45].
Finally, the two-level decomposition given by (42) and (43)
can be successively applied leading to a wavelet-like multiresolution decomposition of the surface. The process of obtaining
the details for a given fine geometry is referred to as the analysis
step

x0

S

4.4. Multiresolution shape optimisation
The introduced subdivision multiresolution editing technique
enables the use of two different resolutions of the same geometry for optimisation and analysis. The two resolutions correspond to different refinement levels in a multiresolution hierarchy. In shape optimisation it is usually necessary to use a
coarse control mesh for geometry updating and a relatively fine
control mesh for analysis. As is known, unwanted geometry
oscillations may appear when the analysis and geometry representations have similar resolutions [2, 10, 11]. These geometry
oscillations are usually a numerical artefact or a result from the

x`
SR

(44)

11

Figure 11: Multiresolution editing of a cylindrical component. The fine resolution mesh shown on the left at level ` = 2 is given.
In the analysis step the geometry is decomposed into the coarse level ` = 0 (second from the left) and details d0 and d1 , cf. (44).
As shown on the third from the left, the coarse level can be edited irrespective of the details. In the subsequent synthesis step (45)
adding the precomputed details to the coarse level yields the control mesh shown on the right.
Algorithm 1 the optimisation level `o is not fixed, it is increill-posedness of the considered optimisation problem. Moremented after a minimum is reached and until a user prescribed
over, in practical applications it might be desirable to optimise
only a very coarse representation out of aesthetic or manufacmaximum optimisation level is reached `o ≤ `o,max .
turability reasons.
The Algorithm 1 describes the proposed multiresolution shape
optimisation approach. The fixed computational level `c and the
maximum optimisation level `o,max are prescribed by the user.
The level `c has to be large enough such that the numerical
solution is accurate enough for practical purposes. The maximum optimisation level has to be chosen `o,max ≤ `c and determines the smallest geometric feature size on the optimised
geometry. The input control mesh x`inp can be a coarse mesh
with `inp = 0 or an already edited fine multiresolution mesh
with `inp > 0. For control meshes with `inp > 0 first a multiresolution decomposition as indicated in (44) is performed.
Throughout the algorithm the optimisation control mesh and
its level are denoted with x`o and `o , respectively. For the immersed finite element analysis the geometry corresponding to a
control mesh x`c with `c ≥ `o is used. The analysis level `c is
usually fixed and the control mesh has elements of similar size
like the cells of the immersed finite element grid. In order to
obtain the analysis control mesh x`c from the optimisation control mesh x`o we use the introduced multiresolution refinement
technique. The immersed finite element analysis yields the cost
function J(x`c , u(x`c )) and the shape kernel g(x`c ), see (22). To
compute the shape gradient at the vertices, the surface normal
vector n(x`c ) is required, which is easily computed with the
known subdivision stencils for tangent vectors, see e.g. [16].
Vertex-wise multiplication of the shape kernel with the normal vector yields the shape gradient vector g`c , which is subsequently projected to the optimisation level vector g`o by successive coarsening with R. For the sake of brevity, in Algorithm 1
the geometry is updated with a steepest descent technique and
no additional constraints are present. In applications it is common to have additional constraints, such as perimeter, area or
volume constraints. In our actual implementation we optimise
the constrained discrete problem with the method of moving
asymptotes (MMA) proposed by Svanberg [46, 47] using the
implementation in the NLopt library [48]. Finally, note that in
12

Algorithm 1 Multiresolution shape optimisation

also the topology of the domain is optimised. To this end, we
make use of the topology derivative, see e.g. [49, 50], to introduce new holes in the domain. In our current implementation
the merging or removing of holes is not considered. In all the
two-dimensional examples we use a plane strain formulation
unless otherwise indicated.

// choose maximum optimisation level `o,max and computational
level `c
// read input control mesh x`inp
1:

if `inp = 0 then

// initialise all detail vectors (d` = 0)

2:
3:
4:

else if `inp > 0 then
for ` ← `inp to 0 do
d` = x` − S(Rx` )

5.1. Shape optimisation
5.1.1. Simply supported plate with a hole
This introductory example aims to highlight the advantages
of multiresolution optimisation over classical approaches using
only one or two representation levels. The problem consists of
a square plate with an edge length L = 2 and a circular hole
with diameter D = 1, see Figure 12. The plate is loaded with
a line load of length of 1. The Young’s modulus and Poisson’s
ratio of the plate are E = 100 and ν = 0.4, respectively.
During optimisation the shape of the hole is to be modified
so that the structural compliance of the plate is minimised. The
area of the hole is constrained with
π
(46)
A(x`c ) ≥ ,
4

// Initialise optimisation level
5:

`o = 0

// Initialise cost function
6:

J=∞

// iterate over optimisation levels
7:

while `o ≤ `o,max do

// update vertex coordinates x`o while the cost function decreases
8:

repeat
// subdivide optimisation level `o up to analysis level `c

9:
10:

for ` ← `o to `c do
x` ← Sx` + d`

// compute cost function J = J(x`c , u(x`c )) and shape kernel
field g(x) = g(u(x`c ))
// compute maximum ascent direction at the vertices x`c with
the outer normals n(x`c )
11:

where x`c are the vertex coordinates of the analysis control mesh
at level `c . This constraint is necessary since the stiffest plate
is the one with a zero hole diameter. The area of the hole is
computed by integrating over its boundary
Z
1
A(x`c ) =
x`c · n(x`c ) d Γ ,
(47)
2 Γ

g`c = g(x`c )n(x`c )

// project shape derivative to optimisation level
12:
13:

for ` ← `c to `o do
g ` ← R g`

// update vertex coordinates of the optimisation level

14:
15:

where n(x`c ) is the normal to the boundary curve Γ. Recall that
we represent the boundary curve with cubic b-splines. In computations, (47) is evaluated using one quadrature point per element of the control polygon. During the optimisation also the
gradient of the area constraint is required, which is computed by
differentiating the discretised version of (47) with respect to the
vertex positions of the analysis control mesh x`c . Alternatively,
it would be possible to use analytic shape derivatives equivalent
to (17). The vertex-wise gradient of the area constraint on the
analysis level `c is projected to the optimisation level `o in the
same way as the shape gradient vector g`c . The detail vectors
created during the decomposition (44) are discarded.
Initially, at level ` = 0 the hole is represented with a cubic
spline with 8 control points. The immersed finite element grid
has 100 × 100 cells of uniform size. Three cases referred to as
C1, C2 and C3 with different geometry and analysis resolutions
are studied:

x`o ← (x`o − αg`o ) with α ≥ 0
until (Jprevious − J) < tolerance
// increment optimisation level

16:
17:

`o ← (`o + 1)
x`o ← Sx`o + d`o

5. Examples
In this section, we present several examples to demonstrate
the robustness and versatility of the proposed multiresolution
shape optimisation technique. In all the examples the domain is
described either by a cubic b-spline curve (in 2D) or a CatmullClark subdivision surface (in 3D). The objective of the optimisation is to minimise the structural compliance (5), which
is equivalent to maximising the structural stiffness. The corresponding adjoint problem (10) has (up to the sign) the same
right hand side as the primal problem (1). Hence, it is sufficient to consider only the primal problem, which is solved with
the immersed finite element technique using quadratic b-spline
basis functions. The resulting smooth stress field in combination with unique boundary normals at the vertices of the control
mesh leads to smooth shape gradients. The optimised boundary
curves or surfaces have in general no corners or sharp edges so
that there is always a unique normal.
Initially, we consider in Section 5.1 only shape optimisation
examples. Subsequently, in Section 5.2 in addition to the shape

- In C1 only one level with `o = `c = 0 is used for analysis
and optimisation.
- In C2 a four times subdivided control mesh at refinement
level `o = `c = 4 is used for analysis and optimisation.
- In C3 the optimisation level starts with `o = 0 and increases until `o = `c = 4 is reached. Throughout the
computations the analysis level is fixed to `c = 4
In case C1 the control mesh that is visible by the immersed finite element grid contains 8 elements and in cases C2 and C3
13

(a) Problem description.

(b) Optimised hole shapes for cases C1, C2 and C3.

Figure 12: Simply supported plate with a hole.
it contains 128 elements. It is clear that in case C1 the hole geometry is poorly resolved on the immersed finite element grid.
In Figure 12 the optimised final hole shapes for the three
cases are shown. In particular, the difference in optimal shapes
for cases C2 and C3 which use the same analysis level `c = 4
is striking. The case C1 is different from the other two cases
because of the mentioned inadequately coarse analysis control
mesh with `c = 0. As indicated in Figure 13, during optimisation only for case C3 the optimisation level `o is successively
increased. The optimisation level is always incremented when
a minimum is reached, cf. Algorithm 1. For the three cases the
reduction of the relative cost function over the number of iterations is shown in Figure 14. The case C2 with fixed fine resolution achieves the smallest cost reduction while the case C3 with
multiresolution achieves the largest cost reduction. The strong
dependence of the optimisation results on geometry parameterisation is well known in structural optimisation and is often associated with the non-convexity of the considered optimisation
problem. We conjecture that by initially using a coarse control
mesh for optimisation the possible number of local minima is
significantly reduced which reduces the possibility of landing
in a non-optimal local minimum. It appears that in case C2 the
optimisation problem is caught in a local minimum which is
significantly higher than the global minimum.

Figure 15: Optimal hole shapes in a two-dimensional domain.
Problem description.

5.1.2. Optimal hole shapes in a two-dimensional domain
In this prototypical example we study the optimal hole shapes
in an elastic domain subjected to bi-axial stress. The problem
setup is shown in Figure 15. As in the previous example the
area of the hole is constrained to be constant during optimisation. The Young’s modulus and Poisson’s ratio are chosen with
E = 100 and ν = 0.4, respectively. According to analytical
results for infinite plates the optimal hole shape depends on the
ratio and sign of the far-field stress [51]. When the two components of the far-field stress are of the same sign the optimal
hole shape is an ellipse with an aspect ratio r x /ry equal to the
bi-axial stress ratio σ xx /σyy . In contrast, when the two far-field
14

stress components are of opposite sign the optimal hole shape
is a quadrilateral with smooth corners. The case with far-field
stresses of the same sign has been widely studied in literature,
see, e.g., [30, 52, 53].
In our computations the initial hole geometry at level ` = 0
is modelled with a cubic spline with 8 control vertices. The
position of the vertices is chosen such that the resulting spline
curve represents approximately a circle with diameter D = 1.
The three times subdivided control mesh with 64 vertices serves
as the computational mesh for describing the hole geometry on
the immersed finite element grid. The optimisation starts with
`o = 0 and is incremented until `o = `c = 3 is reached, cf.
Algorithm 1.
In a first set of computations we quantify the effect of computing with a finite size domain as opposed to an infinite do-

Figure 13: Simply supported plate with a hole. Variation of the optimisation level `o over the number of optimisation iterations.

Figure 14: Simply supported plate with a hole. Reduction of the normalised cost over the number of optimisation iterations. The
initial cost for case C1 is 0.073 and for cases C2 and C3 is 0.065.
main underlying the analytical results. To this end, the length to
diameter ratio L/D is varied between 1.5 ≤ L/D ≤ 7 while the
hole diameter is fixed with D = 1. In all computations the element size on the immersed finite element grid is kept fixed with
1/25. We compute the optimised hole shapes for three different bi-axial stress ratios α = σ xx /σyy ∈ {0.3, 0.5, 0.7} and seven
different length to diameter ratios L/D ∈ {1.5, 2, 3, 4, 5, 7}. As
mentioned, for bi-axial stress components with the same sign
the analytically obtained optimal hole shape is an ellipse with
an aspect ratio equal to the stress ratio α. Figure 16 shows the
error in the computationally obtained ellipse aspect ratios for
different α and L/D. The error is due to the finite size of the domain and the discretisation errors. As can be seen in Figure 16
for stress ratios α ∈ {0.5, 0.7} the computationally obtained aspect ratio converges to the analytic result for sufficiently large
domains. However, for α = 0.3 the obtained aspect ratio does
not converge to the analytic result. This is due to the large discretisation errors close to the two apexes of the relatively tall
ellipse. This error could be reduced by increasing the computational level `c and decreasing the cell size of the immersed
finite element grid.
In the following set of computations the domain size and the
initial hole diameter are chosen with L = 4 and D = 1. The cell
size of the immersed finite element grid is chosen with 1/100.
According to Figure 16 this setup appears to provide sufficient
accuracy while keeping the computation time manageable. The
geometric description of the hole remains the same as in the
previous set of computations. With the described setup we compute the optimised hole shapes for ten different stress ratios

Figure 16: Optimal hole shapes in a two-dimensional domain.
The relative error in the computationally obtained aspect ratio
of the ellipse for different stress ratios α and plate sizes L/D,
with the hole diameter D = 1. The error is defined as the difference between the computationally and analytically obtained
aspect ratios.

α ∈ {−1.0, −0.7, −0.5, −0.3, −0.1, 0.1, 0.3, 0.5, 0.7, 1.0} .
15

1.2

Cherkaev et al. [51] has shown that for negative stress ratios
having several holes gives a lower compliance than having one
single hole. Therefore, it is necessary to prevent the appearance
of multiple holes which cannot be obtained with shape optimisation starting with a single hole. To regularise the optimisation
problem we add to the compliance cost function (5) an integral
penalising perimeter change, i.e.,
!2
Z
Z
J(Ω, u) =
σ(u) : (u) d Ω + ρL
dΓ .
(48)
Ω

1

0.8

0.6

0.4

Γ

0.2
0.2

where ρL is a prescribed penalty parameter. The parameter ρL
controls how much regularisation is applied to penalise the formation of multiple holes. Without this regularisation the solution of the discretised optimisation problem with the method
of moving asymptotes (MMA) does not converge. In computations we integrate the second term in (48) numerically using
the control mesh on the computational level `c . The required
gradient information is obtained by differentiating the resulting
discrete equations with respect to the vertex positions. The obtained gradient vector is added to the shape gradient vector g`c .
With cost function (48) it is now possible to compute the optimal hole shapes for positive as well as negative stress ratios.
The obtained hole shapes and aspect ratios are shown in Figure 17. For positive stress ratios the hole is of elliptical shape
and for negative ratios it is a smoothed quadrilateral. The area
of all the holes is equal due to the prescribed area constraint.
As can be seen in Figure 17, the obtained aspect ratios are in
very good agreement with the analytical result indicated by the
solid red line. The discrepancy for very small stress ratios is
due to the finite size of the domain and the appearance of tall
holes with crack-like stress concentrations requiring a finer discretisation. Finally, the effect of the penalty parameter ρL on
the obtained aspect ratio is relatively mild.

0.4

0.6

0.8

1

1.2

Figure 18: Optimal hole shapes in a three-dimensional domain.
Computationally obtained ellipsoidal hole shapes and their aspect ratios r x /rz for different stress ratios σ xx /σzz . Note that
σ xx = σyy and r x = ry (up to discretisation errors). The analytically obtained aspect ratio is shown in red.
Next, we keep the stress component σzz fixed and independently vary the two stress components σ xx and σyy . To reduce
the effect of domain size on the computationally obtained hole
shapes we chose the domain dimensions in dependence of the
stress ratios α x = σ xx /σzz and αy = σyy /σzz such that
L x = 3|α x |, Ly = 3|αy | and Lz = 3 .
In all computations the cell size is constant (5 × 5 × 5 cells per
unit volume) and is independent of the domain size.
As in the case of the preceding two-dimensional example
in Section 5.1.2, for negative stress ratios the compliance cost
function (5) is augmented with an integral penalising surface
area changes, cf. (48). The penalty modified cost function prohibits the formation of multiple holes, which cannot be obtained
with shape optimisation. The obtained hole shapes are shown
in Figure 19. The applied positive stress ratios result in ellipsoidal hole shapes with semi-axis ratios proportional to the
stress ratios (up to discretisation errors). On the other hand,
the applied negative stress ratios result in hole shapes where the
cross-section in the x − y plane is an ellipse and in the x − z and
y − z planes are smoothed quadrilaterals.

5.1.3. Optimal hole shapes in a three-dimensional domain
The considered computational domain is a cube with a side
length of L = 4 and is discretised with cells of size 1/20. The
initial geometry of the hole is a sphere with diameter D = 1
and is at level ` = 0 approximated with a control mesh of 26
nodes. The optimisation level starts with `o = 0 and is incremented until `o = `c = 3 is reached. The volume of the hole
is constrained to remain constant during optimisation. This is
achieved by computing the volume and its gradient with the
three-dimensional extension of (47).
First we choose the two stress components σ xx = σyy as
equal and only modify the σzz stress component such that

5.2. Shape and topology optimisation
This section introduces several examples with combined shape
and topology optimisation. The topology of the domain is altered by introducing new holes based on the topology derivative. Subsequently, the shape of the holes is optimised with the
presented multiresolution shape optimisation technique. This
approach is in spirit similar to the classical bubble method by
Eschenauer et al. [54]. In our present implementation we do not
consider the merging or removing of holes, hence our approach
is more restrictive than most newer topology optimisation techniques. For excellent recent reviews on various topology optimisation techniques see [35, 55, 56].
Without going in to details, the topology derivative at a
point gives the change of the cost function when a small hole

σ xx /σzz ∈ {0.3, 0.5, 0.8, 1.0}.
According to analytical results the optimised hole shapes are ellipsoids with semi-axis radius r x = ry and have the aspect ratio
r x /rz = σ xx /σzz . Figure 18 shows the computationally obtained
ellipsoidal hole shapes and their aspect ratios. The computational and analytical results agree very well especially for large
stress ratios. For smaller stress ratios the discrepancy is due to
the discretisation errors in resolving the more pronounced stress
concentrations caused by taller holes.
16

Figure 17: Optimal hole shapes in a two-dimensional domain. Computationally obtained hole shapes and their aspect ratios r x /ry
for different stress ratios α = σ xx /σyy . The analytically obtained aspect ratio is shown in red. The multiple curves for σ x /σy < 0
are computed using different penalty values ρL , cf. (48). No penalty is applied when σ x /σy > 0.
• three-dimensional elasticity

is introduced at that point, see e.g. [49, 50, 57]. To this end, in
addition to the original domain Ω a modified domain Ωr containing a hole of radius of r is considered. The hole shape is
either a circle or sphere depending on the dimension of the domain. The cost function J(Ωr , ur ) of the problem defined on Ωr
can be expressed with a series expansion
J(Ωr , ur ) = J(Ω, u) + f (r)DT J(Ω, u) + O(rd )

DT (J(Ω, u) =

(52)

We use the isocontours of the topology derivative DT (x) to
determine the location and shape of the hole to be introduced.
To this end, we introduce a control mesh which approximates
the boundary of the region to be removed from the domain. Although this step is presently performed manually, it is feasible
to automate it. This is particularly straightforward in case of
two-dimensional domains with a polygon as the boundary of
the hole. For three-dimensional problems the isocontour of the
topology derivative can be first extracted with a marching-cube
algorithm and subsequently remeshed in order obtain a uniform
control mesh with mostly regular vertices, see e.g. [60].
After a hole is generated, it is subjected to shape optimisation. During the shape optimisation the size of the hole is
allowed to increase by a prescribed amount. For instance, the
area constraint (46) is modified as follows

(49)

where f (r) is the size of the hole and DT J(Ω, u) is the topology derivative. The hole size in 2D is f (r) = r2 π and in 3D
it is f (r) = 4πr3 /3. The topology derivative can be related
to the shape derivative by considering the expansion of an infinitesimally small hole [50]. In our implementation we use the
expressions given in [58, 59] for the topology derivative of the
compliance cost function. Depending on the dimension of the
domain we obtain the topology derivative for a material with
Poisson’s ration ν with the following equations:
• two-dimensional elasticity, plane-stress
4
σ(u) : (u)
1+ν
1 − 3ν
−
Tr σ(u) Tr (u)
1 − ν2

3 1−ν 
10σ(u) : (u)
2 7 − 5ν

1 − 5ν
Tr σ(u) Tr (u)
−
1 − 2ν

DT J(Ω, u) =

(50)

A(x`c ) ≥ ρA A0 ,

(53)

where ρA is a user prescribed scalar and A0 is the area of the
initial hole. In this context the boundary of each hole is treated
as a separate multiresolution curve or surface.
In passing we note that instead of the topology derivative
it would also be possible to determine the location and shape

• two-dimensional elasticity, plane-strain
DT J(Ω, u) = 4(1 − ν)σ(u) : (u)
(51)
(1 − 4ν)(1 − ν)
−
Tr σ(u) Tr (u)
1 − 2ν
17

of the introduced holes with the density based SIMP technique
widely used in topology optimisation [61, 62].
5.2.1. Cantilever
The cantilever shown in Figure 20a is a widely studied benchmark example in topology and shape optimisation. Most relevant to our study are the results obtained by Eschenauer et
al. [54] using the bubble method, which is as previously mentioned is similar to our approach. In our computations the rectangular domain is loaded with a distributed force ty = 10 close
to the lower right corner and the Young’s modulus and Poisson’s ratio are chosen with E = 100 and ν = 0.4, respectively.
The immersed finite element grid contains 100 × 100 cells of
uniform size.

(a) Positive stress ratios σ xx /σzz > 0 and σyy /σzz > 0.

(a) Problem description. The top boundary is first optimised with
constraints on the end nodes (solid squares) to prevent horizontal
movement.

(b) Two holes are introduced at locations of minimum topology
derivative.

(b) Negative stress ratios σ xx /σzz < 0 and σyy /σzz < 0.

Figure 19: Optimal hole shapes in a three-dimensional domain.
Computationally obtained hole shapes for different stress ratios.
(c) Final optimised geometry.

Figure 20: Shape and topology optimisation of a cantilever
truss. The isocontours indicate the topology derivative.
The final topology and shape optimised geometry is shown
18

in Figure 20c. This geometry is similar to the ones presented
in [54, 63] and has been obtained in three steps.
1. The first step is a shape optimisation step. During optimisation only the top boundary of the plate is allowed
to move based on the computed shape gradients. In the
optimisation level `o = 0 the top boundary is represented
with a control polygon using 2 elements and 3 vertices. In
the twice subdivided computational level `c = 2 the control polygon contains 16 elements. During subdivision
the two end nodes are tagged as corner so that they do
not move horizontally. In addition, during iterative shape
optimisation the right corner node is allowed to move
only vertically. Moreover, we apply area constraints of
the form (53) so that the domain size is reduced. The
geometry obtained with shape optimisation is shown in
Figure 20b.
2. The second step is a topology optimisation step. After computing the topology derivative we manually introduce two holes at locations with minimum topology
derivative, see Figure 20b. The first hole is triangular
shaped and splits the clamped boundary into an upper
and lower part. The second hole is square shaped and
is located inside the domain. The holes have to be large
enough so that they can be represented on the immersed
finite element grid. According to [36] the minimum hole
size has to be larger than
√
(54)
2 d(α + 1)h ,

Figure 21: Simply supported plate with aspect ratio H/L. Each
of the two supports have a width of 0.05L. The two bold vertical
lines indicate the free boundaries which are allowed to move
during shape optimisation. The two ends of each of the vertical
lines (red squares) are constraint to move only horizontally.
optimisation step is followed by several topology optimisation
and shape optimisation steps. In the initial shape optimisation
step the two vertical boundaries of the plate are optimised while
the domain area is allowed to reduce. At the coarse optimisation level `o = 0 each edge is represented with two elements,
which are three times subdivided to obtain the computational
control polygon at level `c = 3. In the subsequent topology optimisation step we semi-manually introduce triangle and square
shaped holes, see Figure 22 middle column. For each of the
polygons the optimisation and computation levels are chosen
with `o = 0 and `c = 3, respectively. The topology optimisation
is followed by the shape optimisation of all the domain boundaries. The movement of the vertices under the distributed force
and supports are constrained to remain fixed. In case of the
plate with small H/L = 0.25 several more topology and shape
optimisation steps are performed, see Figure 22 middle column.
The obtained geometries are shown in Figure 22 right column.
The optimisation is always terminated before any holes start to
merge. As expected, we obtain for a H/L = 1.0 an arch-like
structure and for H/L = 0.25 a truss-like structure.

where d = 2 is the dimensionality of the domain, α is the
polynomial degree of the b-spline basis functions and h
is the cell size.
3. The last step is again a shape optimisation step. The control polygons belonging to the previously introduced two
holes with three and four nodes, respectively, are subdivided twice to obtain the computational control polygon
at level `c = 2. Subsequently the geometry of the two
holes is iteratively optimised using shape optimisation.
The optimisation is terminated before the two holes start
to merge.
5.2.2. Simply supported plates with different aspect ratios
We consider three simply supported plates with different aspect ratios H/L, see Figure 21. It is known from structural analysis that for slender plates a truss-like system and for stocky
plates an arch-like system is more efficient. The plate is loaded
with a symmetrically placed distributed vertical force ty = 1
of length L/5 on its top edge. The Young’s modulus and Poisson’s ratio are chosen with E = 100 and ν = 0.4, respectively.
We consider three different aspect ratios H/L ∈ {0.25, 0.5, 1}
with the corresponding immersed finite element grids containing 80 × 200, 200 × 100 and 150 × 150 cells, respectively. In the
computations the height is fixed to H = 1 and only the length L
is varied.
As in the cantilever example, cf. Section 5.2.1, the final
geometry is obtained in several steps, namely an initial shape

5.2.3. Three-dimensional stool
In this last example we present the combined topology and
shape optimisation of a three-dimensional solid, see Figure 23.
The initial domain is a truncated pyramid and is at its top loaded
with a uniform distributed load tz = 10. At its bottom it is
supported by four distributed roller supports each of size 0.2 ×
0.2. The Young’s modulus and Poisson’s ratio are chosen with
E = 100 and ν = 0.4, respectively.
In the optimisation study only one quarter of the domain is
considered and appropriate bounds and geometry tags are applied at the two planes of symmetry. The corresponding immersed finite element grid is of size 0.7 × 0.7 × 1 and consists
19

Figure 22: Simply supported plates with different aspect ratios. In all snapshots the isocontours indicate the topology derivative. In
the middle column the small triangular and square shaped holes introduced during the topology optimisation step are shown.

20

geometry is always described with a fine control mesh on the
immersed finite element grid, independently from the control
mesh level used for optimisation. Hence, any fine scale geometric details, like fillets or surface undulations, are faithfully
represented on the discretisation grid. As our examples demonstrate multiresolution shape optimisation enables us to find better optima and is exceedingly robust, partly due to the use of
the immersed finite element technique.
The multiresolution editing semantics appears to be particularly appealing for isogeometric analysis because it enables
the full decoupling of the geometry and the analysis representations of the same geometry. It allows to seamlessly map variables and fields between the two representations irrespective
of their resolutions. Beyond optimisation this decoupling can
be useful in a number of applications, such as for computing
fast approximate solutions or multigrid and multilevel preconditioners [32]. The presented multiresolution techniques can
also be extended to non-uniform rational b-splines (NURBS),
which are the more commonly used basis functions in industrial software. It is straightforward to include rational b-splines
through the use of homogenised coordinates, see, e.g., [64]. In
order to consider non-uniform b-splines it is instructive to consult previous works on non-uniform subdivision [65, 66]. Although we focused in the present paper on uniform subdivision
refinement and coarsening, it is conceivable (and desirable) to
develop adaptive multiresolution algorithms in the spirit of hierarchical b-splines [67–69]. The utility of adaptive geometry
representations in shape optimisation has already been demonstrated, for instance, with b-spline surfaces refined by knot insertion [8] and reparameterisation of Bezier surfaces [70]. Furthermore, in the present paper in 3D holes were introduced
by manually fitting control meshes to the isocontours of the
topology derivative. This process can be automated using techniques for extracting subdivision control meshes from isocontours [71]. The related issue of merging of holes can be achieved
with approximate Boolean operations for subdivision surfaces [72].
In closing, it is worth emphasising that basic subdivision techniques recently became available in a number of engineering
design software, including Autodesk Fusion 360, PTC Creo and
CATIA, which most likely will increase their use in future engineering practice.

Figure 23: Three-dimensional stool. Problem description.
Roller supports are applied to all finite element nodes inside
the regions of size 0.2 × 0.2 marked by dashed squares.
of 30 × 30 × 30 cells.
The sequence of the performed topology and shape optimisation steps are shown in Figures 24, 25, 26 and 27. In total
two topology and two shape optimisation steps are performed.
In each topology optimisation step we remove in one go a relatively large amount of material by deleting computational cells
with topology derivative below a threshold. In the first topology
optimisation step illustrated in Figure 24 all cells with topology
derivative DT J(Ω, u) ≤ 0.025 are removed. Subsequently, we
semi-manually generate the coarse resolution subdivision control mesh depicted in Figure 24b for representing the new topology. In the following shape optimisation step, see Figure 25, the
generated control mesh serves as the optimisation level `o = 0
and the computation level is chosen with `c = 2. During the
shape optimisation the volume of the domain is constraint to remain constant. In the second topology optimisation step shown
in Figure 26 all cells with topology derivative DT J(Ω, u) ≤ 0.04
are removed. This is followed by a semi-manual control mesh
generation, see Figure 26b, and the final shape optimisation step
shown Figure 27.
6. Summary and conclusions
A multiresolution optimisation technique based on subdivision surfaces was introduced. The domain boundaries are
described with subdivision surfaces and the domain boundary
value problem is discretised with an immersed finite element
technique. The wavelet-like multiresolution decomposition of
the domain boundary yields a low resolution control mesh and
a sequence of detail vectors. The control mesh at any specific refinement level can be reconstructed on-the-fly with the
introduced synthesis operator. Editing the coarse levels leads
to large-scale geometry changes while editing fine levels leads
to small-scale geometry changes. In addition, the multiresolution editing semantics allows the decoupling of the choice of
the editing level from the size of the geometric features present
in the geometry. In the proposed approach we start optimising
the coarsest control mesh and successively increase the optimisation level each time a minimum is reached. The domain

Acknowledgement
The partial support of the EPSRC through grant # EP/G008531/1
and EC through Marie Curie Actions (IAPP) program CASOPT
project are gratefully acknowledged.
References
[1] 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 (2005) 4135–4195.
[2] V. Braibant, C. Fleury, Shape optimal design using b-splines, Computer
Methods in Applied Mechanics and Engineering 44 (1984) 247–267.
[3] R. T. Haftka, R. V. Grandhi, Structural shape optimization — a survey,
Computer Methods in Applied Mechanics and Engineering 57 (1986) 91–
106.

21

(a) Initial geometry.

(b) Optimised geometry.

Figure 24: Three-dimensional stool. First topology optimisation step. The isocontours indicate the topology derivative. The
wireframe represents the coarse control mesh.

(a) Initial geometry.

(b) Optimised geometry.

Figure 25: Three-dimensional stool. First shape optimisation step. The wireframe represents the coarse resolution mesh used for
optimisation.

(a) Initial geometry.

(b) Optimised geometry.

Figure 26: Three-dimensional stool. Second topology optimisation step. The isocontours indicate the topology derivative. The
wireframe represents the coarse control mesh.

22

(a) Initial geometry.

(b) Optimised geometry.

Figure 27: Three-dimensional stool. Second shape optimisation step. The wireframe represents the coarse resolution mesh used for
optimisation.

[4] K.-U. Bletzinger, S. Kimmich, E. Ramm, Efficient modeling in shape optimal design, Computing Systems in Engineering 2 (1991) 483–495.
[5] N. Olhoff, M. P. Bendsøe, J. Rasmussen, On CAD-integrated structural
topology and design optimization, Computer Methods in Applied Mechanics and Engineering 89 (1) (1991) 259–279.
[6] F. Beux, A. Dervieux, A hierarchical approach for shape optimization,
Engineering Computations 11 (1) (1994) 25–48.
[7] B. A. El Majd, J.-A. Désidéri, R. Duvigneau, Multilevel strategies for
parametric shape optimization in aerodynamics, European Journal of
Computational Mechanics/Revue Européenne de Mécanique Numérique
17 (1-2) (2008) 149–168.
[8] X. Han, D. W. Zingg, An adaptive geometry parametrization for aerodynamic shape optimization, Optimization and Engineering 15 (1) (2014)
69–91.
[9] K.-U. Bletzinger, M. Firl, J. Linhard, W. R., Optimal shapes of mechanically motivated surfaces, Computer Methods in Applied Mechanics and
Engineering 199 (2010) 324–333.
[10] C. Le, T. Bruns, D. Tortorelli, A gradient-based, parameter-free approach
to shape optimization, Computer Methods in Applied Mechanics and Engineering 200 (9) (2011) 985–996.
[11] K.-U. Bletzinger, A consistent frame for sensitivity filtering and the vertex assigned morphing of optimal shape, Structural and Multidisciplinary
Optimization 49 (6) (2014) 873–895.
[12] F. Cirak, M. Scott, E. Antonsson, M. Ortiz, P. Schröder, Integrated modeling, finite-element analysis, and engineering design for thin-shell structures using subdivision, Computer Aided Design 34 (2002) 137–148.
[13] 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.
[14] W. Wall, M. Frenzel, C. Cyron, Isogeometric structural shape optimizisation, Computer Methods in Applied Mechanics and Engineering 197
(2008) 2976–2988.
[15] K. Bandara, F. Cirak, G. Of, O. Steinbach, J. Zapletal, Boundary element based multiresolution shape optimisation in electrostatics, Journal
of Computational Physics 297 (2015) 584–598.
[16] D. Zorin, P. Schröder, Subdivision for Modeling and Animation, SIGGRAPH 2000 Course Notes, 2000.
[17] J. Peters, U. Reif, Subdivision Surfaces, Springer Series in Geometry and
Computing, Springer, 2008.
[18] J. M. Lane, R. F. Riesenfeld, A theoretical development for the computer
generation and display of piecewise polynomial surfaces, IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI-2 (1980) 35–
46.
[19] D. Doo, M. Sabin, Behavior of recursive division surfaces near extraordinary points, Computer-Aided Design 10 (1978) 356–360.
[20] E. Catmull, J. Clark, Recursively generated B-spline surfaces on arbitrary
topological meshes, Computer-Aided Design 10 (1978) 350–355.
[21] D. Zorin, P. Schröder, W. Sweldens, Interactive multiresolution mesh edit-

[22]

[23]
[24]
[25]
[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]
[35]

[36]

[37]

23

ing, in: SIGGRAPH 1997 Conference Proceedings, ACM Press/AddisonWesley Publishing Co., Los Angeles, CA, 259–268, 1997.
M. Lounsbery, T. D. DeRose, J. Warren, Multiresolution analysis for surfaces of arbitrary topological type, ACM Transactions of Graphics 16
(1997) 34–73.
E. J. Haug, K. K. Choi, V. Komkov, Design Sensitivity Analysis of Structural Systems, Academic Press, 1986.
J. Sokolowski, J. P. Zolesio, Introduction to Shape Optimization: Shape
Sensitivity Analysis, Springer, 1992.
J. Haslinger, R. A. E. Mäkinen, Introduction to Shape Optimization: Theory, Approximation, and Computation, SIAM, Philadelphia, 2003.
T. Rüberg, F. Cirak, An immersed finite element method with integral
equation correction, International Journal for Numerical Methods in Engineering 86 (2011) 93–114.
R. Sanches, P. Bornemann, F. Cirak, Immersed b-spline (i-spline) finite
element method for geometrically complex domains, Computer Methods
in Applied Mechanics and Engineering 200 (2011) 1432–1445.
D. Schillinger, L. Dede, M. A. Scott, J. A. Evans, B. M. J., E. Rank,
T. J. R. Hughes, An isogeometric design-through-analysis methodology
based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Computer Methods in Applied
Mechanics and Engineering 249–253 (2012) 116–150.
T. Rüberg, F. Cirak, Subdivision-stabilised immersed b-spline finite elements for moving boundary flows, Computer Methods in Applied Mechanics and Engineering 209–212 (2012) 266–283.
J. Norato, R. Haber, D. Tortorelli, M. Bendsoe, A geometry projection method for shape optimization, International Journal For Numerical
Methods In Engineering 60 (2004) 2289–2312.
G. Allaire, F. Jouve, A.-M. Toader, Structural optimization using sensitivity analysis and a level-set method, Journal of Computational Physics
194 (2004) 363–393.
P. F. Antonietti, A. Borzı̀, M. Verani, Multigrid shape optimization governed by elliptic PDEs, SIAM J. Control Optim. 51 (2) (2013) 1417–
1440, ISSN 0363-0129.
M. Delfour, J.-P. Zolésio, Shapes and Geometries: Analysis, Differential Calculus, and Optimization, Advances in Design and Control, SIAM,
Philadelphia, 2001.
G. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for
Engineering, Wiley-Blackwell, 2000.
N. P. van Dijk, K. Maute, M. Langelaar, F. Van Keulen, Level-set methods
for structural topology optimization: a review, Structural and Multidisciplinary Optimization 48 (3) (2013) 437–472.
T. Rüberg, F. Cirak, A fixed-grid b-spline finite element technique for
fluid–structure interaction, International Journal for Numerical Methods
in Fluids 74 (9) (2014) 623–660.
J. Nitsche, Über ein Variationsprinzip zur Lösung von DirichletProblemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind, Abhandlungen aus dem Mathematischen Seminar
der Universität Hamburg 36 (1971) 9–15.

[38] K. Höllig, Finite Element Methods with B-splines, SIAM Frontiers in
Applied Mathematics, 2003.
[39] H. Prautzsch, W. Boehm, M. Paluszny, Bezier and B-Spline Techniques,
Springer Verlag, 2002.
[40] L. Ying, D. Zorin, Nonmanifold Subdivision, in: IEEE Visualization
2001, San Diego, CA, 325–332, 2001.
[41] A. Finkelstein, D. H. Salesin, Multiresolution curves, in: Proceedings of
the 21st annual conference on Computer graphics and interactive techniques, ACM, 261–268, 1994.
[42] S. Gortler, M. Cohen, Hierarchical and variational geometric modeling
with wavelets, in: Proceedings of the 1995 symposium on interactive 3D
graphics, 35–42, 1995.
[43] C. de Boor, The quasi-interpolant as a tool in elementary polynomial
spline theory, Approximation Theory (1973) 269–276.
[44] N. Litke, A. Levin, P. Schröder, Fitting subdivision surfaces, in: IEEE
Visualization 2001, San Diego, CA, 319–324, 2001.
[45] O. Sorkine, D. Cohen-Or, Y. Lipman, M. Alexa, C. Rössl, H.-P. Seidel, Laplacian surface editing, in: Proceedings of the 2004 Eurographics/ACM SIGGRAPH symposium on Geometry processing, ACM, 175–
184, 2004.
[46] K. Svanberg, The method of moving asymptotes — a new method for
structural optimization, International Journal for Numerical Methods in
Engineering 24 (1987) 359–373.
[47] K. Svanberg, A class of globally convergent optimization methods based
on conservative convex separable approximations, SIAM Journal on Optimization 12 (2002) 555–573.
[48] S. G. Johnson, The NLopt nonlinear-optimization package, http://abinitio.mit.edu/nlopt.
[49] H. A. Eschenauer, N. Olhoff, Topology optimization of continuum structures: A review, Applied Mechanics Reviews 54 (4) (2001) 331–390.
[50] A. A. Novotny, R. Feijóo, E. Taroco, C. Padra, Topological sensitivity analysis, Computer Methods in Applied Mechanics and Engineering
192 (7) (2003) 803–829.
[51] A. Cherkaev, Y. Grabovsky, A. Movchan, S. Serkov, The cavity of the
optimal shape under the shear stresses, International Journal of Solids
and structures 35 (33) (1998) 4391–4410.
[52] E. S. Kristensen, N. F. Madsen, On the optimum shape of fillets in plates
subjected to multiple in-plane loading cases, International Journal for Numerical Methods in Engineering 10 (5) (1976) 1007–1019.
[53] E. Cervera, J. Trevelyan, Evolutionary structural optimisation based on
boundary representation of NURBS. Part I: 2D algorithms, Computers &
Structures 83 (23) (2005) 1902–1916.
[54] H. A. Eschenauer, V. V. Kobelev, A. Schumacher, Bubble method for
topology and shape optimization of structures, Structural optimization
8 (1) (1994) 42–51.
[55] O. Sigmund, K. Maute, Topology optimization approaches, Structural and
Multidisciplinary Optimization 48 (6) (2013) 1031–1055.
[56] J. D. Deaton, R. V. Grandhi, A survey of structural and multidisciplinary
continuum topology optimization: post 2000, Structural and Multidisciplinary Optimization 49 (1) (2014) 1–38.
[57] S. Garreau, P. Guillaume, M. Masmoudi, The topological asymptotic for
PDE systems: the elasticity case, SIAM journal on control and optimization 39 (6) (2001) 1756–1778.
[58] A. Novotny, R. Feijóo, E. Taroco, C. Padra, Topological-shape sensitivity
method: theory and applications, in: IUTAM Symposium on Topological Design Optimization of Structures, Machines and Materials, Springer,
469–478, 2006.
[59] A. Novotny, R. Feijóo, E. Taroco, C. Padra, Topological sensitivity analysis for three-dimensional linear elasticity problem, Computer Methods in
Applied Mechanics and Engineering 196 (41) (2007) 4354–4364.
[60] M. Botsch, L. Kobbelt, M. Pauly, P. Alliez, B. Lévy, Polygon mesh processing, CRC press, 2010.
[61] M. P. Bendsoe, Optimal shape design as a material distribution problem,
Structural optimization 1 (4) (1989) 193–202.
[62] M. P. Bendsoe, O. Sigmund, Topology optimization: theory, methods and
applications, Springer, 2003.
[63] M. Y. Wang, X. Wang, D. Guo, A level set method for structural topology
optimization, Computer methods in applied mechanics and engineering
192 (1) (2003) 227–246.
[64] F. Cirak, Q. Long, Subdivision shells with exact boundary control and
non-manifold geometry, International Journal for Numerical Methods in

Engineering 88 (2011) 897–923.
[65] T. Sederberg, J. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCs,
in: SIGGRAPH 2003 Conference Proceedings, San Diego, CA, 477–484,
2003.
[66] T. Cashman, U. Augsdörfer, N. Dodgson, M. Sabin, NURBS with extraordinary points: high-degree, non-uniform, rational subdivision schemes,
in: SIGGRAPH 2009 Conference Proceedings, New Orleans, LA, 46:1–
46:9, 2009.
[67] E. Grinspun, P. Krysl, P. Schröder, CHARMS: A simple framework for
adaptive simulation, in: SIGGRAPH 2002 Conference Proceedings, San
Antonio, TX, 281–290, 2002.
[68] P. B. Bornemann, F. Cirak, A subdivision-based implementation of the hierarchical b-spline finite element method, Computer Methods in Applied
Mechanics and Engineering 253 (2013) 584–598.
[69] X. Wei, Y. Zhang, T. J. Hughes, M. A. Scott, Truncated hierarchical
Catmull–Clark subdivision with local refinement, Computer Methods in
Applied Mechanics and Engineering 291 (2015) 1–20.
[70] J.-A. Désidéri, B. A. El Majd, A. Janka, Nested and self-adaptive
Bézier parameterizations for shape optimization, Journal of Computational Physics 224 (1) (2007) 117–131.
[71] D. Bommes, B. Lévy, N. Pietroni, E. Puppo, C. Silva, M. Tarini, D. Zorin,
Quad-Mesh Generation and Processing: A Survey, Computer Graphics
Forum 32 (6) (2013) 51–76.
[72] H. Biermann, D. Kristjansson, D. Zorin, Approximate Boolean operations
on free-form solids, in: Proceedings of the 28th annual conference on
computer graphics and interactive techniques, 185–194, 2001.

24

