Shape optimisation with multiresolution subdivision surfaces and immersed finite elements

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.


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][3][4][5][6][7][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; that is, through directly optimising the CAD geometry model [12][13][14][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][19][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, 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][24][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][27][28][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.

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][24][25]33] for a more detailed discussion.

Linear elasticity
The equilibrium equation for a solid body with the domain Ω is given by 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.  We assume a homogenous linear elastic material model with the fourth order constitutive tensor C and linear elastic strain tensor

Shape optimisation 2.2.1. Formulation
The cost function to be minimised 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 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) 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 Ω.

Primal problem
The variation of the Lagrangian L(Ω, u, λ) with respect to λ reads where we used the divergence theorem and δ(∇λ) = ∇(δλ). For arbitrary δλ it is evident that (7) is equivalent to the linear elasticity equations (1).
After introducing the cost function (5) and reformulating the domain term with the divergence theorem we obtain The corresponding boundary value problem, referred to as the adjoint problem, reads 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).

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 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 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 ψ(x t ) at the reference configuration Ω in the direction of δv is defined as Transforming this integral into the reference configuration yields with the determinant of the mapping which has, according to, e.g., [33,34], the derivative After introducing (16) into (15) and applying the divergence theorem we obtain 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 I 1 .
Although not used in this work, we give for completeness the variation of the boundary integral of the scalar function ψ(x t ), The variation of this integral at the reference configuration x in the direction of δv reads 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 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 For structural compliance (5) as the cost function (i.e., λ = −u) we obtain x t tδv(x) 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 (23) such that

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).

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 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 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.

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 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 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 x i and u i 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 changes the low-pass filtering of the geometry provided through the switch to an implicit representation makes the finite element computations more robust.

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.

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 B i (ξ) and the ones on the fine knot sequence withB 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 where S i j is the subdivision matrix with the components · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · 1 2 · · · · · · · · · · · · · · · · · · · · · · · · · · · 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]. Next, we consider a spline curve defined in terms of the coarse b-splines and the corresponding control vertices with the coordinates x i , i.e., The control polygon of the spline is obtained by linearly connecting the control points x i . Because of the refinability property (28), the spline curve (30) can also be represented with the fine b-splinesB j (ξ). To show this, we introduce the refinement relation (28) into (30) Hence, choosing the fine control vertex coordinates with 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 and the other set of weights applies to control vertices with odd indicesx 2i+1 =  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.

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].  , see [20].
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 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., x 0

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 R maps the given control points at level + 1 to the control points at the coarser level , 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., which leads to By comparing with (35) we can identify R as the pseudo-inverse of the subdivision matrix S so that 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., 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 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. 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 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 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 The corresponding synthesis step takes the form 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 d 0 and d 1 , and the original control mesh x 0 . For this specific geometry the detail vector d 0 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.

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 In the analysis step the geometry is decomposed into the coarse level = 0 (second from the left) and details d 0 and d 1 , 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.
ill-posedness of the considered optimisation problem. Moreover, in practical applications it might be desirable to optimise only a very coarse representation out of aesthetic or manufacturability 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 for ← o to c do 10: 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: g c = g(x c )n(x c ) // project shape derivative to optimisation level 12: for ← c to o do 13: g ← R g // update vertex coordinates of the optimisation level 14: x o ← (x o − αg o ) with α ≥ 0 15: until (J previous − J) < tolerance // increment optimisation level 16:

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 Catmull-Clark 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 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.

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 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 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:  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.

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 /r y equal to the bi-axial stress ratio σ xx /σ yy . In contrast, when the two far-field 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-  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  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.
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., 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.

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 σ 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 = r y and have the aspect ratio r x /r z = σ 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. Computationally obtained ellipsoidal hole shapes and their aspect ratios r x /r z for different stress ratios σ xx /σ zz . Note that σ xx = σ yy and r x = r y (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 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.

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 Figure 17: Optimal hole shapes in a two-dimensional domain. Computationally obtained hole shapes and their aspect ratios r x /r y 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.
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 , u r ) of the problem defined on Ω r can be expressed with a series expansion where f (r) is the size of the hole and D T J(Ω, u) is the topology derivative. The hole size in 2D is f (r) = r 2 π and in 3D it is f (r) = 4πr 3 /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 • two-dimensional elasticity, plane-strain • three-dimensional elasticity We use the isocontours of the topology derivative D T (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 where ρ A is a user prescribed scalar and A 0 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 (a) Positive stress ratios σ xx /σ zz > 0 and σ yy /σ zz > 0.
of the introduced holes with the density based SIMP technique widely used in topology optimisation [61,62].

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 t y = 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.  The final topology and shape optimised geometry is shown 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 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.

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 t y = 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 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.

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 t z = 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 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. 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 D T 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 D T 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.

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 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][68][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.