{'fraction_non_alphanumeric': 0.026291134793301325, 'fraction_numerical': 0.013096771675059531, 'mean_word_length': 4.641378550979961, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topologically robust CAD model generation for structural optimisation

Computer-aided design (CAD) models play a crucial role in the design, manufacturing and maintenance of products. Therefore, the mesh-based finite element descriptions common in structural optimisation must be first translated into CAD models. Currently, this can at best be performed semi-manually. We propose a fully automated and topologically accurate approach to synthesise a structurally-sound parametric CAD model from topology optimised finite element models. Our solution is to first convert the topology optimised structure into a spatial frame structure and then to regenerate it in a CAD system using standard constructive solid geometry (CSG) operations. The obtained parametric CAD models are compact, that is, have as few as possible geometric parameters, which makes them ideal for editing and further processing within a CAD system. The critical task of converting the topology optimised structure into an optimal spatial frame structure is accomplished in several steps. We first generate from the topology optimised voxel model a one-voxel-wide voxel chain model using a topology-preserving skeletonisation algorithm from digital topology. The undirected graph defined by the voxel chain model yields a spatial frame structure after processing it with standard graph algorithms. Subsequently, we optimise the cross-sections and layout of the frame members to recover its optimality, which may have been compromised during the conversion process. At last, we generate the obtained frame structure in a CAD system by repeatedly combining primitive solids, like cylinders and spheres, using boolean operations. The resulting solid model is a boundary representation (B-Rep) consisting of trimmed non-uniform rational B-spline (NURBS) curves and surfaces.


Introduction
The application of structural optimisation in industrial design requires the optimised geometries to be converted into CAD models. This is necessary because CAD systems are today an integral part of most product development processes. The prevalent parametric CAD systems are based on boundary representation (B-rep) techniques and trimmed NURBS (non-uniform rational B-splines). Geometries are created using constructive solid geometry (CSG) by recursively combining primitive shapes, such as cylinders, cubes etc. using boolean operations. This recursive process is stored in the form of a CSG tree. Although the input to most structural optimisation is a CAD geometry, its output is usually a finite element mesh of the optimised geometry. The finite element meshes consist of too many elements to be processed and edited in a CAD system. For a geometry to be editable in a CAD system, a compact representation in the form of a CSG tree is needed. The need to edit the optimised part geometry arises when additional geometric features have to be added or when the part is to be combined with other components into a functional product. Moreover, in industrial practice there are usually, in addition to structural efficiency and robustness, many equally important, often not explicitly quantified, design requirements, which require the geometry to be edited after optimisation.
The fully automated process we propose can robustly generate a compact CSG tree representation of topology optimised geometries and convert them to CAD models. The entire workflow is illustrated with the help of the  The rocker arm is a lever mechanism which pivots around the two of the joints shown in (a) when a force is applied at the other four joints. In (e) and (f) the six detailed joint designs, with each containing a roller bearing, have been added in a CAD system after model generation.
topology optimisation of the rocker arm shown in Figure 1. For topology optimisation we use a standard density based SIMP (solid isotropic material with penalisation) approach on a structured hexahedral finite element mesh to be referred to as a voxel mesh, see e.g. [1][2][3]. The optimised density field splits the set of voxels after thresholding into the two subsets solid and void, which gives a three-dimensional binary picture with two grey levels. Subsequently, we extract from the binary image a voxel chain skeleton using a homotopic, i.e. topology-preserving, skeletonisation algorithm from digital topology [4]. The skeleton has the same number of connected components, handles and cavities as the three-dimensional voxel model. Evidently, topology preservation is critical in structural applications because a change of topology may disrupt the load paths discovered during optimisation. The voxel chain model defines a weighted undirected graph which is processed with graph algorithms to obtain a spatial frame structure. Each of the members of the frame structure is a beam and the beams are rigidly connected at the joints. During the conversion of the voxel model to the spatial frame structure the optimality of the structure is usually compromised. This can, however, be recovered after applying a few steps of size and layout optimisation to the spatial frame structure. Size optimisation adjusts the cross-sections of the members and layout optimisation adjusts the coordinates of the joints. The very compact representation which we generate from the original voxel model consists of the connectivity of the frame structure, its joint coordinates and the member cross-sections. These are used to create a binary CSG tree and a solid model of the spatial frame structure in a CAD system. Figures 1e and 1f show two different solid models generated in a CAD system. In the first model the cylindrical members are connected to spherical joints. In the second model the members are still cylinders, which are, however, smoothly blended at the joints. With the process presented in this paper, the first model is generated in a fully automated fashion while the second requires some user intervention. Note although we used a density based approach for topology optimisation, it is straightforward to apply the proposed approach to the more recent level-set based methods [5,6] or the historical homogenisation method [7], see also [8] for an overview on optimisation approaches. There are only a very limited number of structural optimisation approaches which arrive at a geometry in the form of a compact parametric CAD model. These approaches are primarily restricted to shape optimisation. For instance, it is quite common to use the parameters of a CAD or a reconstructed CAD-like spline model as geometric design variables in shape optimisation, see e.g. [9][10][11][12]. Such techniques are especially appealing in an isogeometric analysis context when the same basis functions are used for geometry description and finite element analysis [13][14][15][16][17][18]. The compact spline representation of the optimised geometry can in principle be imported into a CAD system. In the past, there has been some work in topology optimisation on generating a spline representation based on a given optimised finite element geometry. A common approach is to fit the density isocontours obtained with topology optimisation with a spline curve in two dimensions (2d) or spline surface in three dimensions (3d) [19][20][21]. Especially in 3d, these techniques are not robust as they may require the solution of a nonlinear least squares problem and are hard to automate due to the need to manually position the control vertices of the splines. Possibly, therefore, they had no noteworthy impact on commercial software despite being developed around two decades ago. The use of skeletonisation, or thinning, for compact geometric representations has earlier been pioneered in Bremicker et al. [22] for 2d and very recently in Cuillière [23] for 3d. While in [22] the aim was not to arrive at a parametric CAD model, skeletonisation algorithms from digital topology were used to extract a truss structure, which was subsequently size and layout optimised. Homotopic skeletonisation in 3d is substantially more challenging than in 2d, and elimination of possible mechanisms in a spatial truss structure is far from trivial [24]. And so, in our approach, we convert the voxel chain model into a spatial frame structure (with rigidly connected joints) which intrinsically does not exhibit mechanisms.
In [23] a curve skeleton is extracted from the surface mesh representing the density isocontour of the optimised geometry using the skeleton extraction method proposed by Au et al. [25]. This specific extraction technique is rather elaborate as it builds on the successive contraction of a given surface mesh using Laplacian smoothing. The obtained skeleton and estimated cross-sections are used to generate a CAD model with no further postprocessing. It is not taken into account that skeletonisation may have impaired the optimality of the structure. A manufacturing oriented overview of compact parametric representations for topology optimised geometries can be found in the recent review [26]. An essential component of the proposed CAD model generation process is skeletonisation, which is an active research area with applications in computer graphics, animation and volumetric image processing, see the reviews [27][28][29][30]. The curve-skeleton, briefly the skeleton, of a 3d object is closely related to its medial surface, or surface-skeleton. In finite elements medial surfaces are known from mesh generation applications, see e.g. [31]. Informally, the medial surface is the set of all points that have two or more closest points on a 3d object's surface. That is, a sphere centred on the medial surface touches the object's surface at two or more points. In contrast to the medial surface, the skeleton of a 3d object has no rigorous definition. It is expected that the skeleton captures the essential topology of the object and is centred, i.e. lies on or close to the medial surface. There are many algorithms for determining the skeleton of an object starting from different types of geometry representations, such as implicit signed distances, polygonal surfaces etc. In topology optimisation, it is expedient to assume that the object is given as a voxelised binary image consisting of solid and void voxels. Digital topology provides a rigorous basis to study the topology of binary images [32]. Homotopic skeletonisation algorithms based on digital topology can extract a one-voxel-wide voxel chain skeleton with the same topology as the object, i.e. the same number of connected components, tunnels and cavities [27,33]. They usually start from the domain boundaries and proceed by iteratively removing voxels one by one that do not alter the topology of the object. When the algorithm terminates, the entire structure is represented by a network of one layer thick voxel chains. It is sufficient to inspect only a voxel's immediate neighbourhood consisting of 3×3×3 voxels to decide whether a voxel can be removed [4]. The possible states of such a small cluster can effectively be encoded in a look-up table. Crucially, the entire skeletonisation process does not rely on any floating point operations, which makes it exceedingly robust.
The outline of this paper is as follows. In Section 2, the standard density based topology optimisation and the size and layout optimisation of spatial frames are briefly reviewed. Subsequently, in Section 3 relevant aspects of digital topology are introduced and the implemented specific skeletonisation algorithm is described. The robustness and runtime of our implementation are studied with a relatively complex quadcopter frame geometry. The postprocessing of the obtained voxel chain model with graph algorithms and the generation of first a frame structure and then a CAD model are discussed in Section 4. The entire optimisation and geometry conversion process is summarised in Section 5. The application of the proposed approach to a standard 2d benchmark example from topology optimisation and more complex 3d examples is illustrated in Section 6. We study in particular the change of the compliance cost function during the conversion of the optimised voxel model to a CAD model.

Review of structural optimisation
In this Section, we briefly review the standard density based topology optimisation using the SIMP method and the size and layout optimisation of spatial frame structures. In all cases, the cost function is the compliance, and the discussion is restricted to aspects relevant to this paper. For further details see, e.g., the monographs [3,34].

Topology optimisation of solids
The topology optimisation problem for finite element discretised solids reads where J(ρ) is the compliance cost function, ρ is the vector of relative element densities, u is the displacement vector, K is the global stiffness matrix, f is the global external force vector, V(ρ) is the material volume, V is the design domain volume and V f is the scalar prescribed volume fraction. The relative density of each element is constrained to be 0 ≤ ρ i ≤ 1 with i ∈ {1, . . . , n e } and the number of elements is n e . As usual, the global stiffness matrix K and vector f are assembled from the n e element contributions K i and f i in the mesh. In this paper, we discretise the design domain always with a structured grid and hexahedral linear elements. In each element i the material is isotropic and homogeneous, and the Young's modulus E is penalised depending on the relative density ρ i with where E is the prescribed Young's modulus of the solid material, p and E min are two algorithmic parameters. The penalisation parameter p ≥ 3 ensures that elements with densities close to ρ i = 0 (void) and ρ i = 1 (solid) are preferred. The small Young's modulus E min ≈ 10 −9 of the void material prevents ill-conditioning of the global stiffness matrix when ρ i = 0. Each element stiffness matrix K i is computed using the corresponding Young's modulus E(ρ i ) with the relative density ρ i , which is constant within an element. Furthermore, the topology optimisation problem (1) is regularised by filtering the element densities ρ i to prevent checker-board instabilities and mesh dependency of the solution. This is accomplished by convolving the element densities with the kernel function which depends on the coordinates of the centroids x i and x j of the elements i and j, and the prescribed filter length R.
With H(x i , x j ) at hand the filtered densities are given bŷ assuming that all elements have the same volume. For gradient-based optimisation the derivatives of the cost and constraint functions in (1) with respect to the relative densities ρ i are needed. In the following the relative densities ρ i in the topology optimisation problem (1) are replaced with the filtered relative densitiesρ i . The derivative, or sensitivity, of the compliance cost function (1a) reads Differentiating and rearranging the equilibrium constraint (1b) gives and the derivative of the filtered densities (4) is Introducing both derivatives into (5) yields where the stiffness matrix derivative is to be assembled from element contributions considering the penalised Young's modulus (2). The derivative of the volume constraint (1c) is obtained analogously with It is assumed here that all the elements are of the same size, namely the design volume V divided by the total number of elements n e . Finally, with the obtained derivative of the compliance cost function (8) and the volume constraint (9) the optimised density distribution is determined iteratively with the optimality criteria method, see e.g. [3].

Size and layout optimisation of frames
We consider spatial frame structures consisting of straight beam members connected by joints that can transfer forces and moments. The members can deform by stretching, bending and torsion and are modelled as classical Timoshenko beams, see Appendix A. Without loss of generality, size and layout optimisation can be posed as an iterative sequential optimisation problem. In each optimisation step, either the size, i.e. cross-section area, or layout, i.e. joint coordinates, of the members is optimised.
Both the size and layout optimisation problems have the same structure as topology optimisation (1), namely, The vector of design variables s refers either to cross-section areas in case of size optimisation or joint coordinates in case of layout optimisation. The two bounds s l and s u have different interpretations in the two cases. The frame consists of n e beam elements with the element stiffness matrices K i and vectors f i . Each stiffness matrix K i is obtained from a local stiffness matrix K l i formulated in a coordinate system attached to the element with the index i. The local matrix K l i ∈ R 12×12 with the displacements and rotations of the two end nodes of the beam element as degrees of freedom is given in Appendix A. In the local x l i , y l i and z l i coordinate system the beam axis is assumed to be aligned with the x l i axis. The matrix K l i is transformed to the global x, y and z coordinate system according to with the block-diagonal transformation matrix The rotation matrix λ i ∈ R 3×3 can be chosen with which is composed of the two elemental rotations α i around the z l i axis and β i around y l i axis. Similar as in topology optimisation, the derivative, or sensitivity, of the compliance cost function (10a) reads where u j is the vector of nodal displacements and rotations of the element j. In size optimisation the derivative of the element stiffness matrix (11) is given by whereas in layout optimisation it is

Digital topology and skeletonisation
The input to skeletonisation is a structured hexahedral finite element grid with all the vertices within the domain having eight incident hexahedra. For digital topology only the connectivity of the mesh is relevant whereas the coordinates of the grid nodes can be arbitrary so that in the following hexahedron and voxel are used interchangeably. The hexahedral grid encodes the binary voxel model of the topology optimised structure. The binary voxel model is obtained by denoting voxel with a density ρ i above the threshold η as solid and the remaining ones as empty (void). Hence, the set of all voxels v ∈ V = V s ∪ V e is composed of the solid voxels V s and empty voxels V e . Skeletonisation aims to reassign voxels one by one from V s to V e until the solid is represented by a network of one layer thick voxel chains. The decision on which voxels to reassign is guided by digital topology and proceeds starting from the voxels at the boundary between V s and V e .

Review of digital topology
To retain the topology of the voxel model, its Euler characteristic is considered. The solid voxels in V s and empty (void) voxels in V e yield two hexahedral volume meshes, or more abstractly cell complexes, which will be in the following denoted with M s = M(V s ) and M e = M(V e ). The volume mesh of a voxel set consists of the corresponding set of the vertices, edges, faces and voxels. It is also useful to define quadrilateral surface meshes ∂M s and ∂M e formed by the boundaries of M s and M e , respectively. The Euler characteristic of the voxel model can be determined either from the volume mesh M s with or from the surface mesh ∂M s with where n 0 denotes the number of vertices, n 1 the number of edges, n 2 the number of the faces and n 3 the number of hexahedrons in M s or ∂M s . In Figure 2 the Euler characteristics of a few voxel meshes are given. Notice that there is also the relation see, e.g., textbooks [35,36] on topology. For instance, the Euler characteristic χ(M s ) of a sphere is 1, of a torus is 0, of a double torus is −1 and so on, compare also Figure 2. Hence, the total number of different entities in the mesh and the described shape are intrinsically connected.
In determining the Euler characteristic of a voxel model, two different types of neighbourship definitions are possible, see Figure 3. These definitions define how pairs of touching voxels such as shown in Figure 3c and 3d are treated, whether connected and forming one surface or disconnected and forming two surfaces. That is, the number of entities in (17) and (18) depend on the chosen definition of the neighbourship. For a voxel v its 6-neighbours consist of all the cells which share a face with v and its 26-neighbours consist of all the cells which share a common vertex with v. The voxels in the 6-neighbourhood of v are denoted with N 6 (v) and the ones in the 26-neighbourhood with N 26 (v). In our implementation, two solid voxels are considered to be neighbours when they are 26-neighbours. In contrast, two void voxels are considered to be neighbours when they are 6-neighbours. It is essential not to use the same neighbourship definitions for both the solid and void voxels. The use of one single neighbourhood definition leads to ambiguities, for instance, in the planar case to the violation of the Jordan curve theorem [32]. The Jordan curve theorem states that a simple closed curve divides the plane into an interior and an exterior region.
The Euler characteristic (17) and (18) for the entire voxel model can be determined by looping over all vertices in the grid and only considering in each step the voxels attached to one vertex. The 2×2×2 voxels attached to a vertex c are defined as its octant O(c). 1 There are as many octants as vertices in the volume grid (assuming that the voxel domain is padded with empty ghost voxels). The octants are overlapping and the edges, faces and voxels in the to the Euler characteristic. Example of an octant with two voxels.
hexahedral mesh M s appear in several octants. Hence, in computing the Euler characteristic (17) by iterating over the vertices the contribution of each octant has to be suitably weighted such that where n (c) 1 is the number of edges and n

Skeletonisation by successive voxel removal
The Euler characteristic χ(M s ) is critical in determining the voxels at the boundary of the solid mesh M s that can be deleted without changing the topology of the voxel model. A solid border voxel is defined as a solid voxel with at least one void voxel amongst its 6-neighbours. There is a large amount of research on voxels, also referred to as simple points, which can be removed without changing topology [32]. As evident from (20) where ∆ denotes the change of the respective quantities with and without voxel v presents. It is straightforward to determine the Euler characteristic (22a) and it is relatively easy to determine the number of objects (22c) in N 26 (v). The change of the Euler number is determined by summing the tabulated contributions of the eight octants belonging to the eight corners of the voxel v from [4], see (21) and Figure 5. To determine number of objects (22c), we generate an undirected graph with the solid voxels as nodes and introducing edges between voxels that are 26-neighbours. Subsequently, the number of connected components is determined with a depth first search (DFS) algorithm [38]. The skeletonisation proceeds by removing in each step one layer of border voxels which are simple points according to (22). The removed solid voxels are reclassified as void. One removal step is split into six sub-steps and in each sub-step only the voxels approaching from one of the six grid directions are removed. The geometry of the obtained skeleton depends somewhat on the sequencing of these sub-steps, which is not critical for our purpose. It is straightforward to tag some voxels as non-removable irrespective of (22). In skeletonising topology optimisation results, for instance, the voxels at Dirichlet and non-homogeneous Neumann boundaries are tagged as non-removable. Also, end points with only one voxel in their 26-neighbourhood have to be tagged as non-removable, in order to avoid , see [4]. Hence, the voxel cannot be removed because this would lead to a change in the Euler characteristic by ∆χ = 1. a complete elimination of lines. The skeletonisation terminates when none of the remaining voxels is removable without violating (22). The final result is a curve skeleton consisting of a network of one layer thick voxels connected by joints. We refer to this curve skeleton as the voxel chain skeleton.

Illustrative example and timing
We consider the skeletonisation of the quadcopter frame shown in Figure 6a obtained from GrabCAD. The frame has a non-trivial topology and has been designed in a CAD system. To generate a voxel model the implicit, or level set, representation of the frame on a uniform voxel grid is computed first. This is achieved by determining the distance of each voxel in the grid to the frame surface. In doing so, the frame geometry is approximated with the STL mesh exported from the CAD system depicted in Figure 6b. After thresholding the voxel grid the voxel model of the frame in Figure 6c is obtained. Its skeletonisation with the introduced digital topology algorithm yields the voxel chain skeleton in Figure 6d. It can be visually confirmed that the skeleton appears to have the same topology as the voxel model.
To investigate the efficiency and scaling of the introduced skeletonisation algorithm, we consider three different voxel grids for computing the implicit representation and skeletonisation. The size of the three grids, the number of voxels in the different models and the runtimes are given in Table 1. The STL mesh has in all cases 1086791 triangles. The number of solid voxels in the voxel model and the number of removal steps increases with decreasing voxel size because more and more cover the frame. The number of skeleton voxels increases because the smaller voxels can capture the topology of the frame better. The reported extremely short runtimes include only skeletonisation (no implicitisation) and confirm the efficiency of the skeletonisation algorithm. Moreover, notice that the average time per removal step is approximately linear with respect to the number of the voxels in the voxel grid.

Model extraction and conversion
The voxel chain skeleton is used to define a structural frame model. The skeleton provides both the connectivity and the geometry of the frame. Although it is feasible to obtain also the member cross-sections from the voxel model, in our implementation cross sections are obtained from size optimisation of the frame, as will be discussed in Section 5. In converting the voxel chain skeleton to a frame we express it as a weighted undirected graph and    Table 1. make use of its incidence matrix. The frame model provides a very compact representation of the optimised structure, which can be used to generate a parametric solid model using the scripting interface or API (application programming interface) of a CAD system.

Graph model
The voxel chain obtained from skeletonisation contains two different types of voxels, see Figure 7. We refer to the voxels with only two voxels in their 26-neighbourhood as regular voxels. The remaining voxels are the joint voxels, which have other than two voxels in their 26-neighbourhood. As illustrated in Figure 7b, it is not expedient to deduce from every joint voxel a joint for the frame model because this will lead to acutely short members and impractical joint designs.
To robustly merge joint voxels in close proximity to a single joint we resort to a weighted undirected graph representation of the voxel chain skeleton. The nodes of the graph are the joint voxels, and the edges are the connections between the joint voxels. Each of the nodes has an associated coordinate vector which is initially equal to the centroid of the corresponding joint voxel. The edge weights are proportional to the number of voxels between the two joints attached to an edge. The graph model is generated in two steps. First, all the solid voxels in the voxel model are categorised either as regular or joint. Subsequently, starting from the joint voxels we determine the voxel chains and their lengths by marching along the N 26 (v) neighbours until a joint voxel is reached. The obtained graph model is stored in the form of an incidence matrix. In Figure 8a and 8d the graph model and its incidence matrix for the voxel chain skeleton in Figure 7b are given. The rows of the incidence matrix correspond to the nodes of the graph and the columns to its edges. Each edge has two identical non-zero entries corresponding to the length of the edge in voxel units.  The joint voxels in close proximity of each other are connected by very short edges. The corresponding graph nodes are merged by successively collapsing the edges of the graph. The edge collapse is implemented with the help of the incidence matrix. Figure 8 illustrates how the incidence matrix evolves during the collapse of the edge e 3 and the merging of its attached nodes v 3 and v 4 . The coordinate vector associated to the node v 3 is the average of the node coordinate vectors before merging. After collapse, the duplicate edge e 5 is detected and is merged with e 4 . The weight of e 4 is updated to be the closest integer to the mean of the weights of e 3 and e 4 . The voxel chain model usually contains branches which lead to structural frame members with zero stress. These are the chains directly connected to a homogeneous Neumann boundary. Their pruning is again implemented with the incidence matrix. As illustrated in Figure 9, after identifying from the incidence matrix the nodes with only one attached edge, nodes and edges are deleted by simply removing the corresponding rows and columns.

Structural frame and CAD models
The graph model supplemented by the joint coordinates and member cross-sections provides sufficient information to generate a compact structural frame model. We assume that there is a one-to-one correspondence between nodes and edges of the graph model and joints and members of the frame model. Besides, it is assumed that members are only subjected to end forces and moments, i.e. no distributed loads, so that they can be straight and have uniform cross-sections along their lengths. Note that these assumptions are fully justified if there is no distributed loading in topology optimisation. Furthermore, for simplicity, we assume in the following that all cross-sections are circular. These restrictions can be mostly relaxed if necessary.
It is straightforward to create a parametric solid CAD model from the structural frame model. To this end, we use either the commercial CAD system Rhinoceros (Rhino) or the open source FreeCAD. However, the proposed approach can be realised in any CAD system which provides a scripting interface or API (application programming interface). In the solid model, the members are represented by cylinders and the joints by spheres, see Figure 10. They are combined with boolean operations. The radius of the sphere at a joint is chosen slightly larger than the largest attached member radius. In this paper, we use a factor of 1.05. A binary CSG tree of the frame model is generated by starting with one of the members and by successively adding spheres and cylinders to the model. Different variations of this approach can be devised and implemented. For instance, it suggests itself to connect two of the most stressed members continuously across a joint or to add fillets to reduce stress concentrations. See also [39,40] in which a somewhat similar approach was used for generating CAD models for optimised truss structures. After the model is generated it can be edited and postprocessed in the CAD system and be exported, for instance, in STL format for machining or IGES or STEP formats for mesh generation software.

Overall optimisation and conversion process
In this Section, we review the sequence of steps from the definition of a topology optimisation problem to obtaining a structurally-sound parametric CAD geometry. We provide additional details for each step and focus on the interplay  of the different steps. Each step corresponds to one Section in this paper.
Step 1: Topology optimisation. We assume that the finite element mesh for the optimisation problem (1) is a structured hexahedral grid. In most optimisation problems, the design domain is a parallelepiped which can be discretised with a structured grid. Other more complex design domains can be considered by computing their implicit, or level set, representation and embedding them in a structured hexahedral grid, see the quadcopter frame example in Figure 6. From the outset, the voxels outside the design domain are chosen as void by choosing their Young's modulus with E min .
Step 2: Skeletonisation. The input to the homotopic skeletonisation algorithm is a binary image defined on a hexahedral structured grid. The binary image is obtained by thresholding the topology optimised geometry. The threshold value η is chosen such that the prescribed material volume constraint (1c) in topology optimisation is satisfied. The output of the skeletonisation is a curve skeleton in the form of a voxel chain with the same topology as the topology optimised geometry.
Step 3: Structural frame model generation. The topology and joint coordinates of the frame model are extracted from the voxel chain model. It is assumed that all the members are straight, have a circular cross-section with the same diameter and their total volume is equal to the prescribed material volume V f V in topology optimisation. It is possible to obtain frame models with more complex geometries and non-uniform cross-sections. This appears to be, however, usually not desired from an ease and cost of manufacturability viewpoint.
Step 4: Size and layout optimisation. The structural frame model extracted from the skeleton is usually suboptimal. That is, the compliance of the structural frame is larger than the compliance of the topology optimised geometry. To recover the optimality of the frame structure, we apply several steps of sequential size and layout optimisation (10). In size optimisation the diameters d i of the circular member cross-sections and in layout optimisation the coordinates x i of the joints are updated. Both optimisation problems are solved with the SQP (sequential quadratic programming) method [41]. Throughout optimisation the volume of the frame is constrained to be equal to the prescribed material volume in topology optimisation V f V. It is straightforward to consider additional constraints pertaining to the member cross-sections, the positions of the joints, or positions and orientations of the members. The first step is always size optimisation, which is followed by as many as necessary alternating layout and size optimisation steps until the compliance cost function is converged. If the length of any member reduces to zero during shape optimisation, the member is removed, its end nodes are merged, and the iteration continues.
Step 5: CAD model generation. A compact CAD model of the structural frame can be generated essentially in any parametric CAD system using a fully automated process. The members are represented by cylinders and the joints by spheres which are combined by boolean operations. The underlying binary CSG tree representation makes it easy to edit further and to refine the optimised design.

Examples
We consider three examples of increasing complexity to demonstrate the application of the proposed approach. The minimised cost function is the compliance J(ρ) in topology optimisation and J(d, x) in size and layout optimisation, whereρ are the filtered element densities, d are the member diameters, and x are the joint coordinates. In all examples, the Young's modulus and Poisson's ratio of the solid material are E = 2.1 · 10 5 and ν = 0.3. In topology optimisation the penalisation factor is chosen as p = 3 and the minimum Young's modulus as E min = 10 −9 . During layout optimisation, a member shorter than 1/20 of the total length of all members sharing the same node is considered as too short, and its two end nodes are merged.

Cantilever
Cantilevered plates as shown in Figure 11 are one of the most widely studied benchmark examples in topology optimisation, see e.g. [3]. The length, height and width of the chosen design domain are 150 × 50 × 4. The left face of the domain is clamped while all the other faces are free, and a distributed force with a total value of F = 100 is applied along the centre of the right face. The finite element discretisation consists of 150 × 50 × 4 linear hexahedral elements. The maximum volume fraction of the optimised structure is prescribed to be V f = 0.3 and the density filter radius is chosen as R = 3. Figure 12a depicts the optimised cantilever with only the voxels above a relative density η = 0.3. The compliance of this voxel model is J(ρ) = 3.40900, which has been determined after topology optimisation by setting p = 1 in (2). As a first step in obtaining the structural frame model the homotopic skeletonisation algorithm is employed. The skeletonisation yields after 5×6 voxel removal steps the voxel chain skeleton shown in Figure 12b. As visually evident, the voxel and voxel chain models have the same topology. The skeleton is converted into the structure depicted in Figure 13a by first identifying the 16 joints and then connecting them with 21 members. As evident from Figure 13a during the conversion to a frame model the optimality of the voxel model is compromised; note, for instance, the nonstraight top and bottom members. Optimality is recovered by updating the joint positions, i.e. layout optimisation, and the cross-sectional areas, i.e. size optimisation, see Figure 13b. Initially, all members are assumed to have the same diameter of d = 16.24, giving the prescribed total volume of V = 0.3V. In determining the volume only the member cross-sections and member lengths between the joints are considered.
According to Figure 14b, the frame consisting of beams with uniform diameter has a compliance J(d, x) = 4.59841, which is more than the voxel model. The stretch and bending strain energies of the frame are 2.02202 and 0.247108, respectively. The minimum and maximum axial stresses are 0 and 20.9274. The average axial stress is 8.04859 and its standard deviation is 6.26102. After the first size optimisation step the compliance is reduced to J(d, x) = 3.59948 and the subsequent layout optimisation step to J(d, x) = 2.95264. Several more steps of size and layout optimisation do not lead to a significant reduction in the compliance. Notice that the obtained final compliance J(d, x) = 2.90323 is significantly lower than the compliance J(ρ) = 3.40900 of the topology optimised voxel model. The stretch and bending strain energies of the final frame are 1.42114 and 0.0190822, respectively. The minimum and maximum axial stresses are 7.00854 and 8.46234. The average axial stress is 8.09638 and its standard deviation is 0.44709. In the final optimised cantilever in Figure 13b the axial stresses become more evenly distributed and the bending stresses are significantly reduced. The final optimised cantilever consists of 11 joints and 16 members. Finally, in Figure 15

Pipe bracket
As a structure with a truly three-dimensional load path, the pipe bracket shown in Figure 16 is considered. The design domain has a length, height and width of 120 × 40 × 60 and contains two openings with each of radius of 18 for two pipes passing through the domain. Within the openings each pipe is supported at four points, applying at each support a vertical force of F = 100. The four vertical outer edges of the design domain are chosen as fixed. The finite element discretisation consists of 120 × 40 × 60 linear hexahedral elements. To represent the two openings in the finite element model the Young's modulus of the elements within the void region is prescribed with E min . The maximum volume fraction of the optimised structure is prescribed to be V f = 0.1 and the density filter radius is chosen as R = 3. Figure 17a depicts the optimised bracket with only the voxels above a relative density η = 0.5. The compliance of the optimised voxel model is J(ρ) = 2.47602. As easily recognisable from Figure 17a, the optimised structure is a combination of an arch-like and a cable-like structure with vertical tension members between the two. The skeletonisation algorithm yields after 6 × 6 voxel removal steps the voxel chain skeleton in Figure 17b. The skeleton is converted into the structure depicted in Figure 18a by first identifying the 42 joints and then connecting them with 50 straight members. Alternatively, the 50 members can be represented with curved Bézier segments, which would lead to a structural frame with curved beam members. This has not been pursued here because such members may lead to an increase in manufacturing costs and are usually not preferred. Initially, all frame members are assumed to have the same diameter d = 6.296, giving a total volume of V = 0.1V.
According to Figure 19b, the frame with beams of uniform diameter has a compliance J(d, x) = 4.69409, which is significantly higher than J(ρ) = 2.47602 of the voxel model. After the first size optimisation step the compliance is reduced to J(d, x) = 3.03554 and the subsequent layout optimisation step to J(d, x) = 2.26858. In layout optimisation the beams are constrained not to move inside the volume occupied by the two cylindrical pipes. The intersection between a beam and a cylinder is determined by integrating the (positive) signed distance of the cylinder along the beam. It is straightforward to differentiate this constraint integral with respect to nodal coordinates of the beam. The intersection constraints for all the beams are aggregated with the Kreisselmeier-Steinhauser (K-S) function [42]. Several more steps of size and layout optimisation do not lead to a significant reduction in the cost function. Notice that the obtained final compliance J(d, x) = 2.26683 is lower than the compliance J(ρ) = 2.47602 of the topology

Frame-supported plate
In structural design it is very common to combine a plate-or shell-like skin structure with a frame-like support structure. In such composite structures the skin can be modelled as a standard plate or shell structure, and topology optimisation is applied only in the remaining part of the design domain. The design domain shown in Figure 21 combines a thin plate-like skin with a solid bottom part to be optimised. The entire domain is of size 120 × 80 × 20 and contains a small opening of size 40 × 40 × 6. The plate is of thickness 3 and is subjected to a uniform pressure load of 1. The four corners of the design domain are vertically supported. Due to symmetry only a quarter of the design domain is considered in topology optimisation. The finite element discretisation consists of 60 × 40 × 23 linear hexahedral finite elements. Across the height of the design domain there are 6 elements with a size of 1 × 1 × 0.5 in the plate-like top part and 17 unit elements in the bottom part. To represent the small opening in the finite element model the Young's modulus of the relevant elements is prescribed as E min .
In topology optimisation, the maximum volume fraction is prescribed to be V f = 0.25 and the density filter radius is chosen as R = 1.5. The optimised structure with only the voxels above a relative density η = 0.55 is shown in Figure 22a and its voxel chain skeleton in Figure 22b. Although the voxels in the top part of the design domain are present during optimisation and skeletonisation, they are tagged as non-removable. The voxel chain skeleton is converted to the frame structure shown in Figure 23a. The structural model also contains a top plate not shown in Figure 23. The plate is modelled as a Kirchhoff-Love plate and is discretised with 4 × 4 quadrilateral elements and Catmull-Clark subdivision basis functions. The frame has 124 joints and 136 members. Initially, all the members are assumed to have the same diameter of d = 5.010 which gives a total frame volume of V = 4800. The coupling between the plate and the frame is achieved with Lagrange multipliers as described in [43]. The joints between the plate and frame can transfer forces but not moments.
According to Figure 24b, the compliance of the frame-supported plate with members of uniform diameter is J(d, x) = 1210.01, which reduces after several sequential size and layout optimisation steps to J(d, x) = 322.967.
During the layout optimisation, the positions of the joints between the frame and plate are fixed. For comparison, the compliance of the topology optimised voxel model was J(ρ) = 364.596. The final optimised frame structure contains 56 joints and 64 members, see Figure 23b. Its solid CAD model and faceted triangular STL mesh exported from FreeCAD are depicted in Figure 25.

Conclusions
The automated conversion of topology optimised geometries into compact parametric CAD models is a missing link in the wide-scale industrial adoption of optimisation. We introduced a fully automated workflow for synthesising a structurally-sound parametric CAD model from a voxelised solid-void binary image obtained with topology optimisation. Currently, there is no commercial design system, such as Altair Inspire, Abaqus CAE or similar, which offers such a fully automated CAD model generation capability. Homotopic skeletonisation of three-dimensional images is a key component of our proposed workflow, which is an extensively well-studied topic in digital topology. The resulting voxel chain skeleton is interpreted as an undirected graph and processed with standard graph algorithms to yield a frame model. Both the skeletonisation and graph algorithms operate on topology rather than geometry and are not subject to floating-point errors, which makes them exceedingly robust. The so-obtained structural frame models are not optimal because the skeletonisation and graph algorithms do not take into account any mechanical considerations. However, as shown in the numerical examples, after several alternating steps of size and layout optimisation, the structural frames can achieve a similar cost function value like the original topology optimised geometry. The frame model is converted into a compact CAD solid model by recursively combining primitive shapes using boolean operations. Although we used in this paper only cylinders and spheres as shape primitives, it is straightforward to use other primitives.
The generated parametric CAD solid model and its compact CSG tree representation can be edited to refine the optimised design or to combine it with other components into a product. It is worth reminding that in industrial practice, as supported by user studies [44], geometries obtained from optimisation are often the starting point for design exploration and not the endpoint. This usually requires the ability to edit the optimised geometries in a CAD system. A high-level compact representation is also useful for taking into account geometric manufacturing constraints, like minimum thickness or slope or length of the members [26]. Beyond geometry editing, the frame model may also have advantages in structural analysis of the CAD model. The frame model makes it, for instance, easier to check the structure for buckling and inelastic deformations. As a low-order model, it can be analysed orders of magnitude faster than a three-dimensional solid model. This opens up the possibility of instant analysis of CAD models as it has become recently available in, e.g., Ansys Discovery Live or Creo Simulation Live.
In closing, one particular extension of the presented approach is worth mentioning. We assumed that the topology optimised voxel model can be approximated with tubular members. However, in three-dimensional optimisation problems, large flat or curved surfaces may also appear during optimisation. The introduced skeletonisation algorithm in Section 3 will reduce even those into a network of one voxel thick chains. As suggested in Lee et al. [4], the algorithm can be modified to preserve one-voxel thick surfaces which appear during skeletonisation. This is achieved by not deleting the so-called surface points/voxels, which are defined based on the arrangement of solid voxels in an octant. In Figure 26 the topology optimised pipe bracket has been skeletonised with the modified skeletonisation algorithm. The curved surface at the bottom is preserved. In the structural model, this surface can be modelled as a thin-shell and in the CAD model it can be easily converted into a solid by providing a thickness. Currently, we do not have an automated process to extract a shell surface from the skeletonised voxel model, but this seems to be possible and is a promising direction for future research.

Appendix A. Member stiffness matrix
We model the members of the structural frame as Timoshenko beams. As mentioned, the members are assumed to be straight and have uniform cross-sections along their lengths. Hence, each member can be approximated with a single beam finite element and its stiffness matrix can be derived with standard structural analysis techniques from undergraduate textbooks, see e.g. [45]. The stiffness matrix we use takes into account axial stretch, transversal shear, bending and torsion effects. The local stiffness matrix K l i of a beam i is derived in a local coordinate system in which the beam axis is aligned with the x l i axis. Each of the two end nodes of the beam have three displacement and three rotation degrees of freedom, see Figure A.27. To ease the notation, we write in the following for the local coordinate axes only x, y and z.
The element stiffness matrix of the beam is composed out of independent axial, bending and torsion stiffness matrices. The axial stiffness matrix corresponding to the displacement degrees of freedom u (1) x and u (2) x is given by where E, A and L denote the Young's modulus, cross-sectional area and length. The bending stiffness matrix corresponding to the displacement degrees of freedom u (1) y , u (1) z , u (2) y and u (2) z and the rotational degrees of freedom θ (1) y , θ (1) z , θ (2) y and θ (2) z is given by 12EI y (1+b y )L 3 6EI y (1+b y )L 2 0 (4+b y )EI y (1+b y )L 0 (4+b z )EI z (1+b z )L where I y and I z are the second moments of area about the y and z axis. The two dimensionless factors b y and b z take into account the transversal shear and are defined according to b = 12EI/(κGAL 2 ), with G the shear modulus and κ the shear correction factor. In case of a circular cross section we have I y = I z = I, b y = b z = b and κ = 0.9. Note that setting b y = b z = 0 in (A.2) gives the stiffness matrix of an Euler-Bernoulli beam. Finally, the torsional stiffness matrix corresponding to the rotational degrees of freedom θ (1) x and θ (2) x is given by where J is the torsion constant of the cross section. The local stiffness matrix K l i of an element i is obtained by suitably assembling the three stiffness matrices (A.1), (A.2) and (A.3) into a 12 × 12 matrix. This matrix in the local coordinate system x l i , y l i and z l i is transformed into a stiffness matrix K i in the global coordinate system x, y and z according to (11).