Explicit Reconstruction and Shape Optimization of Topology Optimization Results with Mechanical Performance Preservation

1  Introduction

Topology optimization integrates finite element methods and mathematical optimization algorithms to generate optimal material layouts at the conceptual design stage, and has found wide applications across various engineering fields [1]. According to the manner in which optimization results are presented, topology optimization methods can be classified into implicit approaches, such as density-based and level set methods [2–4], and explicit approaches, such as moving morphable components or void methods [5,6]. However, significant challenges remain when applying topology optimization results to downstream applications. Due to mesh discretization and density penalization mechanisms, density-based topology optimization often produces intermediate-density regions and jagged boundaries. These artifacts lead to blurred structural contours and may result in localized stress concentrations [7]. The level set method represents optimization results implicitly, which hinders direct integration with computer aided design (CAD) systems [8]. Similarly, the results of the moving morphable components/void method rely on Boolean operations among components or void and lack a unified parametric description [9,10]. Among these approaches, density-based methods, such as the solid isotropic material with penalization (SIMP), have been widely adopted owing to their conceptual simplicity and numerical robustness [11]. In light of the aforementioned limitations, extensive research efforts have been devoted to developing post-processing techniques for interpreting the optimization results.

In terms of boundary recognition and smoothing of topology optimization results, considerable research has focused on extracting well-defined structural boundaries and enhancing geometric regularity. Hsu et al. [12] extracted contour lines from cross-sections of optimization results and employed sweeping techniques to reconstruct geometric models. Chu et al. [13] utilized support vector machine-based classification to distinguish different density regions, thereby obtaining clear structural boundaries. Tang and Chang [14] transformed the boundaries of topology optimization results into smooth, parameterized B-spline curves and surfaces. Wu et al. [15] performed vectorized boundary modeling based on Freeman chain codes and, by introducing boundary curvature parameters and finite element analysis, automatically identified geometric features and extracted fitting control points, enabling regularized boundary reconstruction and shape optimization. Li et al. [16] proposed a boundary density evolution method, in which unpenalized interpolation and density filtering were used to obtain clear topologies, followed by boundary post-processing based on nodal strain-energy-driven level sets to generate smooth contours. Swierstra et al. [17] employed a radial basis function-based level set approach to automatically extract geometric boundaries and combined it with the finite cell method for high precision boundary shape optimization, significantly improving smoothness. Li et al. [18] developed an efficient boundary smoothing strategy for bi-directional evolutionary structural optimization (BESO) results using pre-constructed lookup tables, generating smooth topologies while strictly preserving volume and key geometric features. Ježek et al. [19] proposed a geometric extraction framework based on signed distance functions and radial basis functions, achieving highly smooth boundaries while preserving topology and volume, thereby enhancing mechanical performance and manufacturing compatibility. Lin et al. [20,21] treated topology optimization results as images and, through geometric feature matching combined with artificial neural networks, automatically identified and fitted holes using parameterized templates, enabling a fully automated transition from topological layouts to shape optimization models. Yildiz et al. [22] employed neural-network-driven image processing techniques to automatically map holes in topology optimization results to predefined, manufacturing-oriented geometric features, thereby achieving an efficient conversion from conceptual layouts to optimizable and manufacturable models. Gamache et al. [23] developed a dedicated skeletonization algorithm for topology optimization results, which converts the density field into a truss-like skeleton while preserving mechanical connectivity, significantly enhancing the interpretability of low-order numerical results in terms of higher-level engineering concepts.

From the perspective of downstream design and manufacturing, extensive efforts have been made to reconstruct topology optimization outputs into parameterized and CAD-ready geometric representations. Joshi et al. [24] proposed an automated reconstruction pipeline that converts voxel-based topology optimization results into NURBS surfaces by combining Dual Contouring with C1-continuous surface fitting, enabling direct import into CAD systems. Liu et al. [25] achieved the automatic generation of high-fidelity yet low-complexity parametric CAD models through skeleton-guided boundary enhancement and curvature-adaptive B-spline fitting. Chacón et al. [26] developed an automatic conversion framework that transforms 2D topology optimization images into IGES-format B-spline models, significantly improving manufacturability within CAD/computer aided manufacturing (CAM) environments. Koguchi and Kikuchi [27] extracted closed iso-surfaces using the marching cubes algorithm and reconstructed them with bi-quartic surface splines, preserving key geometric features while ensuring boundary smoothness for parametric CAD modeling. Yoely et al. [28] embedded explicit B-spline parameterizations directly into the topology optimization process, incorporating hole-size and curvature constraints to ensure manufacturability from the outset. Hsu and Hsu [29] realized a fully automated conversion of three-dimensional topology optimization results into smooth CAD models through sectional decomposition and B-spline contour fitting, combined with density filtering. Cuilliére et al. [30] systematically investigated the integration bottlenecks encountered when transitioning from CAD to topology optimization and back to CAD for unstructured three-dimensional meshes, and compared the applicability and limitations of threshold-based and iso-density surface extraction methods. Wen et al. [31] combined marching cubes with sparse curve fitting to automatically convert topology optimization results into editable CAD models, enabling interactive design refinement via a Rhino-based plugin, particularly for geometries with fine scale features. Larsen and Jensen [32] employed predefined two-dimensional template fitting and sweeping-based modeling strategies to achieve semi-automatic conversion of topology optimization results into parametric CAD models, allowing a controllable trade-off between fitting accuracy and feature complexity. Bacciaglia et al. [33] developed a feature preserving automatic mesh smoothing algorithm that freezes critical geometric regions, preventing hole loss and excessive shrinkage typically associated with conventional vertex-based smoothing, and is suitable for industrial scale topology optimization post-processing. Ren et al. [34] proposed an adaptive multi-resolution CAD reconstruction framework, which integrates instant meshes quadrangulation, improved harmonic mapping, and adaptive sampling to convert topology optimization results into sparse-control-point, low-patch, and highly editable NURBS boundary representation models. By initializing the parameter domain using geodesic distances and adopting a multi-resolution strategy, the framework significantly reduces model complexity.

Although various strategies have been proposed to interpret topology optimization results and reconstruct CAD-compatible geometries, several limitations remain. For topology optimization problems formulated with compliance minimization, most existing reconstruction approaches rely on geometric or graphics-driven techniques and pay limited attention to changes in structural performance, which often leads to reconstructed models that deviate from the volume constraints and performance objectives obtained in the topology optimization stage. This study proposes a two-stage framework that explicitly reconstructs topology optimization results into parameterized geometries while maintaining the structural performance obtained in the topology optimization stage. The flowchart of the two-stage framework is illustrated in Fig. 1. In the first stage, the contour points are sequentially extracted, smoothed, filtered, and interpolated to obtain a parametric geometric representation. In the second stage, the initial geometry is further optimized by modifying the boundary control points, with the optimization driven by finite element analysis, sensitivity caculation, and a gradient-based solver. Finally, topology optimization results are converted into parametric models with mechanical performance preserved.

Figure 1: Flowchart of the two-stage framework with main steps.

To clarify the differences between the proposed framework and representative existing methods, a comparison of key properties of the reconstructed results is summarized in Table 1.

The remainder of this paper is organized as follows. Section 2 introduces the fundamental principles of density-based topology optimization. Section 3 presents the proposed geometry reconstruction procedure. The shape optimization procedure is described in Section 4. Section 5 focuses on the analysis and discussion of the results produced by the proposed framework. Representative numerical examples, together with a practical engineering case of a quadcopter frame, are used to further assess the effectiveness and applicability of the framework. Conclusions are presented in Section 6. Appendix A presents the parametric sensitivity and robustness analysis of structural performance and shape, and Appendix B reports the computational cost analysis for the benchmark case.

2  Topology Optimization

Topology optimization is a class of structural design methodologies grounded in mathematical programming, in which structural design problems are formulated as computable models by optimally distributing material within a design domain to achieve targeted performance objectives. As a representative mathematical optimization problem, topology optimization relies on a well-defined objective function, appropriate constraint conditions, and sensitivity information of the design variables with respect to both the objective and the constraints. Together, these components constitute the theoretical foundation and critically influence algorithmic convergence and the reliability of the obtained results.

The basic principle of density-based topology optimization for compliance minimization problems can be summarized as follows. The design domain is first discretized into a finite element mesh, and an initial density value is assigned to each element, which governs its effective mechanical properties. Using mathematical optimization algorithms, the contribution of each element to the structural stiffness is evaluated. During the optimization process, the densities of inefficient elements are progressively reduced, while those of efficient elements are increased, leading to an optimized material distribution. As a result, the structural topology evolves and the topology optimization process is completed. The corresponding optimization formulation is given as

Find:ρ=(ρ1,…,ρn)∈Rn,min:C=FTU,s.t.{K(ρ)U=F,V(ρ)V0≤V¯,ρe∈[0,1],e=1,…,n.(1)

In Eq. (1), ρ denotes the vector of design variables, with each component ρe representing the density of the e-th element. The objective function C represents the structural compliance. F denotes the external load applied to the structure. U is the nodal displacement vector. K(ρ) is the global stiffness matrix, which depends on the element densities. V(ρ) denotes the volume of the optimized structure, V0 denotes the volume of the design domain, and V¯ denotes the volume fraction constraint.

In topology optimization, the relationship between the elastic modulus of each element and its density is commonly established using the SIMP scheme, which is expressed as

Ee=ρeβE0,(2)

where Ee and E0 denote the elastic modulus of the e-th element and the solid material, respectively. β is the penalization factor, which is typically set to 3 for compliance minimization problems. To ensure numerical stability, a lower bound of ρe=10−3 is commonly imposed [35], thereby avoiding singularity of the global stiffness matrix.

The design variables are updated using gradient-based optimization algorithms, with the optimality criteria (OC) method adopted as the solution strategy. The sensitivities of the objective function and constraints with respect to the design variables are given by

∂C∂ρe=−βρeβ−1ueTk0ue,∂V(ρ)∂ρe=Ve.(3)

Here, ue denotes the displacement vector of e-th element extracted from the global displacement vector U, k0 is the elemental stiffness matrix of the solid material with unit Young’s modulus, and Ve denotes the volume of the e-th element. Accordingly, the scalar quantity ueTk0ue corresponds to the strain energy contribution of element e under the reference material properties, and it measures the contribution of that element to the overall structural compliance.

Fig. 2 shows the topology optimization results of a cantilever beam subjected to a concentrated load at the upper-right corner, obtained using the classic 99-line code [35]. The design domain is discretized into 100×50 elements, with a target volume fraction of 50%. The resulting optimized design achieves a compliance of 72.1995 with an actual volume fraction of 0.4999. However, the structural boundaries shown in Fig. 2b are diffuse and exhibit pronounced jaggedness due to the density-based representation and the underlying mesh discretization. As a result, the optimized design lacks explicit and smooth geometric boundaries, making it unsuitable for direct CAD editing and manufacturing. These inherent geometric deficiencies of topology optimization outputs necessitate a geometry reconstruction procedure to obtain smooth, explicit, and manufacturable structural representations.

Figure 2: Topology optimization results of a cantilever beam. (a) Initial structure; (b) Optimized structure.

3  Geometry Reconstruction

The contour points of both the inner holes and the outer boundary are first extracted from the density field using the marching squares algorithm. Since the extracted points are typically dense and irregular, direct interpolation would result in NURBS-based boundaries with spurious irregularities and an excessive number of control points, thereby reducing geometric controllability. To alleviate these issues, smoothing and filtering operations are applied to the extracted point sets. Finally, NURBS interpolation is performed to produce a smooth and explicitly parameterized geometry suitable for shape optimization.

3.1 Contour Points Extraction

The marching squares algorithm [36] extracts iso-density contours by examining the values within each element and locating contour segments according to a specified threshold τ. Fig. 3a illustrates an example density field, which serves as the input for contours extraction. When the scalar values at the two endpoints of a line segment, denoted by v1 and v2, satisfy the condition v1<τ<v2, the intersection position s along the segment is computed via linear interpolation

s=x1+τ−v1v2−v1(x2−x1),(4)

where x1 and x2 denote the coordinates of the two vertices of the segment.

Figure 3: Illustration of contour points extraction. (a) Illustrative density field; (b) Direct contour points extraction; (c) Extraction after density padding.

However, directly applying the marching squares algorithm may result in incomplete extraction of the outer boundary, as illustrated in Fig. 3b. This issue arises from incomplete sign changes across the threshold in scalar field adjacent to the bounding box. To accurately preserve the overall dimensional characteristics of the structure while automatically extracting the complete outer boundary, a virtual element padding strategy is introduced. Specifically, the material distribution near the design domain boundary is typically either fully solid (ρ≈1) or void (ρ≈0). By extending the optimized density field with an additional layer of virtual elements assigned the value 2τ−1, the scalar values at the two sides of a boundary become v2=1 and v1=2τ−1. Substituting v2 and v1 into Eq. (4) yields

s=x1+τ−(2τ−1)1−(2τ−1)(x2−x1)=x1+1−τ2(1−τ)(x2−x1)=12(x1+x2).(5)

In this expression, x1 and x2 denote the locations of the virtual element and the boundary element in the optimized density field, respectively. The term 12(x1+x2) represents the midpoint between the two elements, ensuring that the intersection point is precisely aligned with the boundary of the design domain. As a result, the complete outer boundary can be correctly extracted, as shown in Fig. 3c.

In summary, the proposed strategy enables contours extraction at arbitrary iso-density levels within inner hole regions, while guaranteeing accurate alignment with the design domain and continuity of the outer boundary contours.

3.2 Contour Points Smoothing and Filtering

Specifically, consider a two-dimensional sequence of points {Pi}i=1N associated with a given contour, where Pi=(xi,yi). To smooth the boundary point sequence, a sliding window of size Nw is introduced. The two endpoints are excluded from the smoothing process and retain their original positions, i.e., P^i=Pi for i=1 and i=N, where P^i denotes the smoothed point. All remaining boundary points are smoothed using a local mean filter strategy to suppress geometric noise

P^i=1|Wi|∑j∈WiPj,Wi=[max{0,i−⌊Nw2⌋}, min{N,i+⌊Nw2⌋+1}],(6)

where Wi denotes the index range of the sliding window and |Wi| is its cardinality. After the points are smoothed, a filtering strategy is further applied to reduce the number of points. Starting from the first point, subsequent points are retained only if their distance from the previously retained one is not smaller than the filter radius rf. In contrast, the outer contour undergoes smoothing only, without applying the distance-based filtering, which results in a relatively dense point distribution to ensure higher geometric fidelity for interpolation.

3.3 Contour Points Interpolation

In this study, a NURBS interpolation method based on global geometric constraints [37] is employed to reconstruct the preprocessed contour points into smooth, explicitly parameterized boundary representations that are well suited for shape optimization and CAD-based design operations. Given a set of points Pi∈R2 (i=0,1,…,N), a p-degree NURBS curve is constructed through three standard steps. The resulting NURBS curve is expressed as

C(u)=∑j=0NNj,p(u)Qj,(7)

where Qj denote the control point to be determined, and Nj,p(u) are the p-degree NURBS basis functions defined over the knot vector [u0,u1,…,um]. This formulation exactly interpolates the point sequence by solving a system of linear equations.

In the first step, a parameter value ui is assigned to each data point Pi, thereby establishing a correspondence between the point sequence and the parameter domain. The parameter values are computed using a chord-length-based scheme,

ui={0,i=0,ui−1+|Pi−Pi−1|∑j=1N|Pj−Pj−1|,i=1,…,N−1,1,i=N,(8)

and are subsequently normalized to the interval [0,1].

In the second step, the knot vector is constructed using the averaging method. With end knots of multiplicity p+1, the total number of knots is given by m=N+p+2, and the internal knots are defined as

uk={0,k=0,…,p,1p∑j=k+1−pkuj,k=p+1,…,m−p−1,1,k=m−p,…,m.(9)

In the final step, the control points are determined by enforcing the interpolation conditions C(ui)=Pi, which leads to the following linear system

AQ=P,(10)

where Ai+1,j+1=Nj,p(ui) is the coefficient matrix, Q=[Q0,…,QN]T is the control point vector and P are the points for interpolation. Solving this system yields a uniquely determined set of control points, resulting in a smooth curve with Cp−1 continuity.

4  Shape Optimization

Since the geometry reconstruction procedure is sensitive to user-defined parameters such as the extraction threshold τ, the smoothing window size Nw, and the filtering radius rf, the selection of these parameters does not explicitly account for structural mechanical performance. As a result, the mechanical performance of the reconstructed geometry may deviate from the volume fraction and the compliance obtained by the topology optimization, thereby motivating the shape optimization step to recover mechanical performance. Moreover, because the structural topology has already been established during the topology optimization stage, preserving topological homeomorphism [38], which maintains connectivity and branching characteristics under continuous geometric evolution, becomes a key requirement to ensure consistency between the optimized shape and the underlying topology.

4.1 Optimization Formula and Design Variables

As shown in Fig. 4, the design domain in the shape optimization stage is defined as the area enclosed by the outer boundary, while the solid domain is obtained through a Boolean operation between the design domain and the void domain. The shape optimization problem is formulated as

Find:r=(r1,…,rm)∈Rm,min:C=FTU,s.t.{K(r)U=F,V(r)V0≤V¯.δ=0.(11)

Figure 4: Structural shape generation via Boolean subtraction of inner holes from the design domain.

Compared with the topology optimization formulation in Eq. (1), both problems share the same objective of minimizing structural compliance. However, the design variables are changed from the element density ρ to the polar radius vector r of the inner boundary control points. In addition to the volume fraction constraint, a topological homeomorphism constraint defined by δ=0 is imposed to preserve topological equivalence, which will be discussed in Section 4.3.

Taking the contour in Fig. 5a as an example, boundary self-intersection may arise during shape optimization if the control points are directly manipulated in Cartesian coordinates, as shown in Fig. 5b. To avoid this issue, the Cartesian coordinates are transformed into polar coordinates according to

ri=(xi−xc)2+(yi−yc)2,θi=atan2⁡(yi−yc,xi−xc),(12)

as illustrated in Fig. 5c. Here, (xi,yi) denote the Cartesian coordinates of the control points, while

xc=1n∑i=1nxi,yc=1n∑i=1nyi,(13)

denote the geometric center of the control points for a single contour. The function atan2⁡(⋅,⋅) represents the four-quadrant inverse tangent. To ensure consistency in the angular representation, the angle θi is normalized to the interval [0,2π) as

θi←(θi+2π)mod(2π).(14)

Figure 5: Polar radius-based design variables. (a) Cartesian coordinates; (b) Self-intersection; (c) Polar coordinates; (d) Shape adjustments.

Thus, the control points can be represented by the triplet

(xc, yc, {(θi, ri)}i=1n).(15)

The polar radii of all holes are then assembled to form the design variable vector r=(r1,…,rm)∈Rm, thereby restricting the movement of each control point to its radial direction. Fig. 5d demonstrates the effectiveness of this approach in enabling shape adjustments while preventing boundary self-intersection.

4.2 Level Set Function and FG-FEM

The optimization objective C is solved using finite element method, frequent remeshing of the grids following each geometric update would significantly reduce computational efficiency and pose considerable implementation challenges. Therefore, FG-FEM is adopted for structural response analysis. In this approach, the structure is represented using a level set function (LSF) Φ, which links the geometric description with the analysis model. A LSF is introduced over the design domain Ω⊂Rd (d=2), defined as

{ΦΩ(x)=0,∀x∈∂Ω,ΦΩ(x)>0,∀x∈Ω∖∂Ω,ΦΩ(x)<0,∀x∈D∖Ω.(16)

Here, x denotes an arbitrary point in space, and the structural boundary ∂Ω is defined as the zero level set of ΦΩ. A positive value of ΦΩ(x) indicates that the point lies inside the solid domain, whereas a negative value corresponds to the void domain. The corresponding composite LSF of the solid domain is defined as

ΦΩ(x)=min{{φk}k=1Nv,ΦD},(17)

where Nv denotes the number of inner holes, with φk and ΦD representing the LSFs of the void domain and the design domain, respectively. The construction of the corresponding LSFs is described sequentially.

4.2.1 LSF for Void Domain and FG-FEM

For an arbitrary point x=(x,y), its polar coordinates (θ,r) with respect to the geometric center c=(xc,yc) of the closed contour are defined as

r(x)=(x−xc)2+(y−yc)2,θ(x)=atan2⁡(y−yc,x−xc).(18)

The radial distance from the boundary to the center at a given polar angle θ is denoted by R(θ). This function is obtained by sampling discrete points along the boundary, computing their corresponding polar angles and radii, and fitting the resulting data using a univariate spline function

R(θ)=𝒮({(θi,ri)}i=1n),(19)

where 𝒮(⋅) denotes a univariate spline interpolation operator, and n denotes the total number of sampled boundary points. The fitted function R(θ) is illustrated in Fig. 6a, based on which the LSF φk(x) of the boundary is defined as

φk(x)=r(x)−R(θ(x)).(20)

Figure 6: LSF representation and region classification. (a) Polar-radius function R(θ); (b) LSF values; (c) Categorized regions.

The resulting LSF values are illustrated in Fig. 6b.

According to the relative positions to the boundary, the elements are classified into three categories, as illustrated in Fig. 6c. To determine the material distribution on the fixed grid, the LSF is projected using a regularized Heaviside function, yielding the effective material indicator

H(ΦΩ)={1,ΦΩ>Δ,3(1−λ)4(ΦΩΔ−ΦΩ33Δ3)+1+λ2,−Δ≤ΦΩ≤Δ,λ,ΦΩ<−Δ.(21)

The corresponding function profile is shown in Fig. 7a, where Δ denotes the width of the transition zone and λ represents the minimum relative density, which is assigned as 0.001 to prevent singularity of the stiffness matrix. By substituting φ into the regularized Heaviside function, the projected LSF values are obtained, as shown in Fig. 7b. Subsequently, the equivalent density of each quadrilateral element is computed as the average of the projected LSF values at its nodes, given by

ρe=14∑n=14H(ΦΩ(xne)),(22)

where xne denotes the coordinates of the n-th node of the e-th element. The resulting element density field is illustrated in Fig. 7c. Accordingly, the structural volume defined in Eq. (11) is evaluated by summing the densities of all elements

V(r)=∑e=1Neρe,(23)

where Ne denotes the total number of elements.

Figure 7: Heaviside projection and resulting density field. (a) Heaviside function profile; (b) Projected LSF values; (c) Density field.

Finally, by incorporating the SIMP material interpolation scheme given in Eq. (2), the equivalent elastic modulus of each element and the corresponding stiffness matrix are computed, thereby completing the transformation from the geometric model to the finite element analysis model for structural compliance evaluation.

4.2.2 LSF for the Design Domain

Since the outer contours of topology optimization results are generally non-star-shaped [39], the formulation of Eq. (19) cannot be directly realized, as non-star-shaped geometries do not admit a unique mapping between the polar angle θ and the radial distance r, with a given angle intersecting the boundary at multiple points. Consequently, the polar-radius parameterization introduced in Section 4.1 is no longer applicable for representing such outer contours. Therefore, to enhance the generality and robustness in handling outer contours, this study proposes a sign–distance decoupled strategy, in which the LSF value at a point is decomposed into two independent components, namely a sign component and a magnitude component. The sign of the LSF reflects whether the point is located inside or outside the boundary, and the magnitude corresponds to the distance from the point to the boundary. Based on this decomposition, the inside–outside relationship and the distance of a point to a boundary are evaluated in a decoupled manner. Specifically, the sign of the LSF at a grid node is determined using the winding number method [40], while the distance is computed through a geometric projection onto the boundary. This strategy avoids the limitations of polar representations and enables reliable treatment of non-star-shaped outer boundaries. Fig. 8 illustrates a representative classification example for fixed grid nodes based on the proposed decoupled strategy. The sign field, shown in Fig. 8a, indicates the inside–outside relationship of nodes with respect to the outer boundary. The corresponding distance field is presented in Fig. 8b. Owing to the use of the Heaviside projection in Eq. (21), distance calculations are restricted to the nodes located within the width Δ to the boundary, while nodes far from the contour are directly assigned a value of 1. This localized distance evaluation significantly reduces the computational cost without compromising the accuracy of the LSF near the boundary.

Figure 8: Computed results of the sign field and distance field. (a) Sign field; (b) Distance field.

Following this decomposition, the sign field and the distance field are multiplied pointwise to construct ΦD associated with the design domain. The resulting ΦD is then combined with the φk terms associated with all inner holes and incorporated into Eq. (17), yielding ΦΩ for the solid domain.

4.3 Topological Homeomorphism Constraints

In the optimization model, the minimum distance required to prevent interference between adjacent holes cannot be prescribed a priori. Consequently, directly imposing simple bound constraints [rmin,rmax] on design variables would restrict the feasible design space and result in overly conservative design limits. To address this issue, an interference detection and elimination strategy is proposed, in which potential boundary interference during geometric evolution is quantitatively identified and resolved, thereby enforcing the topological homeomorphism constraint while allowing full exploration of high-performance structural shapes.

Building on the FG-FEM, consider a domain discretized into Ny×Nx finite elements. Let φk denote the LSF associated with the k-th hole, constructed from the boundary curve Ck and evaluated at an arbitrary grid node (X,Y). According to Eq. (20), φk(X,Y)<0 if and only if the node lies inside Ck. The values of the LSFs for Nv holes evaluated on the discretized grid nodes are stored in a three-dimensional array φ∈R(Ny+1)×(Nx+1)×Nv, such that

φi,j,k=φk(Xj,Yi),i=1,…,Ny+1, j=1,…,Nx+1, k=1,…,Nv.(24)

The coverage multiplicity χ on the grid node (Xj,Yi) is then defined as

χi,j=∑k=1NvI(φi,j,k<0),(25)

where I(⋅) denotes the indicator function. A value χi,j≥2 means that the node is simultaneously covered by at least two holes, signifying the occurrence of geometric interference. The total interference measure, defined in Eq. (11) as the sum of coverage multiplicities over all nodes, is given by

δ=∑i=1Ny+1∑j=1Nx+1χi,jI(χi,j≥2).(26)

The quantity δ serves as a measure of geometric interference and is employed to guide the interference elimination procedure. Specifically, δ>0 indicates the occurrence of boundary interference, whereas δ=0 indicates that no boundary interference occurs, thereby enforcing the topological homeomorphism constraint throughout the structural evolution process.

Fig. 9 illustrates a representative example in which interference between inner holes violates the topological homeomorphism constraint, together with the corresponding interference elimination process. At the initial phase, geometric interference is detected, with 14 nodes concurrently covered by two holes. During the interference elimination process, the polar radii gradually contract, resulting in a progressive decrease in the interference values. Eventually, the interference is completely eliminated, the polar radii stabilize, and the topological homeomorphism constraint is satisfied (δ=0).

Figure 9: Topological homeomorphism constraints and interference handling.

4.4 Sensitivity Analysis

The sensitivity of the compliance with respect to the design variable ri is derived as

∂C∂ri=∂FT∂riU+FT∂U∂ri.(27)

Since the external load F is independent of the design variables, and recalling the equilibrium equation KU=F, differentiation with respect to ri yields

∂K∂riU+K∂U∂ri=0,(28)

from which the displacement sensitivity is obtained as

∂U∂ri=−K−1∂K∂riU,(29)

substituting Eqs. (29) into (27), the compliance sensitivity with respect to ri is obtained as

∂C∂ri=FT∂U∂ri=−UT∂K∂riU=−βE0∑e=1Ne(ρeβ−1∂ρe∂riueTk0ue).(30)

The sensitivity of the volume with respect to ri is expressed as

∂V∂ri=∑e=1Ne∂ρe∂ri.(31)

Based on Eq. (22), the sensitivity of the element density ρe with respect to ri can be written as

∂ρe∂ri=14∑n=14∂H(ΦΩ(xne))∂ri,(32)

which is difficult to evaluate analytically due to the implicit dependence of the LSF on the design variables.

The sensitivity of the interference measure δ with respect to the design variable ri is not analytically differentiable. Thus the sensitivities of both the interference measure and the element density are computed numerically using the central finite difference scheme, which is widely adopted in sensitivity analysis for its simplicity and numerical stability [41]. The finite difference approximation is applied as follows:

∂δ∂ri≈ΔδΔri=δ(ri+Δr)−δ(ri−Δr)2Δr,∂ρe∂ri≈ΔρeΔri=ρe(ri+Δr)−ρe(ri−Δr)2Δr,(33)

where Δr denotes the finite difference step size used in the sensitivity evaluation.

4.5 Shape Optimization Procedure

Fig. 10 illustrates the detailed procedure of the shape optimization framework. The optimization process starts with the initialization of the design variables and algorithmic parameters. Since the reconstructed geometry is obtained from topology optimization results and the boundary is initialized in the intermediate density region, the subsequent optimization process may suffer from premature convergence. To alleviate this issue, small random perturbations are introduced during the initialization stage. Specifically, for each hole, a set of random noise values with the same dimension as the polar radii variables is generated from a uniform distribution within the interval [−2Δr,2Δr], where Δr is the step size used in the central finite difference scheme (Eq. (33)). The perturbations are added to the polar radii variables to slightly modify the initial design configuration.

Figure 10: Shape optimization flowchart.

After initialization, the density field of the design domain is evaluated based on the current boundary representation. The interference measure δ is then calculated to determine whether boundary interference occurs. If δ≠0, an interference elimination procedure is activated. In this stage, the sensitivity of δ with respect to the design variables is computed, and the design variables are updated to reduce the interference until δ=0. This process guarantees that the reconstructed geometry preserves the topological homeomorphism of the topology optimization result.

Once the interference has been removed (δ=0), the algorithm proceeds to the shape optimization stage. The sensitivities of the objective function (compliance C) and the volume constraint V are computed, and the design variables are updated using the Method of Moving Asymptotes (MMA) [42]. After each update, the density field of the solid domain and the interference measure δ are recalculated to verify whether new boundary interference has been introduced during the optimization process.

The entire procedure is embedded within an outer iterative loop. The optimization proceeds until convergence is achieved or the maximum number of iterations is reached. Convergence is evaluated using two criteria: the difference between the current constraint value and its prescribed value, and the maximum change in the design variables between two consecutive iterations. The optimization terminates when both quantities fall below their prescribed thresholds, namely 10−4 for the constraint deviation and 10−2 for the design variable changes.

5  Results and Discussion

This section presents a comprehensive discussion of the proposed geometry reconstruction and shape optimization framework through three categories of examples. First, the topology optimization results shown in Fig. 2b are adopted as a benchmark case to systematically investigate the effectiveness of the proposed framework. Second, to further assess the robustness and generalization capability of the framework, three classical problems, namely the cantilever beam, the Michell beam, and the half-MBB beam with distinct geometric scales are considered as numerical examples. Finally, to demonstrate the applicability of the proposed framework to practical engineering design, a quadcopter frame is presented as an engineering example.

5.1 Benchmark Case

The contour points shown in Fig. 11 are extracted from the topology optimization results using a threshold value of τ=0.6. With the proposed virtual element padding strategy, a closed and continuous outer contour is obtained. Building upon this capability, the strategy enables the extraction of contours at arbitrary threshold levels, while ensuring the outer contour remains properly aligned with the design domain.

Figure 11: Contour points extraction.

The contour points obtained after smoothing and filtering are shown in Fig. 12, where the smoothing window size and filtering radius are set to Nw=3 and rf=3, respectively. Through the sparse sampling procedure, the total number of contour points is reduced from 445 to 99. In addition, the outer contour points are locally densified near the corner regions of the design domain to prevent excessive curvature distortion when high-order NURBS curves undergo sharp turning. As a result, the contour points exhibit enhanced smoothness and geometric consistency, providing a stable and well-conditioned input for subsequent interpolation. For comprehensive investigations into the influence of the reconstruction parameters (τ, Nw, and rf) on mechanical performance, the reader is referred to Appendix A.

Figure 12: Contour points smoothing and filtering.

The interpolated NURBS boundaries and corresponding control points are shown in Fig. 13, where a third-order (p=3) closed interpolation is achieved by appending the midpoint of the endpoints to both ends of the filtered point sequence, resulting in a total of 105 design variables. Compared with the topology optimization results, the reconstructed model provides a clear and parametric description. The main geometric features of both the outer and inner contours are accurately captured, yielding smooth and continuous boundaries that help reduce stress concentrations. Moreover, the sparse and well organized distribution of control points leads to a reduced-dimensional design variable space, which significantly reduces the computational cost of shape optimization.

Figure 13: Contour points interpolation and resulting boundary.

Through the shape optimization step, the inner holes are systematically adjusted. As shown in Fig. 14, the optimized geometry exhibits only moderate changes relative to the initial reconstruction, indicating that the shape optimization primarily serves to recover mechanical optimality and constraint consistency rather than to alter the overall structural layout.

Figure 14: Reconstructed geometries before and after optimization.

The shape optimization history and quantitative results shown in Fig. 15 provide clear insight into the convergence behavior and effectiveness of the proposed optimization strategy. As shown in Fig. 15a, nonzero interference values are detected during the early iterations, which trigger the interference elimination procedure. During this phase, the optimization temporarily prioritizes the reduction of interference rather than the minimization of structural compliance. Once the interference is fully eliminated (δ=0), the algorithm automatically switches to the shape optimization procedure, during which the volume fraction rapidly converges to the target value and remains tightly controlled. Meanwhile, the compliance decreases smoothly, indicating stable convergence without reintroducing geometric interference, with the process terminating after 134 iterations. This behavior demonstrates that the volume constraint is effectively enforced in the early stage of the optimization, after which the algorithm focuses on further reducing structural compliance. Fig. 15a shows that the reconstructed geometry prior to shape optimization exhibits deviations of 3.8% in compliance and 3.2% in volume fraction relative to the topology optimization results, indicating that the geometry reconstruction process introduces noticeable performance discrepancies. After the shape optimization stage, these deviations are reduced to 0.3% and 0%, respectively, demonstrating the effectiveness of the proposed framework in preserving mechanical performance. Fig. 15b visualizes the interference reduction process, showing that the initially intersecting boundaries are gradually separated to satisfy the topological homeomorphism constraints.

Figure 15: Shape optimization history, mechanical performance comparision and interference elimination behavior. (a) Shape optimization history and corresponding mechanical performance comparison; (b) Results of the interference elimination process.

5.2 Numerical Examples

Fig. 16 shows the cantilever beam, Michell beam, and half-MBB beam, each with a distinct geometric scale. The corresponding results of topology optimization, geometry reconstruction, and shape optimization are presented in Fig. 17. It should be noted that, in the Michell beam case, the final optimized structure exhibits a certain degree of asymmetry. This behavior arises from the asymmetric distribution of control points generated during the geometry reconstruction stage. Since the polar radii are adopted as design variables, each control point is constrained to move only along its radial direction, and symmetry in the final optimized structure is therefore not always guaranteed.

Figure 16: Numerical examples. (a) Cantilever beam; (b) Michell beam; (c) Half MBB beam.

Figure 17: Topology optimization results and reconstructed geometries. (a) Topology optimization of Cantilever beam; (b) Reconstructed geometries of Cantilever beam; (c) Topology optimization of Michell beam; (d) Reconstructed geometries of Michell beam; (e) Topology optimization of half-MBB beam; (f) Reconstructed geometries of half-MBB beam.

The detailed mechanical performance values are summarized in Table 2, and the corresponding relative deviations with respect to the topology optimization results are illustrated in Fig. 18. The cantilever beam exhibits the largest discrepancies, with a compliance increase of up to 25.6% and a volume fraction reduction of 16.0% prior to shape optimization. After the shape optimization stage, the compliance deviations are substantially reduced for all cases, remaining within a narrow range of approximately 0.5%–2.6%, while the volume fraction deviations are fully eliminated in every example. The proposed framework effectively mitigates compliance deviation and consistently restores the target volume fraction, demonstrating its robustness and effectiveness across different structural configurations and geometric scales.

Figure 18: Mechanical performance deviations of reconstructed geometries relative to topology optimization results. (a) Cantilever beam; (b) Michell beam; (c) Half-MBB beam.

5.3 Engineering Example

A quadcopter frame is considered as a representative engineering example. Topology optimization is first performed to provide the basis for geometry reconstruction using both heuristic, visually guided methods and the proposed framework. The reconstructed designs are then comparatively evaluated through computer-aided engineering (CAE) simulations, and the superior design is further fabricated by 3D printing.

The initial structure of the frame is shown in Fig. 19a and is subsequently partitioned according to its functional requirements, as illustrated in Fig. 19b, where the rib region is defined as the design domain. The topology optimization results are presented in Fig. 19c, which serves as the basis for geometry reconstruction. It can be observed that the boundaries of the frame are rough, making the results difficult to be directly applied in practical engineering scenarios. Fig. 19d presents the reconstructed geometry obtained through heuristic methods.

Figure 19: Design workflow based on topology optimization. (a) Initial structure; (b) Design domain partitioning; (c) Topology optimization; (d) Reconstructed geometry via heuristic methods.

Using the reconstructed geometry based on heuristic methods as a reference, contour points are uniformly sampled along the straight segments of the contours and used as inputs for the proposed two-stage framework, as illustrated in Fig. 20a. Fig. 20b shows the optimized boundaries obtained under compliance minimization without increasing the structural volume. Fig. 20c illustrates the boundary control points that are imported into CAD software and used to sketch the optimized boundaries onto the initial structure. Fig. 20d then presents the reconstructed geometry generated by subsequent hole-cutting operations.

Figure 20: Reconstructed geometry via the two-stage framework. (a) Input: sampled contour points; (b) Output: optimized boundaries; (c) Sketching; (d) Hole cutting.

As shown in Fig. 21a, under linear elastic assumptions and identical loading conditions, the model reconstructed using the proposed framework exhibits reduced deformation compared with the heuristic reconstruction, with the maximum displacement decreasing from 10.448 to 9.342 mm. This reduction can be attributed to the fact that the shape optimization procedure targets compliance minimization, which redistributes material toward regions of high strain energy and better defines the primary load paths, thereby enhancing structural stiffness. In addition, the stress distribution is improved after the optimization. As shown in Fig. 21b, the optimized structure exhibits more regions with lower stress levels. The maximum stress decreases from 6488.2 to 4573.1 MPa, indicating that the proposed framework also leads to a noticeable improvement in stress performance by generating smoother structural boundaries.

Figure 21: Comparison of deformation and stress fields for reconstructed geometries. (a) Deformation results; (b) Stress field.

Subsequently, the optimized frame is fabricated using 3D printing, and the resulting physical prototype is shown in Fig. 22.

Figure 22: 3D printed frame.

Overall, the numerical examples and the engineering case demonstrate that the proposed framework can effectively reconstruct topology optimization results while preserving the structural performance. Moreover, the reconstructed geometry is fully compatible with standard CAD, CAE, and CAM workflows, thereby facilitating its integration into practical engineering design and industrial manufacturing processes.

6  Conclusions

This study presents an integrated reconstruction—optimization framework to address the challenges of translating implicit topology optimization results into practical engineering designs. The framework reconstructs explicit and editable geometric models from discrete density fields and further refines the structural boundaries through shape optimization to preserve structural performance. Its effectiveness and robustness are validated through representative numerical examples as well as a practical engineering case.

•   The proposed method is first validated using representative numerical benchmarks, where the compliance deviation is reduced to within 0.5%–2.6% while satisfying the prescribed volume fraction constraint.

•   The method is further applied to a practical engineering example of a quadcopter frame, where the reconstructed structure exhibits reduced maximum displacement compared with heuristic reconstruction approaches.

•   These results demonstrate that the proposed framework establishes a practical workflow capable of generating parameterized geometries that are directly compatible with CAD and CAE environments.

By enabling explicit geometry reconstruction while maintaining structural performance, the framework facilitates a seamless connection between topology optimization, engineering analysis, and manufacturing processes. This integration improves the efficiency of the design workflow by reducing the need for manual intervention and additional optimization efforts. While the proposed framework shows promising results, it is currently limited to two-dimensional problems. Extending it to three-dimensional cases remains nontrivial due to increased computational cost and the more complex evaluation of geometric interference. Future work will focus on addressing these challenges while maintaining the performance-consistent reconstruction capability.

Acknowledgement: The authors acknowledge the use of ChatGPT for language polishing of the manuscript. All technical content, analyses, and conclusions are the sole responsibility of the authors.

Funding Statement: This research is funded by the National Natural Science Foundation of China under Grant No. 92371206.

Author Contributions: Conceptualization, Yuting Tang and Yu Li; methodology, Yuting Tang, Yu Li and Xingyu Xiang; software, Yuting Tang and Jiaxiang Luo; validation, Yuting Tang; formal analysis, Yuting Tang and Yu Li; investigation, Yuting Tang and Yu Li; resources, Yuting Tang and Wen Yao; data curation, Yuting Tang; writing—original draft preparation, Yuting Tang; writing—review and editing, Yuting Tang and Yu Li; visualization, Yuting Tang and Xingyu Xiang; supervision, Yu Li; project administration, Yu Li, Weien Zhou and Wen Yao; funding acquisition, Wen Yao. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The authors confirm that the data supporting the findings of this study are available within the article.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

Abbreviations

Appendix A Parametric Sensitivity and Robustness Analysis of Structural Performance and Shape

The threshold value τ, smoothing window size Nw, and the filtering radius rf serve as key parameters of the reconstruction framework, yet their specific impacts on the final reconstruction remain to be fully characterized. Therefore, a parametric sensitivity analysis is conducted to systematically evaluate how variations in each parameter affect the mechanical performance and geometric shape.

Table A1 summarizes the compliance and volume fraction of the reconstructed geometries across varying parameter settings. These quantitative trends are further visualized in Fig. A1, while Fig. A2 illustrates the comparative shape evolution across all configurations, highlighting the geometric differences induced by parameter variations.

Figure A1: Sensitivity analysis of compliance and volume fraction with respect to reconstruction parameters. (a) Threshold value τ; (b) Smoothing window size Nw; (c) Filtering radius rf.

Figure A2: Sensitivity of shape with respect to reconstruction parameters. (a) τ=0.5, Nw=3, rf=3; (b) τ=0.6, Nw=2, rf=3; (c) τ=0.6, Nw=3, rf=2; (d) τ=0.6, Nw=3, rf=3; (e) τ=0.6, Nw=3, rf=3; (f) τ=0.6, Nw=3, rf=3; (g) τ=0.7, Nw=3, rf=3; (h) τ=0.6, Nw=4, rf=3; (i) τ=0.6, Nw=3, rf=4; (j) τ=0.8, Nw=3, rf=3; (k) τ=0.6, Nw=5, rf=3; (l) τ=0.6, Nw=3, rf=5.

As illustrated in Fig. A1a, the extraction threshold τ emerges as the dominant control parameter governing global structural performance. A monotonic increase in τ induces a progressive reduction in volume fraction, accompanied by a corresponding rise in compliance. This trend arises because higher thresholds exclude low-density regions, whereas lower thresholds retain more material, thereby enhancing structural stiffness and deformation resistance. In contrast, Fig. A1b,c indicates that the smoothing window size Nw and filtering radius rf exhibit negligible sensitivity regarding global metrics. Since these parameters operate exclusively on extracted contour points, they serve to refine local geometric features rather than alter the global shape. This observation is corroborated by Fig. A2, which demonstrates that while variations in Nw and rf merely modulate boundary morphology, geometries reconstructed under varying τ values preserve a high degree of self-similarity. These findings validate that τ is the critical lever for tuning the material-stiffness trade-off, while Nw and rf function as secondary tools for geometric post-processing.

The robustness of the proposed framework is further quantified by evaluating performance deviations across diverse parameter settings. As summarized in Table A1, the relative error in compliance for all optimized shapes remains strictly within 1.8%, while the volume fraction deviation is maintained at exactly 0% across all cases, signifying exact adherence to material constraints. A notable exception occurs in the case of τ=0.5, Nw=3, and rf=3. Here, the low extraction threshold causes the initial reconstruction (Stage one) to exceed the target volume. Consequently, the subsequent optimization prioritizes aggressive material removal to satisfy the strict volume constraint, inevitably leading to a increase in compliance. In summary, the framework exhibits remarkable robustness against parameter fluctuations. The combination of zero volume error and minimal compliance variation confirms that the proposed approach yields stable, reliable structural designs regardless of minor variations in τ, Nw, or rf.

Appendix B Computational Cost Analysis for the Benchmark Case

According to the workflow illustrated in Fig. 1, the computational cost of the proposed framework is evaluated using the benchmark case. The central processing unit (CPU) time consumption of major steps is measured to assess the computational efficiency. The detailed time costs are summarized in Table A2.

The workflow consists of a one-time preprocessing stage (Stage one) and an iterative optimization stage. Stage one mainly involves geometric operations for generating the initial NURBS representation and requires approximately 0.89 s, indicating that the reconstruction introduces negligible computational overhead.

In the optimization stage, a single finite element analysis takes about 0.094 s per iteration. The evaluation of interference values and density sensitivities using the central finite difference scheme, costing approximately 1.38 and 1.36 s, respectively. Consequently, the average computational cost per iteration is about 2.41 s.

These results indicate that the computational cost of the key steps remains manageable. As the optimization converges within a limited number of iterations, the overall computational time remains acceptable, demonstrating that the proposed framework maintains practical computational efficiency even when central finite differences are used for sensitivity evaluation. All numerical experiments are conducted on a desktop configured with an AMD Ryzen 7 9800X3D processor and 48 GB of RAM.

References

1. Manguri A, Hassan H, Saeed N, Jankowski R. Topology, size, and shape optimization in civil engineering structures: a review. Comput Model Eng Sci. 2025;142:933–71. doi:10.32604/cmes.2025.059249. [Google Scholar] [CrossRef]

2. Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech. 1999;69(9):635–54. doi:10.1007/s004190050248. [Google Scholar] [CrossRef]

3. Wang MY, Wang X, Guo D. A level set method for structural topology optimization. Comput Methods Appl Mech Eng. 2003;192(1–2):227–46. doi:10.1016/S0045-7825(02)00559-5. [Google Scholar] [CrossRef]

4. Nguyen SH, Nguyen TN, Nguyen-Thoi T. A finite element level-set method for stress-based topology optimization of plate structures. Comput Math Appl. 2022;115:26–40. doi:10.1016/j.camwa.2022.04.001. [Google Scholar] [CrossRef]

5. Zhang W, Yuan J, Zhang J, Guo X. A new topology optimization approach based on Moving Morphable Components (MMC) and the ersatz material model. Struct Multidiscip Optim. 2016;53(6):1243–60. doi:10.1007/s00158-015-1372-3. [Google Scholar] [CrossRef]

6. Zhang W, Chen J, Zhu X, Zhou J, Xue D, Lei X, et al. Explicit three dimensional topology optimization via Moving Morphable Void (MMV) approach. Comput Methods Appl Mech Eng. 2017;322:590–614. doi:10.1016/j.cma.2017.05.002. [Google Scholar] [CrossRef]

7. Subedi SC, Verma CS, Suresh K. A review of methods for the geometric post-processing of topology optimized models. J Comput Inf Sci Eng. 2020;20(6):060801. doi:10.1115/1.4047429. [Google Scholar] [CrossRef]

8. Van Dijk NP, Maute K, Langelaar M, Van Keulen F. Level-set methods for structural topology optimization: a review. Struct Multidiscip Optim. 2013;48(3):437–72. doi:10.1007/s00158-013-0912-y. [Google Scholar] [CrossRef]

9. Li Z, Xu H, Zhang S. A comprehensive review of explicit topology optimization based on Moving Morphable Components (MMC) method. Arch Comput Methods Eng. 2024;31(5):2507–36. doi:10.1007/s11831-023-10053-8. [Google Scholar] [CrossRef]

10. Du B, Yao W, Zhao Y, Chen X. A moving morphable voids approach for topology optimization with closed B-splines. J Mech Des. 2019;141(8):081401. doi:10.1115/1.4043369. [Google Scholar] [CrossRef]

11. Rozvany G. The SIMP method in topology optimization-theoretical background, advantages and new applications. In: Proceedings of the 8th Symposium on Multidisciplinary Analysis and Optimization; 2000 Sep 6–8; Long Beach, CA, USA. doi:10.2514/6.2000-4738. [Google Scholar] [CrossRef]

12. Hsu YL, Hsu MS, Chen CT. Interpreting results from topology optimization using density contours. Comput Struct. 2001;79(10):1049–58. doi:10.1016/S0045-7949(00)00194-2. [Google Scholar] [CrossRef]

13. Chu S, Xiao M, Gao L, Gui L, Li H. An effective structural boundary processing method based on support vector machine for discrete topology optimization. In: Proceedings of the 2016 IEEE 20th International Conference on Computer Supported Cooperative Work in Design (CSCWD); 2016 May 4–6. Nanchang, China. p. 10–5; doi:10.1109/CSCWD.2016.7565954. [Google Scholar] [CrossRef]

14. Tang PS, Chang KH. Integration of topology and shape optimization for design of structural components. Struct Multidiscip Optim. 2001;22(1):65–82. doi:10.1007/PL00013282. [Google Scholar] [CrossRef]

15. Wu X, Shen H, Zhang T, Li H, Liu B, Yan S. Research on a boundary regularization geometric reconstruction method of topology optimization structure based on the Freeman code. Adv Mech Eng. 2023;15(2):16878132231152954. doi:10.1177/16878132231152954. [Google Scholar] [CrossRef]

16. Li X, Qin C, Wei P, Su C. A boundary density evolutionary topology optimization of continuum structures with smooth boundaries. Int J Numer Methods Eng. 2022;123(1):158–79. doi:10.1002/nme.6851. [Google Scholar] [CrossRef]

17. Swierstra MK, Gupta DK, Langelaar M. Automated and accurate geometry extraction and shape optimization of 3D topology optimization results. arXiv:2004.05448. 2020. doi:10.48550/arXiv.2004.05448. [Google Scholar] [CrossRef]

18. Li Z, Lee TU, Yao Y, Xie YM. Smoothing topology optimization results using pre-built lookup tables. Adv Eng Softw. 2022;173:103204. doi:10.1016/j.advengsoft.2022.103204. [Google Scholar] [CrossRef]

19. Ježek O, Kopačka J, Isoz M, Gabriel D, Maršálek P, Šotola M, et al. Smooth geometry extraction from SIMP topology optimization: signed distance function approach with volume preservation. Adv Eng Softw. 2026;212:104071. doi:10.1016/j.advengsoft.2025.104071. [Google Scholar] [CrossRef]

20. Lin CY, Chao LS. Automated image interpretation for integrated topology and shape optimization. Struct Multidiscip Optim. 2000;20(2):125–37. doi:10.1007/s001580050144. [Google Scholar] [CrossRef]

21. Lin CY, Lin SH. Artificial neural network based hole image interpretation techniques for integrated topology and shape optimization. Comput Methods Appl Mech Eng. 2005;194(36–38):3817–37. doi:10.1016/j.cma.2004.09.005. [Google Scholar] [CrossRef]

22. Yildiz A, Öztürk N, Kaya N, Öztürk F. Integrated optimal topology design and shape optimization using neural networks. Struct Multidiscip Optim. 2003;25(4):251–60. doi:10.1007/s00158-003-0300-0. [Google Scholar] [CrossRef]

23. Gamache JF, Vadean A, Noirot-Nérin É, Beaini D, Achiche S. Image-based truss recognition for density-based topology optimization approach. Struct Multidiscip Optim. 2018;58(6):2697–709. doi:10.1007/s00158-018-2028-x. [Google Scholar] [CrossRef]

24. Joshi S, Medina J, Menhorn F, Reiz S, Rueth B, Wannerberg E, et al. Cad-integrated topology optimization. Munich, Germany: Technical University of Munich; 2016. [Google Scholar]

25. Liu S, Li Q, Liu J, Chen W, Zhang Y. A realization method for transforming a topology optimization design into additive manufacturing structures. Engineering. 2018;4(2):277–85. doi:10.1016/j.eng.2017.09.002. [Google Scholar] [CrossRef]

26. Chacón J, Bellido J, Donoso A. Integration of topology optimized designs into CAD/CAM via an IGES translator. Struct Multidiscip Optim. 2014;50(6):1115–25. doi:10.1007/s00158-014-1099-6. [Google Scholar] [CrossRef]

27. Koguchi A, Kikuchi N. A surface reconstruction algorithm for topology optimization. Eng Comput. 2006;22(1):1–10. doi:10.1007/s00366-006-0023-0. [Google Scholar] [CrossRef]

28. Yoely YM, Amir O, Hanniel I. Topology and shape optimization with explicit geometric constraints using a spline-based representation and a fixed grid. Procedia Manuf. 2018;21:189–96. doi:10.1016/j.promfg.2018.02.110. [Google Scholar] [CrossRef]

29. Hsu MH, Hsu YL. Interpreting three-dimensional structural topology optimization results. Comput Struct. 2005;83(4–5):327–37. doi:10.1016/j.compstruc.2004.09.005. [Google Scholar] [CrossRef]

30. Cuillière JC, Francois V, Drouet JM. Towards the integration of topology optimization into the cad process. Comput-Aided Des Appl. 2014;11(2):120–40. doi:10.1080/16864360.2014.846067. [Google Scholar] [CrossRef]

31. Wen Z, Zhai X, Kang H. From density to CAD: automatic and robust CAD model generation of topology optimization results via sparse optimization. Adv Eng Softw. 2025;209:103977. doi:10.1016/j.advengsoft.2025.103977. [Google Scholar] [CrossRef]

32. Larsen S, Jensen CG. Converting topology optimization results into parametric CAD models. Comput-Aided Des Appl. 2009;6(3):407–18. doi:10.3722/cadaps.2009.407-418. [Google Scholar] [CrossRef]

33. Bacciaglia A, Ceruti A, Liverani A. Surface smoothing for topological optimized 3D models. Struct Multidiscip Optim. 2021;64(6):3453–72. doi:10.1007/s00158-021-03027-6. [Google Scholar] [CrossRef]

34. Ren H, Xia B, Wang W, Chen X, Liu Y, Zhou P, et al. AMRTO: automated CAD model reconstruction of topology optimization result. Comput Methods Appl Mech Eng. 2025;435:117673. doi:10.1016/j.cma.2024.117673. [Google Scholar] [CrossRef]

35. Sigmund O. A 99 line topology optimization code written in Matlab. Struct Multidiscip Optim. 2001;21(2):120–7. doi:10.1007/s001580050176. [Google Scholar] [CrossRef]

36. Maple C. Geometric design and space planning using the marching squares and marching cube algorithms. In: Proceedings of the 2003 International Conference on Geometric Modeling and Graphics; 2003 Jul 16–18; London, UK. p. 90–5. doi:10.1109/GMAG.2003.1219671. [Google Scholar] [CrossRef]

37. Piegl L, Tiller W. The NURBS book. 2nd ed. Berlin/Heidelberg, Germany: Springer; 1997. doi:10.1007/978-3-642-59223-2. [Google Scholar] [CrossRef]

38. Mendelson B. Introduction to topology. Washington, DC, USA: Courier Corporation; 1990. [Google Scholar]

39. Bose P, Toussaint G. Geometric and computational aspects of manufacturing processes. Comput Graph. 1994;18(4):487–97. doi:10.1016/0097-8493(94)90061-2. [Google Scholar] [CrossRef]

40. Sutherland IE, Sproull RF, Schumacker RA. A characterization of ten hidden-surface algorithms. ACM Comput Surv. 1974;6(1):1–55. doi:10.1145/356625.356626. [Google Scholar] [CrossRef]

41. Yu D, Cai S, Gao T, Zhang W. A 168-line MATLAB code for topology optimization with the adaptive bubble method (ABM). Struct Multidiscip Optim. 2023;66(1):10. doi:10.1007/s00158-022-03403-w. [Google Scholar] [CrossRef]

42. Svanberg K. The method of moving asymptotes—a new method for structural optimization. Int J Numer Methods Eng. 1987;24(2):359–73. doi:10.1002/nme.1620240207. [Google Scholar] [CrossRef]