TreeTOp: Topology Optimization using Constructive Solid Geometry Trees
Rahul Kumar Padhy Pramod Thombre Krishnan Suresh Aaditya Chandrasekhar

Topological communities in complex networks
    Luís F Seoane
    
===========================================




§ ABSTRACT 

Feature-mapping methods for topology optimization (FMTO) facilitate direct geometry extraction by leveraging high-level geometric descriptions of the designs. However, FMTO often relies solely on Boolean unions, which can restrict the design space. This work proposes an FMTO framework leveraging an expanded set of Boolean operations, namely, union, intersection, and subtraction. The optimization process entails determining the primitives and the optimal Boolean operation tree. In particular, the framework leverages a recently proposed unified Boolean operation approach. This approach presents a continuous and differentiable function that interpolates the Boolean operations, enabling gradient-based optimization. The proposed methodology is agnostic to the specific primitive parametrization and is showcased through various numerical examples.










§ INTRODUCTION


Topology optimization (TO) has emerged as a powerful computational tool for achieving optimal material distribution within a design domain, subject to constraints <cit.>. The maturity of TO methods is evident in their widespread industrial application, particularly facilitated by commercial software.

While various TO methods exist, density-based approaches are widely adopted <cit.>. Popular density-based approaches discretize the design space into finite elements and optimize a fictitious material density within each element to generate organic, free-form designs <cit.>. Consider, for instance, the design domain and boundary conditions in <ref>(a). An optimal design that maximizes stiffness, subject to a volume constraint, using density-based TO is illustrated in <ref>(b). While offering design freedom, the resulting designs can be challenging to interpret <cit.> and modify <cit.>. Furthermore, the density-based TO designs often require extensive post-processing <cit.> that leads to a deviation between the final design and the optimal solution <cit.>. This deviation becomes more pronounced when the structure is manufactured using stock material, with components of fixed shape but variable dimensions <cit.>.


 
To address these limitations, alternative techniques, referred to as feature mapping-based topology optimization (FMTO) <cit.>, have emerged. These approaches utilize high-level geometric descriptors (or primitives), to represent the design. For the design domain and boundary conditions illustrated in <ref>(a), an optimized design obtained using FMTO is presented in <ref>(b). By parameterizing the design using primitives, FMTO facilitates the upfront enforcement of manufacturing rules, and also aids in design interpretation <cit.>. However, FMTO methods often rely solely on Boolean unions, which can limit design flexibility and fail to exploit the complex geometric operations available in modern computer-aided designing (CAD) systems <cit.>. While some FMTO methods incorporate Boolean operations other than union, the constructive solid geometry (CSG) tree is typically predefined <cit.>. The optimization process is then limited to reordering operations or splitting features to introduce new branches <cit.>.





 §.§ Contributions




This work introduces an FMTO framework with the following contributions: 


  * Expanded Boolean Operations: We extend traditional feature-mapping methods, typically limited to Boolean unions, to incorporate subtraction and intersection operations, significantly expanding the design space. Specifically, we employ a unified Boolean operation formulation <cit.>, where an interpolation function represents various Boolean operations.  This formulation enables continuous transitions between different Boolean operations, thus facilitating gradient-based optimization.


  * Concurrent Optimization: Unlike prior research <cit.> that focuses on optimizing only primitive parameters with fixed operators, our framework optimizes both primitive parameters and the Boolean operators simultaneously.


Since we rely on the unified Boolean operation presented in <cit.>, we also inherit the following features:


  * Mitigating Vanishing Gradients: The utilization of unified Boolean operations, facilitated by a continuous and monotonic interpolation scheme, enables gradient-based optimization while mitigating the issue of vanishing gradients.


  * Primitive Parameterization Flexibility: The proposed method is agnostic to the specific primitive parameterization, accommodating a wide range of geometric primitives.






§ RELATED WORK



As mentioned earlier, FMTO parameterizes the topology via high-level geometric features <cit.> or primitives that are mapped onto a mesh for analysis <cit.>. FMTO is driven by the need to integrate features into free-form designs, control specific structural dimensions, utilize stock materials, and generate a geometric representation that can be interpreted by CAD systems <cit.>. Some of the earliest work in this area combines free-form topology optimization with embedded geometric primitives <cit.>; for a detailed review, see <cit.>.

We focus on methods that represent the design exclusively using geometric primitives <cit.>. In particular, we discuss four popular approaches that represent the design as a union of primitives capable of translating, scaling, and modifying their shapes <cit.>:



  * Moving morphable components/voids (MMC/MMV) method: The MMC/MMV method uses geometric primitives, such as B-spline-shaped holes or components, to represent the design <cit.>. This approach allows for control of the design's geometry by explicitly manipulating the boundaries of these primitives <cit.>.


  * Geometry projection (GP) method: The GP method <cit.> uses geometric primitives, such as bars <cit.>, supershapes <cit.> and plates <cit.>, to optimize structural designs by projecting these primitives onto a fixed finite element mesh. This approach has been applied in various contexts, including 3D topology optimization <cit.>, multi-material design <cit.>, and the optimization of unit cells for lattice materials <cit.>. 


  * Method of Moving Morphable Bars (MMB): The MMB uses round-ended bars as primitives, allowing them to overlap and modify their shape and position within the design <cit.>.


  * Moving Node Approach (MNA): Finally, the MNA uses polynomial functions to project geometric primitives, representing the building blocks of the design as mass nodes <cit.>. 



The above approaches rely on the Boolean union of primitives <cit.>, which can restrict design flexibility.  We propose a generalized framework in (<ref>) to enhance design flexibility. <Ref> covers the framework's key components: primitive parametrization (<ref>), projection of primitives onto the density field (<ref>), the use of an expanded set of Boolean operations (union, intersection, subtraction) on the density fields (<ref>), and optimization strategy (<ref>). Several examples demonstrating the proposed method are presented in <ref>. Finally, <ref> concludes this work.





§ PROPOSED METHOD






 §.§ Overview


This study focuses on gradient-based topology optimization, which minimizes structural compliance under a volume constraint. We assume that the design domain, loads and restraints, has been prescribed. While the framework is agnostic to the type of primitives used, we will assume that the primitives are polygons for ease of implementation. The CSG tree is assumed to be a perfectly balanced binary tree with a specified depth. Our objective is to obtain an optimal configuration of the polygons at the leaf nodes and Boolean operations at all non-leaf nodes that, when applied, results in the optimized design; see <ref>.





 §.§ Primitive Parametrization



We populate the leaf nodes of a perfectly balanced binary CSG with polygonal primitives. In particular, a tree with a (given) depth of n_d (discussed in <ref>) would correspond to n_p = 2^n_d primitives at the leaf nodes. The primitives are parameterized using the method proposed in <cit.> where:


    
  * Each polygon primitive p^(i) is associated with a reference point (c_x^(i), c_y^(i)), and is the intersection of S (S ≥ 3) (straight-line) half-spaces, resulting in a polygon with {3, …, S} sides; S = 6 in <ref>(a). 

    
  * Each half-space h_j^(i) is initially oriented at an angle α̃^(i) = [0, π/3, 2π/3, …] (see <ref> (a)) and has an offset distance d_j^(i) > 0 from the reference point (see <ref> (b)).
    
  * To vary the overall orientation of the polygon, we allow for the rotation of all of its half-spaces by an angle θ^(i), resulting in the final orientation angle α_j^(i) = θ^(i) + α̃_j^(i) for the j^th half-space (see <ref>(c)).

In summary, each polygon p^(i) is parameterized as:


    p^(i) = {c_x^(i), c_y^(i), θ^(i), d_1^(i), …, d_S^(i)},     i = 1,…, n_p


Note that, by construction, each polygon is non-empty and can have between 3 and S sides <cit.>.




 §.§ Geometry Projection



The aim of geometry projection is to map primitives, defined by the polygon's parameters (<ref>) onto a density field defined over a mesh  <cit.>. To maintain differentiability, we first map each primitive (<ref>(a)) to a signed distance field (SDF), where the value of the SDF at any point (x, y) is defined as the shortest signed-distance to the primitive’s boundary (inside being negative and outside being positive). We begin by defining the SDF  of each half-space as (<ref>(b)):

 
    ϕ̂_j^(i)(x, y) = (x - c_x^(i))cos(α_j^(i)) + (y - c_y^(i))sin(α_j^(i)) - d_j^(i)



 
Since the SDF value of a point inside the polygon with respect to a half-space is negative, the SDF of the entire polygon at a given point is obtained by computing the maximum of the SDFs from each half-space at that point (<ref>(c)). To ensure differentiability, we use the LogSumExp (LSE) approximation of maximum <cit.>. This yields the SDF of the i^th polygon:


    ϕ^(i)(x, y) = l_0/tlog( ∑_j=1^Sexp(t/l_0ϕ̂_j^(i)(x, y)) )

where, t is the LSE scaling factor (discussed later in <ref>) and l_0 is the length of the diagonal of the domain bounding box.

Next, we compute the density field of the polygon ρ̃(·), from ϕ(·) using the Sigmoid function; see <ref>(d): 

    ρ̃^(i)(x,y) = 1/1 + exp(-γ/l_0ϕ^(i)(x,y))

 
Observe that negative SDF values (representing the polygon's interior) are mapped to a density of one (solid), while positive values (representing the polygon's exterior) are mapped to zero (void). SDF values near zero (representing the polygon's boundary), are projected to intermediate density values ρ̃∈ (0,1), with the sharpness of transition controlled by the hyperparameter γ.

Finally, we impose a threshold filter on ρ̃ to drive intermediate densities towards 1/0 as <cit.>:


    ρ(ρ̃) = tanh(β/2) + tanh(β (ρ̃ - 1/2))/2 tanh(β/2)

where the parameter β controls the sharpness of the threshold function. The obtained primitive densities are subsequently combined using the Boolean operations (as detailed in the following section) to obtain the design density.



 §.§ Unified Boolean Operations



The final design density is obtained by applying Boolean operators (<ref>) to the primitive density fields (operands)  defined in the previous section. A successful Boolean operations approach suitable for gradient-based optimization (<ref>) should possess the following properties:



  * Differentiability with respect to operands (density fields).

  * Differentiability with respect to the Boolean operators (union (∪), intersection(∩), difference(↦), negative difference(); see <ref>).


 

FMTO methods that demonstrate differentiability with respect to the operands <cit.> have been extensively employed. Further, the union and intersection operators are often approximated using smooth maximum/minimum functions <cit.>.  Alternative techniques, such as the smooth blending operator  <cit.> and R-functions <cit.>, have been proposed, exhibiting differentiability with respect to the operands and individual operators.  However, a smooth transition between the operators is difficult to achieve.  Furthermore, while a unified Boolean approach based on the modified R-function <cit.> was introduced to facilitate differentiable transitions between operators, it failed to satisfy key Boolean operator axioms, resulting in unexpected behavior <cit.>.

To address these issues, we adopt a unified Boolean operation approach proposed in <cit.>. Let 𝒳 and 𝒴 represent two primitives (polygons), with their corresponding density functions ρ_x and ρ_y. Then, a generalized Boolean operation between the two primitives is defined as (<ref>): 


    ℬ(ρ_x, ρ_y ; b) = (b_1 + b_2) ρ_x + (b_1 + b_3) ρ_y + (b_0 - b_1 - b_2 - b_3) ρ_x ρ_y

where b = {b_0, b_1, b_2, b_3}  are interpolating parameters. The interpolating parameters are constrained by 0 ≤ b_i ≤ 1 and ∑_i=0^3 b_i = 1. When b is a one-hot vector, we recover the standard Boolean operators (<ref>): 




Observe that <ref> continuously interpolates between individual operators. For the two polygons in <ref>, the continuous interpolation between the intersection and union operators is illustrated in <ref>, while a generic operator state is illustrated in <ref>. 

 
 
 

Furthermore,  the function ℬ is differentiable with respect to both the operands (ρ_x, ρ_y) and the operators (b_i). For example, the derivative w.r.t operand ρ_x is:


    ∂ℬ/∂ρ_x = (b_1 + b_2) + (b_0 - b_1 - b_2 - b_3)ρ_y

while the derivative with respect to operator coefficient b_1 is

    ∂ℬ/∂ b_1 = ρ_x + ρ_y - ρ_xρ_y

 
Finally, with the interpolated operators defined at all non-leaf nodes and the primitive densities at the leaf nodes, we can evaluate the CSG tree bottom-up to obtain the  density field ρ(x,y) of the evolving design at any instance.



 §.§ Finite element analysis



For finite element analysis, we use here a structured mesh with bilinear quad elements. Having obtained the design density field ρ(x,y), let ρ_e be the density at the center of element e.  The corresponding Young's modulus E_e is obtained using the Solid Isotropic Material Penalization (SIMP) model <cit.> as:


    E_e = E_min + (E_0 - E_min)(ρ_e)^p

where E_0 is Young's modulus of the solid material, E_min is a small constant added to prevent a singular global stiffness matrix, and p is the SIMP penalty. One can evaluate the element stiffness matrix as

    [K_e] = E_e ∫_Ω_e [B]^T[D_0][B] d Ω_e


where [B] is the strain matrix, and [D_0] is the constitutive matrix with a Young’s modulus of unity and an assumed Poisson's ratio (see Table <ref>). This is followed by the assembly of the global stiffness matrix K. Finally, with the imposed nodal force vector f, we solve the state equation for the nodal displacements u = K^-1f. 



 §.§ Optimization



We now summarize various aspects of the optimization problem:



  
Objective: We consider here a compliance minimization problem. where the compliance is computed as J = f^T u.



  
Volume constraint: The design is subjected to a total volume constraint. With v_f^* being the maximum allowed volume fraction and v_e being the element areas, the volume constraint g_v is defined as:


    g_v ≡∑_e=1^N_eρ(x_e) v_e/v_f^* ∑_e=1^N_e v_e - 1 ≤ 0




  
Optimization variables:  The optimization variables are as follows:

    
  * Recall that the design is defined by polygon parameters, including center coordinates, angular offset, and plane distances, resulting in n_p(S+3) free parameters (see Section <ref>). Additionally, the design requires n_b = 2^n_d - 1 boolean operations (b⃗ = b^(1), …, b^(n_b)) for a tree of depth n_d. 
    
  * For optimization, we define an augmented normalized design vector z = [z_c_x, z_c_y, z_θ, z_d, z_b] that lie in [0,1]. With c_x, c_y, θ, d being the lower bound and c_x, c_y, θ, d being the upper bound on the parameters of the polygons, we can retrieve the unnormalized x-center as c_x ←c_x + (c_x - c_x) z_c_x and so on. The design variable z_b^(i), corresponding to the i^th boolean operator is transformed into one-hot encoding using a softmax function <cit.>.
    
  * All design variables z are uniform-randomly initialized (using ) with a default seed value of 2. 






  
TO problem: The final TO problem is posed as:  


	
	
    zminimize      J = f^T u(z) 
       subject to      K(z)u = f
              g_v (z) ≤ 0 
               0 ≤ z_i ≤ 1   ∀ i





  
Optimization method:  We employ the method of moving asymptotes (MMA) <cit.> to perform the design updates. Specifically, we use a Python implementation with all default parameters corresponding to the version of MMA presented in <cit.>. The choices of the MMA move limit, the step tolerance, the Karush-Kuhn-Tucker (KKT) tolerances, etc are specified later under numerical experiments.




 §.§ Sensitivity Analysis



Sensitivity computation, a crucial aspect of gradient-based optimization, involves determining the derivatives of the objective function and constraints with respect to the optimization parameters. Traditionally, this is conducted manually, which can be labor-intensive and error-prone. However, by utilizing the automatic differentiation (AD) <cit.> capabilities of frameworks such as JAX <cit.>, this step can be fully automated, ensuring accurate and efficient computation. In practical terms, we only need to define the forward expressions, and the derivatives of the objective and volume constraint with respect to the optimization variables is automatically computed with machine precision. Finally, we summarize the  algorithm for the proposed framework in <ref>.






§ NUMERICAL EXPERIMENTS



In this section, we conduct several experiments to illustrate the proposed framework. The default parameters in the experiments are summarized in Table <ref>. All experiments are conducted on a MacBook M2 Air, using the JAX library <cit.> in Python.  Additionally, all design parameters z are uniform randomly initialized.






 §.§ Validation




We compare the designs from our proposed framework with those obtained using the SIMP-based optimization <cit.>. Using the default parameters outlined in <ref>, and with n_d=4, we optimize both the Boolean operators and polygon parameters for an MBB beam (<ref>(a).  In the SIMP implementation, the density field is optimized with a filtering radius of 1.3, targeting a volume fraction of v_f^* = 0.5 ( <ref>(b)). For the proposed method, we display the final design in <ref>(c).  We observe that both methods yield similar designs and performance.  For illustration, we present the complete resulting CSG tree in <ref>.

Additionally, the convergence is illustrated in <ref>. The resulting density fields at the root nodes for the 10^th, 50^th, 75^th, and final (111^th) iterations are shown as insets. In comparison, the SIMP-based implementation took 126 iterations to converge. In our framework, the percentage of computational time is as follows: 0.6% for geometry projection, 8.1% for CSG tree evaluation, and 91.3% for FEA and sensitivity analysis; each iteration takes 1.6 seconds.



 







 §.§ Effect of Tree Depth


In this experiment, we examine the effect of the CSG tree depth n_d on the performance of the design. We revisit the MBB beam in <ref>(a), with v_f^* = 0.5, and investigate the impact of different values of n_d. The CSG tree for  n_d =2 and n_d =3 are illustrated in <ref>, while the final designs for n_d =4,5,6,7  are illustrated in <ref>. The performance improves as n_d increases; no significant improvements were observed beyond the depth of n_d = 4. Based on similar experiments, we recommend using a depth n_d ≥ 4; we use a depth of n_d = 6 for the remainder of the experiments.

 









 §.§ Flexibility of Framework



In this experiment, we demonstrate the framework's ability to obtain optimal designs with specific Boolean structures. We again consider the MBB problem (<ref>(a)) with v_f^* = 0.5. First, we set the root node operator as a difference operator, resulting in the design shown in <ref>(a). Observe that the optimization can sometimes result in empty nodes. These nodes and their descendants are detected and pruned  <cit.>. Next, we only allowed union operators; the resulting topology is illustrated in <ref>(b). 






 §.§ Mesh dependency



Next, we study the effect of the FE mesh on the computed topology using the MBB problem with v_f^*=0.5 (<ref>(a)), while keeping all other parameters at their default values. <Ref> presents the topologies for varying mesh sizes: 60 × 30, 80 × 40, and 100 × 50 elements. No significant difference in performance was observed across these mesh sizes. However, we note that the designs exhibit less noise as the number of mesh elements increases. The boundary exhibits small features for a coarse mesh, since the latter fails to capture the impact of small features on the performance. However, we have observed that these undesirable features tend to reduce as the mesh size increases.






 §.§ Effect of Initialization




Recall all design variables z are uniform-randomly initialized with a default seed value of 2. We now investigate the influence of the seed value on the optimal designs. While all topology optimization techniques are inherently sensitive to the initial design <cit.>, it has been observed that this dependency is particularly pronounced in feature-mapping techniques <cit.>. In this specific example, we consider the design of an MBB beam (<ref>(a)) with 100 × 50 mesh elements.  The resulting designs and their corresponding performances for various initial seeds are compared in <ref>. Observe that while we obtain diverse designs, the performances are similar. This suggests that, as expected, the landscape is highly non-convex with numerous local solutions.







 §.§ Pareto



An essential consideration during the design phase is the exploration of the Pareto front; evaluating the trade-offs associated with various design choices. Let us consider the mid cantilever beam problem illustrated in <ref>(a). With 100 × 50 mesh elements,  n_d = 6 and S = 6, we investigate the trade-off between the structure's compliance and volume fraction. <ref>(b) illustrates the resulting designs (and the intermediate densities at the first depth) at various volume fractions. As anticipated, we observe an increase in compliance as the allowable volume fraction is decreased. We did not impose explicit symmetry requirements on the CSG tree.


 




§ CONCLUSION



We introduce an FMTO framework that concurrently optimizes both the primitive parameters and the Boolean operators (union, intersection, and subtraction). Notably, by leveraging a unified Boolean operation approach, the framrwork improves upon existing FMTO methods that typically rely on predefined CSG trees and/or are limited to optimizing only primitive parameters operated with unions. The efficacy of the proposed method was demonstrated through various numerical examples.

The generated tree structure appears unnatural and dissimilar to human-designed trees. This presents an opportunity to combine the proposed method with machine learning-based frameworks <cit.> to mimic human-designed CSG trees. The method also suffers from some of the shortcomings inherent to FMTO methods. For example, the optimizer is highly susceptible to being trapped in local optima; it is also sensitive to initialization. Additionally, the optimizer may fail to converge to an optimal solution when fewer primitives are used <cit.>.

Nevertheless, several avenues for future research remain. For instance, manufacturing constraints such as symmetry, feature size control must be incorporated.  One promising direction is to incorporate CNC operations into the TO process, ensuring that designs are both optimal and manufacturable <cit.>. Further, optimizing the primitive type and depth of the tree in conjunction with the operators and primitive parameters requires investigation <cit.>.  




§ ACKNOWLEDGMENTS

The University of Wisconsin, Madison Graduate School supported this work. 



§ COMPLIANCE WITH ETHICAL STANDARDS

The authors declare that they have no conflict of interest.



§ REPLICATION OF RESULTS

The Python code is available at https://github.com/UW-ERSL/TreeTOpgithub.com//UW-ERSL/TreeTOp


unsrt  
  


