{'fraction_non_alphanumeric': 0.02646304499883423, 'fraction_numerical': 0.010666822103054325, 'mean_word_length': 4.669695968274951, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Practical isogeometric shape optimization: parametrization by means of regularization


 Shape optimization based on isogeometric analysis (IGA) has gained popularity in recent years. Performing shape optimization directly over parameters defining the computer-aided design (CAD) geometry, such as the control points of a spline parametrization, opens up the prospect of seamless integration of a shape optimization step into the CAD workflow. One of the challenges when using IGA for shape optimization is that of maintaining a valid geometry parametrization of the interior of the domain during an optimization process, as the shape of the boundary is altered by an optimization algorithm. Existing methods impose constraints on the Jacobian of the parametrization, to guarantee that the parametrization remains valid. The number of such validity constraints quickly becomes intractably large, especially when 3D shape optimization problems are considered. An alternative, and arguably simpler, approach is to formulate the isogeometric shape optimization problem in terms of both the boundary and the interior control points. To ensure a geometric parametrization of sufficient quality, a regularization term, such as the Winslow functional, is added to the objective function of the shape optimization problem. We illustrate the performance of these methods on the optimal design problem of electromagnetic reflectors and compare their performance. Both methods are implemented for multipatch geometries, using the IGA library G+Smo and the optimization library Ipopt. We find that the second approach performs comparably to a state-of-the-art method with respect to both the quality of the found solutions and computational time, while its performance in our experience is more robust for coarse discretizations.


Introduction
Isogeometric analysis (IGA) introduced in [22] is a Galerkin method that uses splines to approximate both the geometric domain and the solutions to partial differential equations (PDEs). Splines are commonly used in computer aided design (CAD) and IGA is an attempt to bridge the gap between simulation and design [9]. This makes it beneficial for shape optimization as the optimization can be performed directly over parameters defining the CAD geometry, for example the control points of a spline parametrization, and it opens up the prospect of seamless integration of a shape optimization step into the CAD workflow.
One of the key challenges when using IGA in general, is that one needs a parametrization of the interior of the physical domain, on which the PDE is posed [42,20,17,32]. This parametrization 1 is used to pull back the weak form of the PDE to the parameter domain where the basis splines (B-splines) are defined. The choice of parametrization can affect the accuracy of the resulting IGA discretization [17,46] and, at the very least, the parametrization should be valid (a bijective map), that is, its Jacobian determinant should be non-zero. One approach to constructing a valid parametrization in 2D is to search for the one which has harmonic inverse. In [20] this property is reformulated as a nonlinear PDE and the parametrization is found by solving this PDE. In [19] this PDE based parametrization technique is used for a gradient based shape optimization algorithm with IGA. In [17,33] the same property is attained by minimizing the Winslow functional [45]. The method can be made more flexible by the use of adaptive splines [11], that allow to enrich the feasible region near to complex boundaries. In the recent works [42,39] the approach of parametrizing a complex domain by deforming a given template is explored. In [35] the focus is on producing parametrizations with low-rank with respect to the coefficient tensor.
When using IGA for shape optimization the challenge of finding a valid parametrization is even more important, since the shape of the physical domain changes during an optimization process. This means that a valid parametrization needs to be maintained during this process.
To guarantee that the parametrization remains valid during the optimization process, shape optimizations methods based on IGA often rely on constraints on the Jacobian of the parametrization [14,33]. These can be enforced either by using injectivity cones or by using the spline coefficients of the Jacobian determinant [47]. However the number of constraints needed quickly becomes very large, especially in 3D. Furthermore, when using the coefficients of the Jacobian determinant for the constraints, as we will do in this work, it may be necessary to expand the Jacobian determinant on a finer spline space, which increases the number of constraints even further.
In this work we will compare an existing approach to IGA shape optimization, relying upon reparametrizing the domain, with a simple approach to maintaining a valid parametrization without the use of explicit validity constraints. Namely, one lets the positions of all the control points that define the parametrization and the shape of the domain enter the formulation as independent optimization variables and adds a regularization term to drive the optimization towards a design with a valid parametrization. Such an approach has been considered in the context of shape optimization in mechanics in [40], and has to the best of our knowledge, only been considered very briefly in the context of shape optimization with IGA in [13]. It remains a question whether this approach performs comparably to state of the art methods, and the aim of the work is to investigate exactly this question. Different frameworks for shape optimization with IGA exist in the literature. For example in [29] a tetrahedral mesh is used to represent the interior of the computational domain, while still representing the boundary using splines. The authors rely on external mesh generator software, while we in this work aim to avoid mesh generation by instead maintaining a parametrization of the interior. In [15] the shape is represented implicitly as the level set of a function. The PDE is then posed on a design domain, which remains constant during the optimization, and the level set function enters the formulation by the ersatz material approach. This framework has the advantage that it allows for changes in the topology of the shape, however at the same time the final shape is represented as a level set and therefore post processing is required to represent this shape using splines, which is necessary for importing the shape into CAD software. In our work the optimization is performed directly on the spline representation and thus the result can be readily imported into CAD software after the optimization. Another possible approach is to use the isogeometric representation of geometry combined with the boundary element method (IGABEM) [28,6]. In IGABEM the PDE is reduced to an integral formulation on the boundary of the domain, and therefore maintaining a parametrization of the interior of the domain is avoided altogether. However in this case one has to deal with the standard complications of BEM, in particular the fact that the system matrices are dense, non-symmetric, and costly to compute. Additionally the Greens function for the considered PDE has to be known, which is not always the case.
In this work we illustrate that the simple regularization based approach is able to handle complicated geometries, by comparing its performance to a shape optimization approach based on using a linearization of the Winslow minimization problem as a parametrization strategy, and employing locally refined splines to represent the Jacobian determinant. The method closely resembles the one in [33]. The main difference is that we use Truncated Hierarchical Basis splines (THB-splines) [16], which possess the partition of unity property, to refine the spline space in which the determinant is expanded locally. This reduces the number of constraints needed compared with tensor product global refinement.
We will apply the two methods and compare their performance on the shape optimization problem of designing electromagnetic reflectors. In this problem we have two metallic reflectors in a dielectric medium and search for a shape that maximizes the electrical energy close to a chosen point. The same problem has been studied with topology optimization in [1,44,7], and with IGA in [33].
The methods are implemented for multipatch geometries, using the IGA library G+Smo 1 and the optimization library Ipopt 2 . A benchmark study on the performance of different optimization algorithms in the context of structural optimization can be found in [38]. For a review of implementation aspects of IGA in general we refer the reader to [34] and for details of the G+Smo library we refer the reader to [26,30]. The code used for this work can be found at https://github.com/gismo/shapeopt.
The paper is organized as follows. In section 2 we outline the relevant notation, and in sections 3 and 4 we describe the two methods we are going to compare. In section 5 we apply these methods to the aforementioned shape optimization problem, and discuss and compare the performance of the two approaches. We end the paper with some discussion and conclusions. Some of the more technical details are presented in the appendices.

Preliminaries and notation
Let us consider the following PDE-constrained shape optimization problem: where O ad is a set of admissible shapes, E is the objective and (1b) is the governing PDE in the weak form.
Within the IGA framework both Ω and u will be approximated numerically using splines. Namely, we have Ω = G(]0, 1[ d ), with the parametrization where c i ∈ R d are the control points, N g is the number of control points, and R g i are the basis (B) splines. In this work, unless specifically stated otherwise, we will use tensor product B-splines. The superscript g indicates that the B-splines R g i (that is, their degrees and knotvectors) are specific to the geometry representation. Within the shape optimization framework it will sometimes be necessary to distinguish between boundary and inner control points. We will therefore introduce the notation where c b are the boundary control points and c i are the inner control points. The Jacobian J = ∂G ∂c will also play an important rôle in the forthcoming development.
Similarly to (2), we approximate the state of our system as a pulled back spline where R g i , i = 1, . . . , N are B-splines. The expansion coefficients u i , i = 1, . . . , N will be found by solving a system of linear algebraic equations As standard in the Galerkin approach, the elements of the stiffness matrix K c and the load vector . Note that the dependence of K c , f c on the control points is encapsulated in (2) and (4).
Already at this point the importance of the geometry parametrization should be apparent. Indeed, at the very least it should be an invertible map, which is used to pull back the weak form of the PDE defined on the physical domain Ω into the parameter domain ]0, 1[ d . In particular, for all ξ ∈]0, 1[ d it is necessary that det(J(ξ)) > 0. 3 A sufficient condition, which guarantees the validity of the parametrization, is discussed in Appendix A.
3 Boundary-driven approach to IGA shape optimization In this section we will consider one possible approach to IGA shape optimization, which follows the ideas developed in [33,17]. Within this framework we formulate the optimization problem in terms of boundary control points c b . The collection of interior control points c i for the geometry parametrization is treated as an implicit function of c b , see Appendix B. Additionally, the parametrization validity constraints det(J) > 0, or a sufficient condition for these (cf. Appendix A), have to be explicitly included into the problem formulation.
In order to compute domain parametrizations of high quality we rely upon minimizing the Winslow functional (Section B.2). However, to avoid solving a non-linear optimization problem at each shape optimization iteration, we construct a quadratic approximation to the Winslow functional around a reference parametrization, and update the reference parametrization when it becomes necessary to do so. Specifically, given a reference parametrization G 0 defined by the control points c 0 , to find a new parametrization we consider the quadratic programming problem where W is the Winslow functional, and H is its Hessian. The minimizer of this problem can be found by solving a linear system H(c 0 )∆c = −∇W (c 0 ).

Using (3), this can be restated as
where The new parametrization is then defined by the control points given by c = c 0 + ∆c.
With this in mind, to approximate (1) numerically we solve a sequence of subproblems where c 0 is the reference parametrization, ∆c = (∆c b , ∆c i ), and (7e) is the sufficient condition for the validity of the parametrization discussed in Appendix A. Each time we solve the subproblem (7) we update the reference parametrization. In the model problem considered in this work we saw no further progress after 5 to 10 reparametrizations. Reference parametrizations can be computed as follows. We minimize the Winslow functional as described in Appendix B, and check if the sufficient condition d > 0 is violated. If it is, then this condition is too strict and should be relaxed. To facilitate this we refine the spline space S det where we compute expansion coefficients d of det J. To reduce the number of constraints resulting from such refinement steps, we utilize local refinement. Specifically, we use Truncated Hierarchical B-splines (THB-splines) as basis functions. Note that it is important here to use the truncated version of hierarchical splines, since the partition of unity property (cf. [16]) implies that the spline control polygon converges locally to function values. Note that other locally refinable splines with this property are available, for example Polynomial splines over Hierarchical T-meshes (PHTsplines) [10] or Locally Refinable (LR B-splines) [24]. For a comparison of these methodologies we refer the reader to [27].
The refinement strategy we employ is as follows. For all indices i that have negative spline expansion coefficient d i ≤ 0 of det J, we refine the support of the associated basis function R det i . This is repeated until d i > 0 for all i = 1, . . . , N det . 4 In the subproblem (7) we then put ε = ρ·min i d i with ρ = 0.25, see (11).
The full optimization loop is illustrated in Figure 1. Note that we have chosen to carry out the spline space refinement described above only when the reference parametrization, and therefore also the subproblem (7), is updated. This allows us to keep the number of constraints constant when solving (7) numerically and therefore employ off-the-shelf optimization software. The initial guess for the non-linear Winslow optimization problem is generated using Coons' patches, see Appendix B.
4 Regularization-driven approach to IGA shape optimization In this section we discuss an alternative approach to shape optimization using IGA, which does not involve explicit constraints on det J. The positions of the inner control points c i enter this formulation as independent optimization variables, in the same way as c b . Consequently, we do not need to explicitly compute a domain parametrization, as this will be part of the outcome of the optimization process.
To this end we add the Winslow functional W as a regularization term to the objective function. Its role is to penalize configurations of control points that result in poor parametrizations. This Solve min c i W to find reference parametrization c 0 .
Compute d.
Did we converge? Stop.

Yes
No Figure 1: Flowchart of the optimization algorithm.
idea has been used previously in the context of shape optimization in mechanics [13,40]. Thus for a regularization parameter τ > 0 we consider the optimization problem We put W (c) = ∞ if det J ≤ 0 at one of the quadrature points used for the integration when calculating W (c). In this way, when the optimization algorithm enforces e.g. the standard sufficient decrease condition, such a choice ensures that the chosen step will always have positive determinant at the quadrature points. This does not guarantee that it is positive everywhere, but it means that the numerics will not collapse due to a division by zero. The regularization parameter τ needs to be tuned for the specific problem at hand. If it is too large the minimization will find a design with a small value of the Winslow functional but disregarding the objective E(c, u). If it is too small the optimization will find positions of the control points that have a low objective E(c, u), but with a poor parametrization, which might give a large discretization error of the discretized PDE. The appropriate values of τ would lead to a compromise between these two extreme situations. One simple strategy for choosing such a value is to solve a sequence of problems (8) for decreasing values of τ . 5 The two points above constitute the main drawback of this method. Namely that, in contrast to a constraint based method, we cannot guarantee that the parametrization is valid in between the quadrature points. Additionally one has to find a suitable value of the regularization parameter τ , which is problem-dependent.

Reflector (Gold)
Air Point with maximum field Electromagnetic wave Truncated Domain ? ? Figure 2: Sketch of the shape optimization problem. The goal is to find a shape of the reflector that maximizes the field close to a point.

Case study: Optimization of electromagnetic reflectors
In this section we will consider a 2D shape optimization problem originating from the field of electromagnetism. Our goal is to design a reflector that concentrates electrical energy in a desired area. This problem will serve as a model problem for comparing the two optimization approaches outlined in Sections 3 and 4.

Physical model
We consider a two dimensional scattering problem where a plane wave with frequency f travels in a dielectric (air) and is scattered by two symmetric metallic (gold) reflectors, as depicted in Figure 2. Let cr and µ r denote the complex permittivity and permeability of the medium. Using the first order absorbing boundary condition [23] at the boundary Γ t of the truncated domain, the electromagnetic fieldû should satisfy the following PDE: In the equations above, k 0 = 2π √ 0 µ 0 is the wave number and 0 , µ 0 refer to the permittivity and permeability, respectively, of free space. The imaginary unit is denoted by j, the radius of the truncated domain is given by r t , and u i is the incident plane wave, given by The objective function of the shape optimization will be given by where δ is a Gaussian bell-function   Figure 3: Patch layout.
with α = 0.1. Thus we aim to focus the incoming energy in the vicinity of the origin (0, 0). The physical parameters that we use are given in Table 1. The complex permitivity of the reflector is calculated as s cr = r,gold − j σ ω 0 . The weak statement of the PDE (9) is to findû ∈ H 1 (Ω) such that for all test functionŝ v ∈ H 1 (Ω) the following equality holds Due to the symmetry we only consider the upper half of the geometry shown in Figure 2. To accommodate the change of material parameters between the metallic reflector and the surrounding dielectric medium, we will split the domain into five patches, one for the reflector and the other four for the surrounding air. The layout is shown in Figure 3. 6 Each patch is parametrized using splines of degree p = 2 as described in Section 2. We will use strong patch coupling to enforce C 0 continuity at the patch interfaces, while noting that other alternatives exists, for example weak coupling [21] and strong C 1 coupling [5].
The shape of the reflector is represented via the patch interfaces, which consist of 4 spline curves with C 0 continuity at the 4 corners. Such a description, for example, allows shapes that are only piecewise smooth such as e.g. the classical bowtie antenna [8]. There are other applications of IGA where even more complex geometries are relevant, for example when modelling crack propagation as in [37,36].
Using the parametrizations we can pull back the equation (10) to the parameter domain and apply the Galerkin method to it, which ultimately results in a system of linear algebraic equations; see Appendix C for details. To evaluate the integrals involved we use element-wise Gauss-Legendre quadrature, owing to their immediate availability in the G+Smo library, while noting that more efficient alternatives exists [4,2].

Results with boundary-driven approach
In this section we apply the method described in Section 3 to our model problem. We start with an initial design where the reflector has the shape of a circle. We consider two different spline spaces in which to approximate the PDE (9), namely where the knotvectors used for representing the geometry are refined uniformly either 3 or 4 times. 7 Both spline spaces have degree p = 2, and the number of degrees of freedoms are N coarse = 2548 and N fine = 9300, respectively. We will refer to these as the coarse and fine meshes.
We will use a tolerance tol = 10 −3 when solving the subproblems (7) and a fixed number of reparametrizations, namely 10 when using the coarse mesh and 5 when using the fine mesh. We observed that using more reparametrizations did not lead to significant improvements in the design. In our experiments the results with this method are sensitive to the number of quadrature points used when calculating the Winslow functional. To produce the results presented here we use 12 quadrature points per knot interval, to avoid under-integration 8 .
In our implementation we use the interior point solver Ipopt to solve the subproblems (7). As we solve a sequence of subproblems a warm start is available, namely the minimizer from the previous subproblem; so an interior point algorithm might not be the best choice of optimization algorithm [25]. However with the following parameter tuning we found the method to work well with Ipopt. Namely, after solving one of the subproblems (7) some of the design bounds (7f) will be active; however, since Ipopt is an interior point algorithm, the starting point for the subsequent subproblem will be pushed away from the boundary as controlled by the parameter bound push. We found that this parameter needs to be lower than the default value since the constraints on det J are quite sensitive and a relatively small perturbation of the control points might violate these constraints, which is undesirable. The default value is 0.1, but in the experiments we set it to 10 −5 instead.
Another key parameter is the barrier parameter mu init. Specifically, we use the monotone strategy, where the barrier parameter is monotonically decreased as the optimization algorithm progresses. However, if this parameter is too large in the beginning of the algorithm we found that it will push the design towards configurations with large det J. To remedy this, instead of the default value 0.1 we use 10 −4 . For more information about the optimization algorithm implemented in Ipopt and its parameters see [43].
In Figure 4 the design at different stages of the optimization is presented, when using the coarse mesh.
We observe that the design becomes increasingly hard to parametrize as the objective increases each time we change the reference parametrization. In Figure 5 the objective is plotted against the number of iterations. We see that the objective function increases after the reference parametrization is changed, but relatively quickly reaches a plateau. We already use a fairly large tolerance of 10 −3 for the stopping criterion when solving the subproblems (7), however this behaviour indicates that it might help to relax the stopping criterion even further in these subproblems to improve the overall efficiency of the method. However, to allow for a fair comparison between the two methods, we do not investigate this further and use the same tolerance for both methods. The final objective, after 10 reparametrizations when using the coarse mesh, is E h = 1.803. However if we calculate the objective with a mesh that is refined uniformly twice we get E h/4 = 1.556, that is, a 16% difference. This means that there is actually a significant discretization error at this refinement level, which leads to an artificially large objective value.
As described in Section 3 we use local refinement to adapt the constraints on det J to the current design. This is done by using local refinements in the areas where the spline coefficients are non positive. The resulting meshes are plotted in Figure 6 to illustrate where the refinement is needed. We see that it is primarily inside the reflector and near the reflector-air interface, that this  Figure 4: The designs at different stages of the optimization process, when using validity constraints and the coarse mesh. The reflector is outlined with a black line, and the control points of this boundary is colored black. The grey lines are parameter lines mapped with the geometry map, to illustrate the parametrization.  (d) After 10 reparametrizations, N det = 572. Figure 6: The mesh (knot lines), used for representing det J during the subproblems (7) for the constraints (11). The number of constraints is given by N det and the designs are obtained using the coarse mesh. Figure 7: The final design using validity constraints and the fine mesh, E h/2 = 1.638 refinement is applied. We see that the number of constraints increases during the optimization process, but by no more than a factor of two. If we were to use uniform refinement the number of constraints would increase by more than a factor of 9.
The final design, after 5 reparametrizations, when using the fine mesh, is shown in Figure 7. We see that it is similar to the design obtained using the coarse mesh. The final objective here is E h = 1.638. After the mesh is refined uniformly the electrical energy is E h/4 = 1.628, which is only a 0.6% difference. The objective function during the optimization process is plotted in Figure 5. We see that the algorithm converges faster when using the fine mesh as no progress was observed after 5 reparametrizations.

Results with regularization-driven approach
In this section we will present the results obtained with the method described in Section 4. With this method we perform the optimization with all control points as optimization variables while using the Winslow functional as a regularization term. Again we use the interior point solver Ipopt with a tolerance of tol = 10 −3 for solving the problem (8). The regularization parameter is set to τ = 1 8 . One can compare the design evolution shown in Figure 8 with those obtained previously, see Figure 4. The designs obtained using the regularization approach seem to have more regular parametrizations compared to those in Figure 4.
The final objective is E h = 1.684. If we calculate the electrical energy for this design on a twice refined mesh we get E h/4 = 1.546, that is, a difference of 9% . This is less than the 16% we observed when using the boundary-driven method. This increase in accuracy might be due to the parametrization being of higher quality.
In Figure 9 we plot the electric energy E h , the regularization term τ W and the objective function τ W − E h .
When using the method with the fine mesh we get the final design shown in Figure 10. We see that the shape of the reflector is very similar to the one shown in Figure 8d. The main difference is that when using the fine mesh the parametrization is more regular, since the error in the discretization of the PDE is smaller, and therefore the optimization cannot exploit it to the same extent. This is especially notable at the bottom of the reflector were the inner control points where moved away from the point of interest when using the coarse mesh, as seen in Figure 8d   result E h/4 = 1.545 with the difference at the 5th digit. The objective during the optimization process is plotted in Figure 9. The behaviour is similar for the two meshes, however the tolerance tol = 10 −3 is reached with fewer iterations when using the fine mesh.

Comparison and discussion
In Table 2 we summarize the performance of the two methods. We report the objective after a fixed number of iterations, in this case after 100 iterations, the final objective computed on 3 different refinement levels and the average running time per iteration. We observe that the average execution time per iteration is the same order of magnitude for the two methods. The main portion of the running times is spent on solving the state equation and computing the gradient of the objective function. The difference in the running time that we observed between the two methods might be due to a different number of function evaluations per iteration needed for trial steps of the algorithm. If we were to consider a larger problem, for example in 3D, the large number of validity constraints would likely lead to an increase in running time for the boundary-driven approach. For the boundary-driven approach the computation of reference parametrizations by minimizing of the Winslow functional only accounts for 2.4% and 0.25% of the total running time for the coarse and the fine grid respectively.
Regarding the quality of the designs we find, we note that shape optimization problems are prone to having many local optima, so it could be that the two methods find two different local optima. Therefore it can be futile to directly compare objective values. That being said, we observe that the boundary-driven approach happens to find solutions with slightly higher (better) objective value. On the other hand, the regularization based approach seems to estimate the objective value more accurately on coarser meshes, since our results were more reliable using this method, probably due to the better quality of the parametrizations it produced.
Also, observe that with the regularization based approach we only need the objective E h , the Winslow functional and their first order derivatives. In addition we can solve a single optimization problem with design bounds as the only constraints. This means that the method is significantly easier to implement.

Conclusion
We described and compared two methods for shape optimization on spline-based representations. One uses validity constraints to enforce the validity of the geometry parametrization. The other uses a regularization term, and thus avoids both the validity constraints and the need for an explicit parametrization strategy altogether. We demonstrated how this simple approach performed comparably to the more complicated approach in terms of the final design, while requiring similar running times for the 2D problem we considered. The regularization based approach seems to produce more reliable results and it is in addition much simpler to implement, since we only need the objective, the Winslow functional and their first order derivatives. These results are encouraging and we plan to use the regularization based approach for 3D problems, where we expect that the efficiency advantages of the regularization approach will be more prominent, since the number of validity constraints for 3D parametrizations grows quickly. For obtaining competent overall running times, one can couple the approach with low-rank approximation techniques, which have been recently introduced in IGA [31,41] and have proved efficient in the frame of PDE-constrained optimization [3]. Now we can use this expansion to derive a sufficient condition to replace the pointwise parametrization validity constraints det(J) > 0, by requiring that where ε is a small positive algorithmic parameter. This condition guarantees that det J > 0 since B-splines are non-negative and form a partition of unity. This is not a necessary condition, meaning that we might have d i ≤ 0 for some i even though det J(ξ) > 0, ∀ξ ∈]0, 1[ d . But if the spline space S det is refined then the spline expansion coefficients will move closer to values of the spline. So if ε is small enough and det J > 0 then the constraint is likely to be satisfied for a sufficiently refined spline space S det .

B Domain parameterization techniques
In this section we review some techniques for finding a parametrization of the interior given the boundary. In IGA this comes down to finding the position of the inner control points given boundary control points. We do not aim to give a thorough review of all the techniques that are available, as this is out of the scope of this work. We will only introduce the methods that are related, or directly used, by the two shape optimization methods considered in this work

B.1 Coons patch and spring method
Two simple methods for constructing a grid of control points are the Coons patch [12] and the spring method [17]. They both produce inner control points that depend linearly on the boundary control points, however they only produce valid parametrizations for geometries that are not too complicated. In this work we will use these methods for finding an initial guess for the optimization based approach that is described below.

B.2 Optimization-based techniques
A more complex, and in general nonlinear, class of parametrization methods consists of optimizationbased methods. Here the geometry map G is chosen such that it minimizes a quality metric w(ξ) where γ is a given boundary curve.
There are several different quality metrics to choose from, cf. [17,48]. In this work we will consider the Winslow functional which is given in terms of the Jacobian matrix J as with w = tr J T J det J .
In 2D the Winslow functional has the nice property that its minimizer has a harmonic inverse [17]. This guarantees that the minimizer of (12) is unique and bijective, i.e. det J = 0. It should be noted that this minimizer is not necessarily a spline, so looking for a spline parametrization on the form (2) by minimizing (12) with the Winslow functional will only guarantee a valid parametrization if the spline space used for the parametrization has high enough resolution.