Realization of a Framework for Simulation-Based Large-Scale Shape Optimization Using Vertex Morphing

There is a signiﬁcant tendency in the industry for automation of the engineering design process. This requires the capability of analyzing an existing design and proposing or ideally generating an optimal design using numerical optimization. In this context, efﬁcient and robust realization of such a framework for numerical shape optimization is of prime importance. Another requirement of such a framework is modularity, such that the shape optimization can involve different physics. This requires that different physics solvers should be handled in black-box nature. The current contribution discusses the conceptualization and applications of a general framework for numerical shape optimization using the vertex morphing parametrization technique. We deal with both 2D and 3D shape optimization problems, of which 3D problems usually tend to be expensive and are candidates for special attention in terms of efﬁcient and high-performance computing. The paper demonstrates the different aspects of the framework, together with the challenges in realizing them. Several numerical exam-Communicated by Zenon Mróz.


Introduction
Design optimization has been an indispensable part of the engineering design process across many disciplines; prominent fields include aerospace and airplane design. In this context, topology optimization has been used extensively to generate optimal designs under given conditions [2,5,39]. With the developments in recent times, shape optimization has also gained importance in simulation driven design and optimizations. [4,15,27,29] present some interesting applications of shape optimization. Traditionally in shape optimization, a parametrization of the shape of the object and corresponding parameters are used in optimization as design variables. For complex geometries, this parametrization procedure is tedious and can limit the freedom in shape optimization. In contrast to the above node-based shape optimization provides a higher degree of freedom, by taking the nodal coordinates of the discretized geometry as design variables.
The advantage of using the nodal coordinates as design parameters is twofold; first, it provides the largest possible design space for optimization, and second, it eliminates the necessity of using an explicit parametrization of the geometry. Previous works together with the original paper of [8] and the works [8,35] introducing vertex morphing regularization in the context of node-based shape optimization have shown promising results in different applications. Today, thanks to the improvements in the computational power and hardware, the number of design variables used in shape optimization have increased by many folds. This increase, along with an expanding number of application areas, is pushing the performance limits of the optimization techniques and the frameworks providing optimization and related capabilities.
In the wake of increasing usage and importance of node-based parametrization, an efficient and flexible framework offering respective functionalities becomes the need of the hour. The two essential characteristics of such a framework are its ability to include new algorithms and techniques, and the ability to extend the existing ones. Apart from these, it is vital that the framework is able to work with external solvers which provide the primal solution and sensitivity information required for shape optimization.
In this work, we present a variety of multi-disciplinary applications of shape optimization achieved with the framework along with a description of its important capabilities. The above-listed requirements of the framework are realized using implementations in C++, which then are exposed to Python. The implementations in C++ deal with computationally expensive tasks, whereas the python routines steer the optimization process. The python scripting layer gives the necessary flexibility and ability to prototype different optimization problems and algorithms rapidly.
Previous works on similar tools for a general optimization like pyOpt [31], GEMS [17], DAKOTA [14] though have several optimization algorithms readily available, do not provide functionalities necessary for computing geometrical constraints, which arise in shape optimization problems. Apart from this, the lack of possibility to deal with different external solvers and parametrizations in a unified way makes them less suitable for shape optimization. Recently developed tools put forth like ABCD tool box [32] come close to the current work, but are developed with focus on specific applications and especially lack capabilities to handle node-based shape optimization. Though OpenMDAO [21] frameworks' objectives are parallel with the current work, it relies on applications and tools with extensive access to inner data structures to compute necessary quantities for optimization. modeFRONTIER from ESTECO is a commercial tool with a user-friendly GUI platform which has the similar capabilities as the current work. Ebes + jMetal a tool presented in [38] though can do multi-objective optimization, lacks capability to customize and specialize the provided objectives or constraints and depends heavily on external tools for this. On the other hand, tools for solving coupled multi-physics problems like EMPIRE [36] and preCICE [10] contain the ability to work with different solvers. However, their focus is on a different set of problems and are less suitable for advanced and sophisticated shape optimization problems.
The current framework used in the current work differs from the tools mentioned above in terms of its focus on requirements for shape optimization using vertex morphing parametrization. It also provides a flexible and generic interface to compute and include geometrical and other special types of constraints. It also uses a flexible methodology for the exchange of data required for optimization. These features also allow the usage of the current framework together with other tools, thus enabling a multi-physics optimization; more details are presented in Sect. 4.3.1.
The remainder of the paper first presents the theory of shape optimization and vertex morphing technology used in the optimization. The discussion of software concept and the details framework follows the theory. The paper concludes with the discussion of different application cases and results obtained using the framework presented, thus demonstrating its capabilities.
where ⊂ ∂ is the design boundary of the computational domain ; unless noted otherwise, X , X ∈ R 3 denotes the vector of nodal coordinates of size 3m ; A is a generic mapping matrix which performs transformation between the design geometry X and the design variables s; J , g and h, respectively, denote the objective function to be minimized, the vector of inequality constraints and the vector of equality constraints.
g and h are the subsets of the boundary which are subject to geometric inequality constraints g and equality constraints h , respectively. r denotes the residual of the state/primal governing equations which may be nonlinear and Q denotes the state/primal vector.
The geometric constraints can be enforced on the optimization problem in multiple ways. A direct approach is to formulate them point-wise, which means for each mesh point X , j , scalar-valued constraints h j and g j are evaluated and satisfied individually. In context of node-based shape optimization with vertex morphing, this has been investigated by [4] for no penetration constraints. As an alternative, multiple geometric constraints of the same type can also be aggregated into one single constraint and therefore be treated in a weak sense. In such cases, they will not be fulfilled on every point, but in an average sense. This approach is used in the example presented in Sect. 4.2.1. A third option is to directly include the geometric constraint into the design parametrization, and apriori exclude infeasible designs [34]. A similar approach is used in the example from Sect. 4.1. In addition to these geometric constraints, the optimization problem can be subjected to any number of physics-based equality and inequality constraints. Literature provides robust methods for enforcing these constraints [7,11,16].

Vertex Morphing Technique
Gradient-based first-order optimization methods can be used to solve the node-based shape optimization problem presented above. Apart from being one of the fastest methods, they are also robust with minimal effort. Adapting these methods and considering that for the node-based approach, the solution of optimization problem requires the gradients, also called sensitivities, of the objectives and constraints with respect to the design variables s. Since the nodal coordinates are considered as the design variables, the sensitivity with respect to the nodal coordinates at each node on the design boundary is required. For many classical objectives considered in engineering applications, these nodal sensitivities tend to be noisy and an optimization with noisy sensitivities produce an undesirable shape of the geometry, though it is optimal. Figure 1 shows and example noisy adjoint sensitivities calculated for strain energy and the corresponding optimal shape generated. To deal with this problem, in the current framework for optimization we use the vertex morphing parametrization technique proposed in [8] to consistently transfer the shape optimization problem from the spatial design space X (the geometry space) to a new design space s , called control space, using a mapping operation. This will also filter the sensitivity field, thus smoothing it. Unless otherwise noted, in the following, the control field lives on the same discretization as the surface of the geometry. This allows the largest design space possible with minimal As defined in Eq. 1, the association of the discretized geometry with the discretized design space is expressed via the operator matrix A. This operator maps s onto X as follows: In a continuous space, the three-dimensional geometry at point X 0 = (X 1 0 , X 2 0 , X 3 0 ) of the optimization surface is generated from the surface control field s = (s 1 , s 2 , s 3 ) via a smoothing filter operation: where F could be any reasonable filter (kernel) function. is the portion of which lies within a sphere of radius r centered at X 0 , where r is the filter radius (assumed to be constant); X − X 0 is the Euclidean distance to the center of the filter X 0 ; n is the unit normal vector to the surface; d = nd is the unit normal component of the surface element. See Fig. 2 for a schematic of the used notation. For more detailed information on different filter functions and their effect on the nature of the A matrix, the readers are referred to [8].
As discussed, the key component required for any gradient-based shape optimization method is the objective and constraint function gradients with respect to the design variables d( J,g,h) ds which will be used to generate the geometry update s. In case, the derivatives of the functions J, g, h are computed w.r.t to the spatial coordinates dX 0 , it is required to compute the shape derivative, i.e., the derivative of the surface coordinates w.r.t. control field dX 0 ds , see Eq. 7. Various established methods for calculating the sensitivities are in practice; for example [5] shows a discrete adjoint method, [6] introduces a semi-analytical method, and [19] among others showcases a continuous adjoint method this adjoint sensitivity calculation is extended to unsteady problems in [13]. Some of these methods are already available for usage in commercial and open-source software packages. An in-depth discussion about the available methodologies for calculation of sensitivities is out of scope for the current contribution. Since the current contribution focuses on different aspects the framework necessary for shape optimization and the applications themselves, in the following discussion, we consider the shape sensitivities are given and the computational costs involved are not considered though they cannot be neglected. In the examples presented in Sect. 4, the methods used for calculating the sensitivities are explicitly mentioned for each application as multiple methods are used. This in turn establishes the versatility of the framework much strongly.
Given the control field and the objective gradients at the k th optimization iteration, the new set of design variables can be calculated as: where s k is the update of control field. For a simple steepest descent algorithm, the computation of the new design can be written as: Here α is the optimization step length and ∇ s k is the gradient of the objective J with respect to the control variables s. Other more advanced algorithms like sequential quadratic programming [20], front propagating algorithms [28] can also be used to generate the value of s k . In the current contribution, unless otherwise mentioned, a steepest descent algorithm is used to generate the update of the control field. ∇ s k can be calculated using the chain rule of differentiation. Thus, the objective derivative with respect to the shape controls reads: where d J dX k and dX k ds k are the objective derivative with respect to mesh coordinates and shape derivatives of vertex morphing, respectively. From Eq. 4, the value of dX k ds k can be given by the expression [4,8,9,24]: Following the above derivations, the shape update can be computed as: The mapping operation in Eq. 9 ensures the smoothness of the shape update calculated. Finally, it is worth to mention that vertex morphing can preserve the main features of the design surface. More precisely, it allows shape changes which are not affecting the "design character" aesthetic and geometrical features [8,23]. As mentioned previously, F and r are the design handles which control the feature preservation. A good example is the preservation of feature lines (sharp edges), which could be achieved by a proper choice of the filter radius r [23]. The reason for this is that all the features smaller than the radius r are only subject to a bulk and rigid motion, without considerable shape deformation. Another example would be the thickness control which results in the no self-penetration property. For applying the constraints defined by g, h their derivatives d(g,h) dX 0 are transferred to the control space by the following operation dX k (10) and the constraints can be applied using one of the many methods available. In this work, a version of the gradient projection method introduced in [22, Chapter 5] is used to apply the constraints. This method satisfies constraints by projecting the descent direction onto the subspace tangent to the active constraints. Furthermore, a correction is performed to bring the design update back to the feasible domain. As a conclusion of this section, Algorithm 1 outlines the steps followed and computations required for the constrained node-based shape optimization by means of vertex morphing.

Software Framework
In the fast-evolving design process involving teams from different backgrounds, the ability to test various ideas with several constraints and objectives is fundamental. For optimization, this means the ability of the framework to achieve this with minimal effort from the user. For efficiently solving the optimization problem described in Algorithm 1, the following vital capabilities of the framework are noteworthy: Usage of external tools and solvers as black-box In the optimization workflow described in Algorithm 1, the objective value and its gradient (statement 6), constraint values and their gradients (statement 8) with respect to the shape are calculated by the external solvers and provided to the optimization framework. This enables a modular and black-box treatment of the solvers for evaluating the responses and if necessary constraints. This approach also allows the framework to perform optimization based on different physics and design parameters. Based on this rationale, a black-box treatment of the solvers is used in this work. For such a design, the communication module with different external solvers to transfer the quantities as mentioned earlier becomes a vital part of the optimization framework.

Parametrization techniques
In shape optimization, parametrization of geometry is the procedure to define a set of parameters to control the shape that is optimized. Choice of a the parametrization can result either in a reduction or an increase in the number of design variables, thereby effecting the computational effort. Different parametrization techniques have been previously proposed and used [33]. As not every technique is suitable, efficient, and possible for every optimization problem, the optimization framework should have the ability to extend and assimilate different parametrization techniques quickly. The design structure of the framework allows for different parametrizations. The current work implements a discrete node-based parametrization approach described in Sect. 2.1. The versatility of the presented framework is also demonstrated in its ability to nest between two different parametrizations during the optimization processes. An example is presented in Sect. 4.1.
Computing geometric and other specialized responses As mentioned previously, the external physics solvers perform the calculation of state-based responses, constraints, and their gradients. In addition to such responses, other types like geometric [4], stamping constraints presented in Sect. 4.2.1 are also crucial in industrial applications of shape optimization. Computation of some of these specialized responses may not be possible with the external solver. In such cases, the possibility to compute them within the framework is necessary.
These features are critical in enabling the researchers to experiment with algorithmic procedures and allow the usage of the framework for industrial applications by practitioners.

Implementation
The implementation of the optimization framework is conceptually divided into two parts: first, the computation and memory intensive part, and second is the algorithmic part which organizes the data and workflow. Considering the nature of these two parts, C++ and Python are chosen, respectively, for them. This combination simultaneously enables efficient computations and improved readability of the algorithmic procedure. The Python layer also helps with fast and easy prototyping and setup of new and experimental optimization procedures. Figure 3 shows the thematic arrangement of different components of the framework and steps in the optimization workflow. It also illustrates the critical aspects of the framework, of which the distinction between Python and C++ layers and the interface with the external solvers are noteworthy. The C++ implementation deals with computationally intensive tasks of formulating the mapping between different parametrizations (Eq. 1), node-search, linear algebra, and the data structure. The second part, implemented in Python, deals with orchestrating the functionalities in an optimization algorithm.

Parameterizations and mappers
In the vertex morphing technique described in Sect. 2, though the design parameters are the nodal coordinates of the discretized surface, the optimization happens on a different control space s. A transformation between the geometry coordinates and control parameters is performed using a mapping matrix (Eq. 8). This same procedure can be used for transformation between arbitrary parameterizations. Such generalization allows the possibility of having different parameterizations. Section 4.1 presents a series of nested parameterizations using this generalization.

Communication with black-box solvers
One of the essential features of a general optimization framework is its ability to work with solvers treating them as "blackbox" and communicate the objective and constraint values, their gradients and the shape updates. Such black-box treatment of solvers enables the framework to work on different optimization problems with different physics involved. A novel detachedinterface approach enabling modularity with minimal and framework independent changes to the solver is developed. The features of this approach are: -The interaction between the framework and the external "black-box" solver(s) is done via a solver-wrapper developed to adhere to the interface of the base solver in the framework. This allows uniform treatment of different solvers-wrappers and easy switch between them. Together with functions for exchanging data, this interface also contains functions to control the solver. A solver-wrapper will only implement the necessary functions depending the degree of control the "black-box" solver allows. This wrapper can also contain specific routines for communicating  Figure 4 shows the important functions in the interface. -The actual data communication between the framework is done via a different input/output (IO) class object. The functions from Fig. 4, ImportDataField, ExportDataField, ImportModel, ExportModel are responsible for the data exchange in the solver-wrapper. These functions are then delegated to the IO class object. This delegation of data exchange to an interfaced IO class enables decoupling the IO from the solver-wrapper and thus enables re-usage of existing IO methodologies and implementations between different solver-wrappers. -The external "black-box" solver is completely independent from the developed framework. This implies that the solver can choose and implement any routines to export the objective or constraint sensitivities together with their values. This will enable reuse of already existing routines developed for other tools, thus simplifying the deployment of the solvers on different computing environments.
The following solvers-wrappers for external solvers are implemented in the framework using the above-described detached interface approach and are readily available for usage.
• OpenFOAM [37] • SU2 [30] • KratosMultiphysics [12,25] • Altair OptiStruct The usage of the above solver-wrappers is also stated in the numerical examples presented in Sect. 4. The following are the steps to setup a new solver-wrapper, using the above-described detached interface, and use a new black-box solver together with the optimization framework: -Implement the following routines in the solver -To output the objectives and their sensitivities.
-To accept the design update from the optimization framework.
-To accept a convergence signal from the optimization framework.
This implementation in the solver has no dependency on the optimization framework whatsoever. Thus, the developers are free to use any tool of their choice and convenience to achieve this functionality.
-Use the interface provided by the optimization framework to create a wrapper to the solver either in C++ or in Python. This wrapper will be responsible for: -Interacting with the solver and give necessary instructions to the solver to perform different computations. For example, to compute the primal or sensitivity values. -Provide the design update computed by the optimization framework to the solver.
Though these operations depend upon the data structure of the optimization framework, these can be customized to the needs specific to the solver. For example, data transfer via file I/O or in memory transfer using MPI protocols.
The geometric constraints and the procedures presented in Sect. 4 are implemented following the procedure mentioned above and are included in the framework.

Optimal Shape Design Applications and Numerical Results
In the following section, we present different application cases which are optimized using the above-described optimization framework. Since the objective of the presented examples is to showcase the capabilities of the framework, description of the specialized optimization procedures adapted is provided when possible and references are provided otherwise. Unless explicitly mentioned, all the simulations below have a numerically stable behavior. The application examples are grouped in different focus areas, demonstrating capabilities of the framework and the special optimization procedures applied.

Focus Area: Multi-physics and Multi-objective Shape Optimization
This application example is a multi-physics problem from BMW Motorsport: Aerodynamic and structural shape design of the BMW M8 GTE wheel (Fig. 5). As can be seen in Fig. 6a, the spinning structure is put in a virtual wind tunnel to evaluate the aerodynamic performance. A moving reference frame (MRF) is used to account for the rotation of the wheel. A finite-volume discretization with 8 million cells forms the fluid domain and is modeled using in-compressible Navier-Stokes equations. The structural model is made up of approximately 1 million tetrahedral finite elements and rigid body elements. The fluid and structural models including boundary conditions are shown in Fig. 6. OpenFOAM and Altair OptiStruct are used to solve the primal quantities fluid and structural models, respectively. The outer surface of the spokes is chosen as the design surface, which is also the interface between the fluid and the structure models. The surface shape sensitivities of the fluid are calculated using the adjoint solver of OpenFOAM: adjointOptimizationFoam, and the structural surface shape sensitivities are calculated using adjoint technology of OptiStruct. The surface mesh from fluid and structure do not match due to the different mesh requirements for the flow and structure. The structure interface mesh consists of 28,000 nodes, while the fluid inter-   6 Computational models of the rim face mesh consists of 145,000 nodes. To deal with the non-matching meshes at the fluid-structure interface, we take advantage of the important property of the vertex morphing technique presented in Sect. 2.1 of being able to discretize the geometry and control spaces with different mesh resolutions. For this purpose, the interface mesh with less resolution, here the structure mesh with 84,000 design variables, should be chosen as the control mesh. The motivation of this choice comes from the fact that the representation of design modifications on the coarser mesh on a finer mesh is possible but not vice versa. Then, the objective functions' gradients w.r.t. the design variables (∇ s J F ,∇ s J S ) and the geometries' updates ( X F , X S ) are calculated, respectively, as where J F and J S represent objective functions whose spatial gradients on the fluid surface mesh and the structure surface mesh, respectively. The matrices A F S and A SS are the operators which define the association between the non-matching interface meshes and the discretized design space. Based on the derivation of vertex morphing in Sect. 2.1, the entry in row i and column j of A F S is computed as where S is the portion of S which lies within a sphere of radius r and center X F ,i . Multiple response functions quantify the performance of the described problem. There are multiple ways to deal with multiple objectives in an optimization problem. Some of them include weighted-sum method, weighted min-max method, lexicographic and Pareto optimization methods. In the framework presented, owing to its highly effective and simple formulation, a weighted-sum method is used. In this method, a compromise function is formed as a weighted sum of existing objectives, the total drag force (J F ) acting on the wheel, and the total strain energy (J S ) of the wheel. This forms a new objective, thus transforming a multi-objective into a singleobjective optimization problem. To give both the objectives equal importance in the optimization, equal weights are chosen for them in forming the combined objective. The following geometric constraints act on the shape optimization: (a) an equality constraint on the inner volume of the wheel (mass) and (b) an n-fold cyclic symmetry condition on the design surface. While the former is a single scalar-valued constraint,  Figure 7 explicitly focuses on the cyclic symmetry property of the wheel. As can be seen, the wheel surface is divided into five identical surfaces i , i = 1, . . . , 5, each of which is generated by rotating 1 by θ i = (i−1)π 5 around the Y-axis. In the current implementation, making the shape gradients rotationally symmetric enforces the cyclic symmetry constraint. The following linear transformation makes the shape gradients rotationally symmetric: where AP i,1 , i = 1, . . . , 5 is the transformation matrix from surface i to 1 . T i is the rotation matrix for each section, and it may be calculated as Then the shape variations in the surface are associated with the shape variations of the surface 1 as follows: dX dX 1 = AP (16)   Finally, the rotationally symmetric form of objectives' gradients can be calculated by using the following equations: The above-described process of making the surface gradients symmetric is implemented in the current frame work taking advantage of it modular nature of implementing at C++ level shown in Fig. 3. Table 1 and Fig. 8 summarizes the outcome of the described optimization problem. It shows a clear improvement in drag and compliance. One can also notice the cyclic symmetry in the resultant design geometry presented in Fig. 9. This example, in addition to demonstrating the ability to work with different solvers for physics and using multiple objective functions, also demonstrates the possibility to introduce geometrical constraints into the parametrization.

Focus Area: Flexibility and Adaptability for Manufacturing Constraints
Real-world industrial applications rarely produce unconstrained shape optimization problems. Physical and geometric constraints arise naturally due to the limitations of the physical properties or those of the manufacturing process of the mechanical component. In such applications, individual manufacturing requirements have to be fulfilled by the optimized geometry. Those are relevant to the respective manufacturing process, such as casting, molding and stamping. In the newly emerging fields like 3D printing, new types of constraints arise on the shape of the 3D printed parts.
These different types of constraints and their inclusion in the optimization process requires a flexible framework. In the following, we present two examples of such specialized constraints implemented in the current framework to demonstrate its the flexibility and adaptability.

Stamping
This example presents the inclusion of the stamping constraint. The optimization problem consists of a load-bearing member of a car subjected to respective load conditions, as shown in Fig. 10c. The member is shape-optimized for minimal strain energy. In addition to a mass constraint, the stamping constraint is applied to the optimization problem. Considering the surface of the component as design geometry, Altair OptiStruct is used for structural calculation of primal and shape sensitivities. A mesh of 10500 nodes is used for the computation.
There are multiple ways to impose the stamping manufacturability constraint. This framework approaches the problem from a geometric point of view. The following condition describes a perfectly stampable geometry: for every node, when a ray trace is shot from it in the stamping direction, the ray should not collide with any element of the mesh. In other words, looking at the mesh in the stamping direction, all nodes should be visible. The volume V h formed by the hidden parts of the mesh forms the constraint objective. Such a constraint is called a soft constraint as it is not satisfied node wise but on average. Figure 10a shows the hidden volume V h visually. Numerically, V h is calculated using ray-tracing. A novel radius search in the mesh is employed to reduce the computational costs induced by ray tracing. Figure 10b shows the optimized cross section of the component with and without stamping constraints. The described geometric constraint is specifically designed to ensure that the optimized geometry can be produced using a stamping manufacturing process. Here, the geometric formulation replaces a much more elaborate simulation of the manufacturing process. This application case shows the ease of integration of such specific constraint formulations into the optimization process.

Stackability
The emerging and established techniques in the field of additive manufacturing and 3D printing impose a vital requirement of stackability on the components being printed. This requirement is more pronounced when using powder bed-based techniques in metal printing. A stackable component allows printing of more components per unit area. This transforms into a geometrical requirement of placing two consecutive parts as close as possible to each other. In the current framework, this is realized as a geometric no penetration constraint. The current framework operating at BMW Group employs this technique successfully to different 3D-printed geometries.
An exemplary application is presented here. It consists of a structural component subjected to pseudo-load conditions. As shown in Fig. 11a, the component is to be Fig. 11 Stackability test case stackable in the given direction. Altair OptiStruct is used for simulation of the structural model with 22,000 nodes. The shape of the component is optimized for minimal strain energy subjected to constraints on stackability. The stackability constraint is implemented as a solver using the solver interface "SolverWrapper" of the framework. This solver instead of solving for the primal values will only perform the necessary geometric operations required for calculating the sensitivities for the constraint. Figure 11b shows the optimized component after 20 optimization iterations. As seen here, the shape of the component respects the no penetration with the two adjacent components and thus is more suitable for mass production using 3D printing technologies.
This application uniquely showcases the re-usability of existing features, in this case, geometric no penetration constraints, to formulate and realize new optimization problems to include diverse application fields.

Focus Area: Communication and Usage of Multiple Solvers
Multi-physics simulations have been an essential tool in understanding complex physical phenomena. In this context, multi-physics optimization is also equally important. An effortless way for performing such optimization is to use the existing primal solvers for calculations. This approach not only is modular but reuses the expertise of physical solvers independent of each other without introducing more complexity.
These multi-physics optimization problems introduce a new level of complexity to the enabling framework as these problems need to work with multiple solvers simultaneously and choreograph them in a optimization algorithm. The following two application cases make evident the capability of the software design of the optimization framework to work with such complex problems.

FSI Shape Optimization: Flexible ONERA M6 Wing
The first example in this focus area is a multi-objective and multi-disciplinary shape optimization of a flexible ONERA M6 wing immersed in a compressible inviscid fluid flow. Both the fluid analysis (CFD) and the structural analysis (CSD) assume steady cruise conditions. For the CFD analysis, the compressible solver from SU2 [30] is used and the CSD analysis is carried out in KratosMultiphysics framework [25]. The wing structure is clamped at the wing root, as shown in Fig. 12. The steadystate transonic flow over the ONERA M6 wing at Mach 0.8395 and angle of attack of 3.06 • is computed using nonlinear Euler equations. The wing structure is purposefully modeled to produce large deformations and uses 4-node tetrahedral nonlinear solid elements. The flexible structure for the wing introduces fluid-structure interaction into the model. So the corresponding shape sensitivity analysis becomes an aeroelastic problem and requires a coupled sensitivity analysis to be performed. The optimization problem is formulated with the lift to drag ratio (efficiency) from fluid and strain energy as a weighted objective and the inner volume of the wing constraints this optimization problem. The shape sensitivities for these objectives are also calculated in the respective solvers. Naturally, the fluid and the structural domains have two different levels of refinement. As in the previous application of Sect. 4.1, the coarser wing surface mesh (structure surface mesh) is used to parametrize the surface using vertex morphing technique. Therefore, the filtering operations presented in Eqs. 11 and 12 are used in the optimization process.
The optimization resulted in a 32.4% increase in the lift to drag ratio and a 52% decrease in the total structural strain energy. The optimization history is presented in Fig. 13. Furthermore, Fig. 14a compares the optimized (scaled deformations) and initial configuration of the wing. As seen in Fig. 14b, the strong shock wave that existed along the span is reduced significantly. A detailed discussion about the methodology followed and results is presented in a parallel publication by one of the authors [3].
This example shows the applicability of the framework even for highly complex shape optimization with a fully coupled multi-physics problems involved. The shape update of the non-matching computational meshes is controlled using the vertex morphing method without additional modeling effort.

Multi-disciplinary Shape Optimization: Spiral Water Jacket
Further expanding upon the use of vertex morphing in multi-disciplinary applications, the optimization of a spiral water jacket is investigated. Optimal design of an electric motor water jacket, seen in Fig. 15b, is a de facto multi-disciplinary problem requiring the consideration of many different physics. The primary role of the water jacket is to cool the motor. This cooling effect in turn improves power output and the life of components in and around the motor. The initial geometry, as shown in Fig. 15b, represents an internal flow region within the motor housing. This application of shape optimization with vertex morphing considers a simplified optimization considering only the flow and static structural properties. The objective of the optimization is to minimize the pressure drop while respecting a maximum stress constraint. Besides, an internal geometric restriction was placed on the geometry to maintain a minimum internal wall thickness between the water jacket and the stator located inside the house as previous runs had shown this to be a problem.
The flow and thermal problem are solved using a monolithic conjugate heat transfer (CHT) problem using the commercial solver Star-CCM+, while the structural analysis is carried out in the commercial solver Optistruct. The consideration of the thermal load in the structural solver was the source of the coupling between the models. The temperature field calculated in Star-CCM+ is applied to the Optistruct model via a mapping operation also implemented in optimization framework described. The design surface is the surface of the water jacket and the 22400 nodal coordinates on this surface act as control variables.  The optimization was carried out for 10 steps before it stops because of meshrelated issues. Though geometric changes are small, the improvement in the pressure drop amounted to a reduction of about 30%. The main changes to the geometry are around the sand core support geometries. This is shown in Fig. 15a. These features exist due to the sand core casting technique used to manufacture the part. Although there are improvements in this case, a proper adjoint coupled problem may be more successful at maintaining stresses below the permissible limit. Additionally, further constraints and objectives representing the manufacturing technique may be required as the changes seen here are likely to have a direct impact on the manufacturing process as well. In this application case, in addition to the communication with two external solvers, the transfer of the temperature field from fluid to structural models is also performed using the current framework. This application is a prime example of the flexibility and generality of the interfaces and the implementation.

Computational Performance
Quantification of the computational performance of an optimization process is not straightforward as the components that consume the lion's share of the computation time are the individual physics solvers computing the primal, possible response gradients. A discussion on performance of such solvers is a discipline on its own and is out of scope of current contribution as they are treated as black-boxes. Therefore, Fig. 17 Scaling for vertex morphing in the optimization framework this subsection focuses on the performance of the framework itself; this includes the operations of 8 till 15 of Algorithm 1. To this extent, an exemplary test case of a flat square surface with pseudo-sensitivity information is taken to demonstrate the computational capabilities of the optimization framework. An additional factor that strongly influences the performance of the implemented vertex morphing technique is the filter radius illustrated in Fig. 2. This determines the average number of neighbors, for each node, which are inside the radius. Figure 17a shows the strong scaling plot of the framework while keeping the average number of neighbors constant. Figure 17b also shows the performance of the framework for the case with 1,000,000 nodes and an increasing number of neighbors. As expected, the computational cost increases as the number of neighbors increases, but only linearly. This shows that the framework is efficient in both regards.

Conclusion and Outlook
In this work, we introduced a concept of a framework for performing shape optimization using vertex morphing. The framework offers a flexible interface to use different external physics solvers. For this, an innovative detached-interface approach is introduced. Multiple commercial and open-source tools used in the applications give an insight into the generality of this approach. It also has the possibility to include geometric and other special constraints in the optimization process. Wide range of application cases with different levels of the complexity, especially involving multi-physics simulations together with geometrical constraints, prove the flexibility and robustness of the framework presented. This design of the framework drawing on the object-oriented features of C++, and the scripting nature of Python programming languages provide a high degree of flexibility for the developers and practitioners to realize new optimization problems. The efficiency and the capability of the framework in handling a large number of design variables are also showcased with the applications.
As an outlook, coupled adjoint sensitivities are planned be used in the example presented in Sect. 4.3.2 to achieve more improvement in the objective. A constant effort in including more commercial solver is done to increase the number of applications. Further additions to the specific industrial production constraints are in progress together with industry partners.