Creating shapes optimized for a particular function is one of the main tasks in computer-aided design. Modern fabrication methods have greatly simplified the creation of complex free-form shapes at an affordable cost and opened up new possibilities for improving functionality and customization through automatic optimization.
Methods for constructing optimal shapes have enjoyed considerable attention in a variety of settings, such as large-scale architectural forms, engine parts, footwear, medical prosthetic devices and metamaterial structures. Notably, most of the work focuses on designing continuous structures, fabricated from one material, or a set of materials fused or glued together.
In this work, our focus is on supporting shape optimization for assemblies, and more specifically for assemblies that are held together by contact and friction, a setting which has received relatively little attention. At the same time, assemblies are ubiquitous, as most manufactured objects around us are assembled from separate parts, often made from different materials. For example, the steel legs of a table may be inserted into openings of a wooden or MDF top; a phone may have a snap-on plastic protector or cover; a prosthetic device is attached to the body with friction. While the specific mechanisms for holding objects together may vary broadly, they are all based on combining deformation with contact and, in many cases, friction.
In all these cases, contact plays a major role in the function and mechanical behavior of the assembled object. Shape optimization helps to achieve better performance or save on the costs of material for fabrication (for additive fabrication these two are closely related). Some of the important measures of performance of assemblies include total volume or weight, total deformation energy, maximal stress and more complex measures such as permeability.
In this paper, we present a formulation and a robust numerical method for computing optimized shapes in the presence of contact and friction. Compared to shape optimization tasks not involving contact, the problem is significantly more difficult to solve, as it involves complex inequality constraints required for handling collisions and friction. The resulting problem is non-smooth, and often hard to solve sufficiently accurately. Our overall approach is to use a smoothed version of the problem [Eck et al. 2005] amenable to standard optimization techniques on the one hand, and allowing us to approximate the desired solution as close as possible on the other hand.
Contributions. In summary, the contributions include
• A shape derivative-based formulation for optimization problems with contact and friction, building on [Maury et al. 2017]; • A novel shape-optimization framework based on FEM discretization of this formulation, capable of handling contact regions between two deformable objects as well as a deformable and a rigid object. It provides sufficiently accurate elastic deformation computations to support, e.g., max stress reduction. • The frameworks supports conventional functionals (stressbased and volume) and new, contact-specific functionals (assembly and disassembly, parallel alignment) that ensure that connection strength is maintained and that at the same time the parts held together by contact can be assembled. The framework also supports optimization involving multiple load scenarios. • We demonstrate a range of 2D and 3D examples of shape optimization, and qualitatively evaluate these examples using laser cutting and 3D printing to fabricate them and demonstrate the expected behavior.
There is a broad range of work on shape/topology optimization and related methods, but relatively few works were trying to solve problems with contact. We focus on related work on shape optimization with contact, and briefly mention other shape/topology optimization research that we rely on.
Shape optimization with contact. Some previous works have considered contact of a soft body with rigid surfaces, for example, [Beremlijski et al. 2014], [Haslinger et al. 1986] and [Herskovits et al. 2000]. While other some other works have studied the interaction of two or more bodies in contact, like recent works from [Maury et al. 2017], [Desmorat 2007] and [Stupkiewicz et al. 2010]. Most papers do not consider friction, and those which do often consider simplified (compared to the standard Coulomb) friction models as discussed in [Maury et al. 2017].
For contact models, there are two families of algorithms, Lagrangian and the penalization methods. The first type of methods adds Lagrange terms for the model constraints to the objective function and uses sub-gradient-based optimization to deal with the fact that the problem is non-smooth. Examples can be found in [Herskovits et al. 2000] and [Stupkiewicz et al. 2010] , where a bilevel approach and an augmented Lagrangian method are used. The second type of methods, penalization methods, use smooth approximations to the problem, which add terms to the variational formulation. Our method belongs to this category. For this type of method, the objective function is smooth, which considerably simplifies optimization. In our work, we choose to use penalization method following the mathematical model presented by Eck et al. [Eck et al. 2005]. In [Maury et al. 2017], a similar overall approach with level-set discretization. While level-set modeling has many advantages (e.g., allowing for easy topology changes) in our experience, it is not well-suited for handling important types of problems, in particular, those involving stress reduction, which is our focus in this work, as stress is harder to resolve precisely.
Stress minimization. Maximal stress minimization was considered in a number of papers; the formulation closest to ours is [Panetta et al. 2017], where worst-case optimization for periodic metamaterial structures is considered; it uses parametric periodic shapes, which are meshed for shape derivative computation, and does not consider contact. Earlier work on minimization of maximal stress is [Allaire et al. 2004], which applies topology optimization to design lightweight minimal-stress objects built from sequentially laminated composites. Another similar work is [Allaire and Jouve 2008], which applies the level-set topology optimization method to minimize the 𝑝-norm of stress. None of these works consider contact, and level-set methods (with Eulerian discretizations on a fixed grid) require impractically fine meshes for accurate optimization of high 𝑝 norms of stress. Other works considering max stress include [Lian et al. 2017;Polajnar et al. 2017;Sonmez 2009;Van Miegroet and Duysinx 2007;Xia et al. 2012]. In the computer graphics community, [Stava et al. 2012] was one of the first works to introduce heuristic shape correction techniques that effectively result in stress reduction. [Zhao et al. 2016;Zhou et al. 2016] consider problems involving bounding von Misses tress.
Contact and friction modeling. The literature on contact is extensive, and we mention only few most closely related works here; for general theory see e.g., [Stewart 2001]. Typically, contact problems are viewed as constrained optimization problems, with per-elementpair constraints. In particular, contacts between deformable objects that we consider in this paper, are defined as constraints between surface primitives (triangles, edges and vertices); some examples include [Belytschko et al. 2000;Bridson et al. 2002;Harmon et al. 2009;Otaduy et al. 2009;Verschoor and Jalba 2019]. Penalty-based methods for handling these constraints are among the oldest, but were largely supplanted by constraint formulations and LCP (linear complemenarity)-based or SQP solution methods. [Harmon et al. 2009] developed a method for which progressively high penalties are applied as the distance decreases, growing arbitrarily large as the distances to the object decreases. A recent work [Geilinger et al. 2020] provides a differentiable method for solving dynamic problems with contact and Coulomb friction using a penalty based model combined with equality constraints for static friction, which can be used with gradient-based optimization. Also dealing with Coulomb Friction, [Ding and Schroeder 2020] proposes a penalty based solution that can be used coupling rigid to deformable bodies and material point method (MPM). Another recent method [Li et al. 2020] shares some aspects with the approach we use, specifically, smoothed version of constraints is used, in particular for friction. However, we use finite penalties, rather than infinite barriers for constraints, and add smoothing both to contact and friction formulations.
It is well-known (see, e.g., [Moreau 1973]) that friction introduces significant non-smoothness to solutions, and due to their dissipative nature, solutions of problems with friction cannot be obtained using energy minimization. A variety of iterative methods for solving resulting problem were developed [Alart and Curnier 1991;Daviet et al. 2011;Jean and Moreau 1992]. Nonsmooth solution methods were applied in recent work with some success [Macklin et al. 2019], but are very difficult to use in the optimization context. Smoothed version of friction were proposed both in simulation and optimization context in the works we have mentioned.
Contact-based Assemblies. It is also important to mention works designing contact-based assemblies using geometric techniques and not necessarily relying on elasticity simulations to achieve their goal. Works like [Panozzo et al. 2013;Vouga et al. 2012;Wang et al. 2019] consider stability for self-supporting surfaces. Other works [Sun and Zheng 2015;Ureta et al. 2016;Yao et al. 2017] consider the design of joints for more general objects, including furniture.
We start with a high-level overview of the general problem we are solving. The input to our algorithm is a collection of 3D meshed objects, some of which may be in contact (the meshes will be updated in the process of optimization).
We solve the problem of the general form
The main components of this formulation include the following.
• The unknowns in the optimization are shape parameters 𝑝, in our case displacements of mesh vertices on the boundary, defining the domain Ω(𝑝), and displacements 𝑢 of the points of Ω(𝑝) resulting from elastic deformation with contact. (We use 𝑢 for the continuous solution of the elasticity problem, and ū for the vector of displacements of vertices of a meshed Ω(𝑝)). Note that while 𝑢 corresponds to the deformations of the shape Ω(𝑝) in our simulation, displacements 𝑝 (on 𝜕Ω(0)) define how the rest-state object shape is changed by optimization. • The PDE constraint 𝐹 ( ū) = 0 is a FEM-discretized elasticity equation for ū on Ω(𝑝)with contact and friction, leading to a non-linear system of equations in ū. • The objectives in our optimization are of the one of two forms below, as we have objectives with computation on the whole object domain (e.g., stress) or only on its boundary (e.g., disassembly objective)
where both the domain of integration and the integrands depend on the unknowns, and 𝑒 represents a pointwise measure (e.g., stress), computed as a function of FEM solution defined by the discrete vector ū. Each objective can be either an optimization target to minimize as in (1) or a part of an inequality constraint (upper bound on an objective) described below.
As optimization targets, we consider 𝐿 𝑝 -norms of stress, yielding compliance for 𝑝 = 2 and a close approximation of the maximum stress for large 𝑝 and volume, of a part or the whole assembly. Other objectives are used primarily as constraints. • The inequality constraints, critical for our formulation, are of the form 𝐵( ū) ⩽ 𝐵 0 , where 𝐵( ū) is an integral of the type (2). 𝐿 𝑝 norm of stress (effectively allowing to bound maximal stress) and volume objectives can also be used as constraints, e.g., we can impose a bound on the maximum volume allowed.
Other objectives used primarily in constraints are the assembly constraint, ensuring that two parts can be assembled together; disassembly constraint, ensuring that once assembled, the object does not fall apart; and parallel alignment constraint, which makes it more difficult to disassemble the object using any direction other than the disassembly one.
The key element of the optimization process is computing objective and constraint gradients with respect to ū, required for any first-or second-order optimization algorithm.
There are two main approaches to this problem: one can discretize the problem first, fixing the mesh for Ω(𝑝) for all 𝑝, and a FEM basis. This converts the problem into a finite-dimensional nonlinear algebraic problem, and then compute gradients of the objectives and constraints with respect to 𝑝, the positions of boundary mesh vertices. The alternative is the "differentiate-first" approach. Specifically, for each objective and constraint 𝐽 (𝑢), before it is discretized, we derive its shape derivative, a continuous analog of the gradient with respect to ū. The shape derivative is a functional 𝑑 𝐽 [𝑣], where 𝑣 is a velocity field of the deformations of the domain Ω(𝑝) (i.e., the displaced position of a point 𝑞 is 𝑞 + 𝑡𝑣). 𝑑 𝐽 [𝑣] yields the rate of change of 𝐽 as 𝑡 → 0. 𝑑 𝐽 is defined by a function 𝜌 on Ω, which can be obtained by solving a 𝑃𝐷𝐸, similar in structure to the elasticity equation (the adjoint equation). The advantage of the latter approach is that it naturally allows for remeshing and refinement (the adjoint equation solver is just a standard elasticity solver), which is essential for large changes in the object shape: any method used to solve large-deformation elasticity PDEs can be used to compute 𝜌, without fixing the discretization in advance.
The cost of each optimization step is approximately equal to solving two elasticity problems, the nonlinear elasticity problem, and the linear adjoint problem, to obtain the solution and shape derivative respectively.
In the rest of the section, we provide more specifics on the components of the formulation described above and their discretizations. The complete derivations are included in the supplementary material.
We use the following notation, also illustrated in Figure 2 • Ω: optimization domain (may consist of multiple objects) Ω Ω ΓD ΓN ΓC ΓS Fig. 2. Notation for boundary parts for a simple assembly: Γ 𝐷 is attached to the vertical wall (red), Γ 𝐶 can slide along the ground (yellow), Γ 𝑁 is where the load is applied, and Γ 𝑆 is the contact area of two parts. We view both parts as a single domain Ω consisting of detached parts.
In this paper, we only consider static problems, so there is no time dependence in the equations we use.
The basic form of static equations of linearized elasticity (without contact or friction yet) is
The variational form of this problem, that we rely on to formulate the contact constraints in a computationally practical way, as well as for the finite element discretization, is given by the equation
satisfied for any 𝑤 in an appropriate function space with 𝑤 = 0 on Γ 𝐷 , and 𝑢 satisfying the Dirichlet boundary conditions on Γ 𝐷 .
Contact constraints. We consider two forms of contact constraints, rigid-deformable (RD) and deformable-deformable (DD), defined on the parts of the boundary Γ 𝐶 and Γ 𝑆 respectively.
The former type of constraints involves a fixed boundary to which we refer as obstacle, on one side, which can be considered rigid, e.g., legs in contact with the floor.
The first equation on Γ 𝐶 says that the displacement should move the object points away from the obstacle; the second equation says that the normal force on this boundary should be 0 or point towards the other object; and finally, the last equation, the complementarity condition, ensures that if the force is nonzero, the displacement in the normal direction is zero, i.e., there is contact. The DD contact constraints are similar, however, because we have deformable material on both sides of the boundary, and the differences in stress need to be considered. Normals 𝑛 -and 𝑛 + and displacements 𝑢 - and 𝑢 + correspond to the parts on the opposite sides of contact, and
Friction constraints. Similarly, there is a set of equations for friction:
where 𝜆 > 0, and 𝜇 is the friction coefficient; 𝑡 and 𝑛 refer to tangential and normal components of the force. The first inequality captures the main aspect of Coulomb static friction model (the force is bounded by 𝜇 times the normal force). The second equation states that if the force is below maximal, no displacement happens, and the third one that the displacement at maximal friction force is parallel to it and opposite in direction.
A modification (6) (replacing one-sided constraints with differences), applies to friction (please see the supplementary document), yielding deformable-deformable contact.
While the basic elasticity problem (4) is quite straightforward to solve and well-understood, contact constraints, especially friction constraints, result in numerous difficulties: (1) the system becomes highly nonlinear and non-convex; (2) the solutions may be nonsmooth; (3) due to dissipative nature of friction, the problem cannot be cast as a (constrained) energy minimization problem. This makes even the direct solution of the problem difficult to make robust, and presents a particular challenge for shape optimization.
The key mathematical ideas for resolving these difficulties can be found in [Eck et al. 2005], further developed in [Maury et al. 2017]. (Recent work [Li et al. 2020] also follows a related approach for dynamic deformable contact).
The main elements of the approach include:
• Use smooth approximations of contact with friction (these are used to prove solution existence in [Eck et al. 2005], but are also valuable computationally). In particular, the norms in ( 7), that result in non-smoothness of the constraints are replaced by a smoothed norm that will be described in detail in this section. • Replace the inequality constraints by suitable (also smooth) penalty functions, added to the variational formulation of the problem (4); a number of works (see [Maury et al. 2017] for discussion) show that the results converge to the true solution, as the penalty weight 1/𝛼 → ∞. Like penalty approaches commonly used in graphics, using these methods has the benefit of replacing constrained optimization with unconstrained, and the downside of potential constraint violations which we discuss below. We largely follow the smooth formulation of [Eck et al. 2005], which demonstrates that the solution to the smooth problem exists under typical assumptions, and with smoothing parameters approaching zero, converges to a solution of the original problem.
Smoothed contact functional. We introduce a smoothing function [Eck et al. 2005]
approximating 𝑚𝑎𝑥 (0, 𝑦) as 𝜂 → 0. This is the same function used in [Maury et al. 2017].
The contact equations (5) lead to the following objective
As 𝛼 → 0, the solutions obtained with this term added to (4) converge to solutions satisfying the constraint 𝑢 • 𝑛 ⩽ 0. Similarly, we have
for DD contact.
Friction functionals. The Coulomb friction constraints ( 7) and the analogous constraint for DD friction, are expressed in terms of norms. For the smoothed functional, we use the following function to approximate |𝑦|, as 𝜂 𝑛 → 0:
Then friction constraints are captured by the following objectives [Eck et al. 2005]:
Smoothed elasticity with contact and friction. Four objectives ( 9)-( 12) are added to (3) to model contact with friction both for deformabledeformable and deformable-rigid contacts on parts of boundaries.
for any 𝑤. The finite element discretization of these equations is standard and leads to a non-linear system of equations, because of the nonlinear functions. This system requires computing a Jacobian of the left-hand-side to solve efficiently. We refer to the supplementary document for the derivation of the Jacobian.
Next, we describe the optimization problem we solve, specifically, the set of objectives we use in the functional 𝐽 as optimization targets or inequality constraints, and how their shape derivatives are computed.
In our model, we consider only the simplest type of assembly, namely, moving parts together in a specific direction 𝑦. This can be generalized to nonlinear trajectory settings (e.g., a screw motion trajectory) but this will make these constraints substantially more complex.
𝐿 𝑝 stress and volume. These objectives are standard, and expressed as
The stress objective for large 𝑝 approximates the non-smooth maximum-stress functional well. Following [Panetta et al. 2017], we consider the average stress at each element. Also, we use value 𝑝 = 20 in most cases. We note that von Mises stress can be used just as easily. In addition, we have options for setting target stress and target volume, which are expressed as follows:
where 𝑆 𝑡 and 𝑉 𝑡 represent the stress and volume targets and 𝜑 is the following function
We normalize our objective approximating maximal stress by the ratio between the load (in our boundary conditions) and the area of our shape's surface, which has units of stress. For the volume objective, we use 1/𝑉 2 𝑡 as a normalization constant. Fig. 4. Assembly, disassembly and parallel alignment objectives. Top: left shape is easy to assemble in -𝑦 direction, while shape on the right cannot be assembled without effort (high objective value, due to normals pointing in the direction opposite to 𝑦). Middle: For the same disassembly direction 𝑦, and user-defined load 𝑙, the assembly on the right falls apart under the load; for the shape on the left, the parts cannot move apart in direction 𝑦 under load 𝑙. Bottom: left shape has a large region parallel to disassembly direction (marked in blue), making it hard to disassemble in any direction other than 𝑦, while shape on the right allows for multiple easy disassembly directions (shown in light red).
We assume an disassembly direction 𝑦 per connection, with assembly happening in opposite direction -𝑦. We define two objectives, an assembly objective, which ensures that parts can be put together without much deformation in direction -𝑦, and dissasembly objective, which penalizes parts moving apart along 𝑦, under the loads specified by the user.
Assembly objective. A fully physical treatment of assembly would require simulating the assembly process, leading to a large number of nonlinear solves for different time steps for a single gradient evaluation. Instead, we use a geometric heuristic with negligible cost to allow parts in contact to be assembled without a large deformation while moving a part along an a priori fixed direction in space. Assembly is ensured if the contact surface is a height field when we view 𝑦 as a vertical direction. The angle between the direction of connection/disconnection and surface normal needs to be less than 𝜋/2. In our experiments, we always start with assemblable non-optimized parts, in order to facilitate the optimization. This is not a hard constraint for deformable objects, as it can be violated by a small amount, to squeeze in a connector. We do not attempt to bound maximal possible deformation, but empirically it can be easily controlled by increasing the weight of this objective.
where 𝑦 is the disassembly direction vector (and -𝑦 the assembly direction), 𝑛(𝑥) is the normal at the contact position 𝑥, and
Disassembly objective. This objective ensures that under userdefined loads, or for a collection of different user-defined loads, the parts should not move away from the assembled positions. We formalize this by requiring that the optimized shape in equilibrium does not move in the disassembly direction, i.e., the forces holding contact points together are sufficient.
A tolerance 𝑢 𝑡𝑜𝑙 is added, allowing a small displacement in the assembly direction. While the assembly term is purely geometric, the disassembly objective relies on actual simulation responses and it bounds how much optimized parts can move. As long as it is finite, the assembly does not fall apart.
Although the last two objective terms may work against each other in the cases when the loads have a significant component along the disassembly direction, in other situations they work together to provide structures that are, at the same time, both assemblable and are not prone to accidental disassembly under user-defined load. See Figure 33 for an example.
Parallel alignment objective. Finally, this objective is a secondary heuristic that is not essential for solving the problems but empirically makes the structure more robust with respect to being disassembled by forces close to assembly direction, since it favors solutions with larger contact zones tangential to the assembly direction. When such zones are present, forces deviating from the disassembly direction have normal components resulting in deformations and friction keeping the parts from accidental disassembly. In the bottom of Figure 4, we can see how the shape without parallel alignment allows for easy disassembly in different directions (shown in red), without any reaction force (normal or friction), while the shape on the left allows for a single disassembly direction 𝑦.
This term is defined per continuous contact area Γ 𝑆 with respect to a disassembly direction 𝑦:
where 𝑊 target is the target percentage area in the connection which should be parallel to the disassembly direction, which can help increasing/decreasing reaction forces avoiding disassembly. The function 𝑔 is parameterized by 𝜖, defining the interval for which 𝑔 is positive and 𝛽 setting the linear decay of the function when 𝑎𝑏𝑠 (𝑧) is larger than 𝜖. It has the following form:
where 𝜖 is a smoothing constant and where 𝑎 and 𝑏 are defined as:
Summary. To summarize, our solver solves problems of the form (1), for which
where each of 𝐽 𝑖 is one of the objectives, enumerated above, combined with weights. The equality constraint 𝐹 (𝑢) = 0 is given by ( 13), the static elasticity equation with smoothed contact and friction.
The inequality constraints are of the form 𝐽 𝑚 < 𝐽 𝑏𝑛𝑑 𝑚 where 𝐽 𝑚 are some of the remaining objectives for which upper bounds are imposed (e.g., a bound on 𝐿 𝑝 norm of stress). These are imposed as soft constraints by adding 𝑖 1
In this section, we describe our approach to computing the discretization of the problem and the gradient of the functional projected to the constraint space 𝐹 (𝑢) = 0 (i.e., space of displacements satisfying the elasticity equation with contact and friction).
For each objective, regardless of it is used as an optimization target, we need to compute the gradient with respect to 𝑝. As explained in Section 3.1, we use shape derivatives, which are computed using the solution of an adjoint equation. More specifically, the solution is used to construct a form 𝑑 𝐽 𝑑 on the surface, which, when applied to changes 𝛿𝑝 of vertex positions on the boundary, produces the (linearized) change of the functional 𝐽 :
where 𝑢 𝑡 is the solution of the elasticity equations obtained on the domain with deformation 𝑡𝛿𝑝. The theory of shape derivatives is well-established, and we refer to [Bonnetier and Dapogny 2020].
Here, we only present a brief summary of the steps for computing these derivatives. A detailed derivation can be found in the supplementary material. Our formulation follows [Panetta et al. 2017], in that we consider volume, rather than surface integrals, which proved to be critical for derivative accuracy, at a moderate additional cost, compared to the elasticity solves involved.
Adjoint equation. The first step in computing 𝑑 𝐽 is to solve the adjoint equation, which in our case has the following form:
for any 𝜓 , where the terms 𝑇 ′ 𝑋 , corresponding to the contact and friction smoothed penalties are of the form ∫
• 𝜓𝑑𝑆, and 𝜌 is the unknown. We observe that this is a linear elasticity equation which is relatively inexpensive to solve compared to the primal nonlinear elasticity equation. The quantity 𝜏 is defined as the sum of 2𝑒 ′ 𝜎 : 𝐶 for all objectives of the form (2) where 𝑒 (𝑠 ( ū), 𝑥) is a function on a stress measure, as our case. Then, 𝑒 ′ means the partial derivative of 𝑒 with respect to our stress measure.
Shape derivative. If the boundary deformation is defined by a set of basis functions 𝜆 𝑚 , with the deformation expressed as 𝑣 = 𝑚 𝛿𝑝 𝑚 𝜆 𝑚 , where 𝛿𝑝 is the vector of changes of our variables, it is possible to express the shape derivative as a dot product
where 𝑆 [𝑢, 𝜌] is a vector of the same length as 𝑝 of vector valued functions depending on 𝑢 and 𝜌. The expressions for these functions are included in the appendix, and derived in the supplementary document. We emphasize that no discretization, other than discretization of the deformation of the boundary was performed so far.
The vector 𝑆 [𝑢, 𝜌] computed by numerical integration of FEM solutions of two elasticity problems for 𝜌 and 𝑢 is the gradient of the functional 𝐽 with respect to 𝑝; this is what we use in the optimization process as discussed in Section 5.
To make our discussion of the formulation complete, we summarize the discrete form of the problems we solve. The exact expressions are straightforward to derive but tedious (please see the supplementary material). We use quadratic Lagrangian elements on tetrahedras for discretization of (13) to obtain a system of the form
The first two terms come from the standard elasticity equation, and the rest correspond to the remaining terms in (13). The remaining terms, while simple, are not linear, so the system after discretization is a general algebraic system, and requires a non-linear solver. We also derive expressions for the Jacobians 𝐷𝑁 and 𝐷𝐹 , which are needed for efficient optimization.
In contrast, the adjoint PDE ( 23) is linear and has the form
i.e. involves exactly the Jacobians of the constraint functions. The right-hand side 𝐷 is expressed in terms of 𝜏. Formulas for the entries of these matrices are provided in the supplementary document. Finally, once ū and ρ are available, the coefficients of the shape derivative are computed from these values by integration, following formulas in the appendix, using numerical quadrature described in Section 5.
Input. The input to our algorithm is a collection of tetrahedral meshes, with some parts of their surfaces in contact, and with boundary conditions applied to other parts of surfaces. Each object can be designated as deformable or rigid.
Our overall method is conceptually straightforward (Algorithm 1). The functions used in the pseudocode are:
• ElasticitySolve solves the nonlinear elasticity system (3), to obtain ū; • AdjointSolve solves the adjoint system (26) to obtain 𝜌; • DiscreteShapeDerivative, given ū and 𝜌, computes the gradient 𝑆 [𝑢, 𝜌] with respect to vertex positions on the boundary using ( 24). • Converged is the outer iteration stopping criterion discussed below.
The inner loop works on a fixed connectivity for Ω(𝑝), and is close to the standard BFGS algorithm: at each step, a descent direction is computed, and a line search is performed to determine our step size. There are three important differences: (1) we check for any inversions of tetrahedra and choose a step that maintains a bound on mesh element shape quality; (2) after each update of the boundary vertices, we call the SLIM smoothing algorithm [Rabinovich et al. 2017] on interior vertex positions 𝑝 𝑖𝑛𝑡 with boundary vertices 𝑝 fixed, to move the interior vertices so that the quality of the mesh is improved; (3) with a valid mesh, we run simulation and if it doesn't converge in a maximum number of iterations, we reduce step by half in the line search.
If the step becomes too small, the inner loop is terminated, and the domain is remeshed in the outer loop.
A natural stopping criterion for the algorithm consists of three parts: (1) the objective reduction obtained in a step of (outer) iteration is below a threshold 𝜖 𝑟 computed relative to the initial objective value; (2) the step size of the line search falls below a threshold (3) the maximal number of iterations 𝑚𝑎𝑥_𝑜𝑖 is exceeded. Due to remeshing, the energy however may oscillate slightly, and for robust behavior we require that sufficient energy decrease does not happen over 𝑚 steps. While many other options are possible (e.g., relative or absolute gradient norm threshold), we consider the rate of change in the objective to be most appropriate: our goal is to obtain a reduction in the objective, and a slow rate of reduction indicates that optimization will not improve the target by much more in a reasonable number of iterations; this happens either due to being close to a local minimum value or step size going to zero for geometric reasons, typically thin regions, resulting in distorted elements.
We discuss the choices for 𝜖 𝑟 and 𝑚 in Section 6, as well as examine convergence for specific test cases.
For remeshing, we use Triangle in 2D, and CGAL and fTetWild [Hu et al. 2020] for 3D tetrahedral meshing. More information on meshing is presented in the Appendix E.
and 𝑁 𝑣 is the neighborhood of 𝑣. The value of power 𝑝 can be adjusted to obtain smoother surfaces at the cost of less optimal shapes; we use value 2 for most cases, increasing it to 4 for some objects. This term is scale-invariant and pushes the triangles/tetrahedra of the mesh toward equilateral.
Elasticity solver. In ElasticitySolve we use the standard Newton's method with line search to solve the nonlinear elasticity system with contact. We consider the simulation solved when residual is lower than a given tolerance. For 2D, we use 10 -10 for the tolerance, while the value of 10 -8 is used for our 3D examples. Moreover, for solving linear systems at each iteration of the Newton's method, we use CHOLMOD (for frictionless scenarios) and UMFPACK (when friction is present) in 2D. In 3D, we use MKL Pardiso library.
Adaptive Quadrature. An extremely important aspect of our implementation is the quadrature used to compute the integrals in the FEM system discretization, as well as in the shape derivative coefficient formulas. For problems with friction, precise computation of these integrals proved to be very important. At the same time, due to functions like ℎ 𝜂 , the functions we integrate, while smooth, have higher derivative discontinuities. We use a combination of adaptive refinement on triangles and high-order (order equals to 10) Gaussian quadrature to integrate all functions to high precision. For more information, see Section D in the Appendix.
Convergence behavior. Figure 5 shows the objective as a function of iteration number for a 2D connector (Figure 9) and a 3D stool (Figure 24). For both cases, we optimized 𝐿 𝑝 norm of stress, while keeping volume below or equal the initial value. For the connector, during the last 200 iterations, the energy decreased only by 0.87%, while for the stool, the minimum decreased by 0.27% in the last 30 iterations. We provide additional data in Section 6.
Observe that there are increases in energy at some of our iterations, corresponding to the remeshing of the current solution. However, after that the objective function quickly decreases, which motivates our formulation of the convergence criterion.
Simulation Validation. As we use approximations to the standard physical models of contact and friction we evaluate the accuracy of these models. We use two examples. The first one is box in contact with a rigid surface. The second one is joint in contact with a holder as shown in Figure 6.
The results for both cases are largely the same. Thus, we will only show the results for box in contact with a rigid surface in this section and leave out the results for joint in contact with We first change the value of 𝛼 and 𝜂 The reference we use for comparison is obtained by simulation with 𝛼 = 𝜂 = 10 -6 . We set the smallest value of 𝛼 to be 10 -7 and scale it by the power of 2. The relative value of the difference in relative displacement |𝑢 -𝑢 𝑟𝑒 𝑓 |/|𝑢 𝑟𝑒 𝑓 | is shown in Figure 7.
Then, we test the same scenarios but with friction at contact regions. First, we consider the dependence of accuracy on 𝜂 𝑛 , using 𝜂 𝑛 = 10 -4 as the reference value, as for friction we typically need more smoothing in the constraints. To test how friction coefficient affects simulation, we add in a small horizontal load that is 10% of the vertical load for the box example. Test with 𝜇 ranging from 0.1 to 0.16 and increase 𝜇 every time by 0.01 for the two examples. We compare the tangential displacement of the examples with respect to different 𝜇 to the tangential displacement of the examples with respect to 𝜇 = 0.1 and take the ratio. We get the plots shown in Figure 8.
Effects of optimization targets and constraints. We use a simple twodimensional connector example to demonstrate the effects of various optimization objectives and constraints (Figure 9) We compare stress distribution and maximal stress in different cases.
For most of the examples related to this scenario, we used a similar optimization configuration, running a maximum of 1000 iterations and remeshing at least every 100 iterations (𝑚𝑎𝑥_𝑖𝑖 = 100 and 𝑚𝑎𝑥_𝑜𝑖 = 10).
The most radical difference is between minimizing the volume while bounding stress and minimizing maximal stress.
We also compared the scenario of minimizing maximal stress with and without a bound on volume (equal to the initial volume). In this comparison, stress results were actually very close to each other, with a slight advantage to the version with the constraint.
Multiple loads. Our framework allows optimizing for multiple separate acting forces, meaning that the energy related to stress will be a combination of the values from each separate scenario. An example is shown in Figure 10, where forces of the same intensity are applied to the right and to the left of the top piece. The stress results of the optimized shape are very close to each other and the final shapes present considerable symmetry, even though we are not enforcing any geometry symmetry through constraints. Differences in material properties. Figure 11 shows how the results of optimization are affected by combining a highly flexible and a stiff material on our two different pieces. (In both cases the target is the baseline case of minimizing stress with no constraints). We observe that the stress reduction is similar, although shapes required for this are substantially different.
In the next comparison (see Figure 12), we observe that in the presence of friction (𝜇 = 1.0), we achieve a similar maximum stress value with the assembly not needing a more extreme protrusion to stay together, due to additional forces resulting from friction. This is demonstrated in Figure 12. The rightmost example shows the drastic effect of fixing only one part. Assembly/disassembly and parallel alignment constraints. We show the effects of our connection-related constraints presented in Section 3.5.
Figure 13 shows the baseline optimization (minimizing maximal stress with a bounded volume) compared to different cases when our assembly constraint in all cases, and different target parallel alignments. Observe that, as expected, the contact area aligned with the disassembly direction increases.
All three examples are possible to assemble in vertical direction without deformation. Also, the parallel alignment objective guarantees that the desired proportion (30% and 50%) of the contact walls are parallel to the input direction.
Figure 14 shows the baseline case compared to one using a disassembly constraint with a very small tolerance. Compared to the the baseline optimization, two pieces do not detach when the disassembly constraint is applied. Another example of the importance of this constraint is shown in Figure 20.
Assembly/disassembly tradeoff. Finally, we also studied the effect of choosing different balance of assembly and disassembly objectives, in a setting when these counteract each other. Consider the case of Figure 15, where a force is pulling the top part of the connector up and we initially optimize the shapes using the same weight (100.0) for both terms and weight 1.0 for stress. Then, we reran our optimization using three lower weights for the assembly term. Figure 15 shows that, by reducing assembly's importance, almost no movement is observed in the last shape where assembly weight equals 0.1. When the disassembly term is dominant, the optimized shape has small protrusions on both sides that keep it firmly in place, at the expense of greater effort required for (dis)assembling it along the (dis)assembly direction. Table 1 presents a comparison of our main energies in this study. Note how assembly objective values increase when it's weight is reduced, while, at the same time, disassembly and stress values decrease considerably.
To validate our framework, we computed optimized shapes for a number of 2D and 3D realistic scenarios. In all illustrations in this section, parts of the object boundary with Dirichlet conditions are shown in red, and contact regions with external supports in yellow. We use 𝑝 = 20 for stress optimization in all examples. Except for 2, where we add information about the simulation, as well as objective functions used in each example, instance size and initial and final maximum values of stress.
Lever. In this example shown in Figure 16, one part of the boundary of the black piece of the assembly is clamped to the table and in simulation assumed to satisfy the Dirichlet boundary condition; the part touching the vertical wall has an RD contact condition. Maximal stress is optimized, with no constraints imposed. As seen in Figure 16, we are able to reduce maximum stress more than 14 times compared to the initial shape. The optimized shape naturally evolves into an interlocking assembly that can support a far higher load without large deformation.(see the video in supplementary material). We fabricated resulting shape (the process is described in Appendix C, and performed a stress test on it, by loading it with increasing weight. Our setup is shown on Figure 17. The unoptimized lever breaks with a weight of 400g, whereas the optimized lever can endure more than 6400g, consistently with simulated results.
Bridge. Our second 2D example is a simple bridge model (Figure 18), with Dirichlet conditions on fixed rectangular parts, 3 optimizable parts, with boundaries of the supports partially fixed, and DD contact conditions between parts. Due to symmetry, we run the optimization on one half of the shape, resulting in 3 times stress reduction.
Table 2. Summary of experiments. In constraints column, "/" means that the weight is modified after each remeshing and varies between these values linearly. Column Relative gradient norm presents the ratio of the final gradient norm and the initial one from the first iteration. The last column (Final 𝐿 𝑝 -norm reduction) shows the reduction in the stress objective due to the last outer iteration, expressed as a percentage of the initial objective value. * For both bridge and coat rack (hooks), running 5 more outer iterations reduces stress to 1.80 and 346.51, relative gradient norm to 2.6e-2 and 2.81e-2, and final 𝐿 𝑝 -norm reduction to 0.95% and 0.00%, respectively. For fabrication, we used optimization results at outer iteration 5 for lever, hook, crane and at iteration 10 for bridge and coat rack. The left and right support parts of the bridge were clamped to the table and weights were attached to the middle of the bridge. The experiments show that unoptimized bridge breaks at the left or right assembly with a load of 1500g, whereas the optimized bridge can hold more than 5800g of weights. See photo of our bridge in Figure 19.
Hook. Figure 20, where we built an optimized hook that fits on the top of a door, with a single large contact area Γ 𝐶 , and 𝑅𝐷 boundary conditions on that part. In this scenario, we use multiple separate loads, one emulating, e.g., a coat hanging from the hook, and the other a force exerted when it is taken off the hook. In addition to multiple loads, this example also demonstrates the importance of disassembly energy: we make sure that accidental push from below does not result in the hook getting detached. The disassembly direction is shown in Figure 20. We perform stress minimization with a volume bound equal to the original volume, and compare the results with and without the disassembly constraint. Note that when the hook optimized without disassembly constraint is loaded, the left part lifts up, resulting in detachment from the support (Figure 20, middle) It remains stable once the disassembly constraint is enabled (Figure 20,right).
Figure 21 shows the experiments with 3 versions of 3d printed hook (unoptimized and optimized with and without disassembly Tower Crane. Our last 2D example is the tower shown in Figure 22, demonstrating an assembly consisting of four pieces with a Dirichlet condition on the bottom one, and DD contact conditions on the remaining parts.
In this example we demonstrate how our method can be combined with a simple ESO-like topology optimization technique (e.g., [Huang and Xie 2010]). to decrease the weight beyond what is possible with shape optimization only. We used a filtering technique to remove triangles with average stress lower than 10% of the maximum value. The filtered mesh was then again optimized, obtaining the result shown on the bottom-right. The laser-cut tower model is shown in Figure 23. In our experiment, the base of the tower was clamped to the table. Experiments show that the unoptimized tower breaks with 700g load, whereas the optimized tower can sustain more than 4600g load.
In the next 3 examples, we apply our algorithm to 3D assemblies of parts made of different materials: wood (𝐸 𝑤𝑜𝑜𝑑 = 10000,𝜈 𝑤𝑜𝑜𝑑 = 0.3), MDF (𝐸 𝑀𝐷𝐹 = 4000,𝜈 𝑀𝐷𝐹 = 0.25) and PLA plastic (𝐸 𝑃𝐿𝐴 = 3500,𝜈 𝑃𝐿𝐴 = 0.36). Moreover, in these examples, only some pieces of the shape (those made of PLA, which can be 3D printed) were optimized. Stool. Our first example was a stool assembly with a fixed MDF top, fixed wood legs and connectors made PLA (Figure 24). In this example, use volume The stress for optimized and noncan seen on Figure 25.
The optimization reduces stress 2.5 times without changing the total volume. If 10% increase in the volume is allowed, the stress decreases to 0.88, and to 0.85 at 20%.
Our second example is a bench with fixed parts (seat,legs) made of wood and optimizable PLA connectors ( Figure 26). We consider two different leg shapes and show how connectors are optimized for each. We observe that, depending on the type of legs, we obtain a very different level of stress reduction: 30% in one case, and 2.5x times in the other (see Figure 27).
Coat rack. Our final multimaterial example is a coat rack made of wood and plastic (Figure 1). Differently from both previous examples, this object is modular, in the sense that you can always use a longer cylinder and add identical plastic parts to increase the height and the amount of hangers in your object. Here, we optimized our two different plastic parts for separate sets of loads, the leg connector and the hook. The leg connector is optimized with (fixed) legs in contact with the ground and the load applied at the top (Figure 28).
The hook attachment was optimized with loads applied to hooks. As shown in Figure 29, our framework was able to lower maximum stress more than 3 times.
We 3D printed the result using PLA (Figure 30). For our last set of experiments, we implemented a framework for generating pipe assemblies from graphs, where each vertex of our graph becomes a sphere and each edge becomes a pipe. See an example on Figure 31. We used this tool to generate three different simple instances of our problem: a tent, a truss-like bridge and a dodecahedron. Below we detail more about each of these assemblies.
Tent. This object (Figure 31) is composed of 5 ball connectors and 8 pipes. A load is applied on the top of the structure pointing down and the bottom part of the object is in contact with the ground. We minimize stress, while keeping the same volume impose the assembly constraint. We were able to reduce stress by a factor of 2. We again ran our framework for reducing the volume of the shape, while keeping the stress bounded by the result of the previous optimization for stress (Figures 31 and32). Then we add disassembly and parallel alignment objectives and the results are shown in Figure 33. Note that we optimize both the bars and connectors.
While the resulting structure can support this particular load well, it is likely to fall apart if any lateral load is applied, as resulting connectors are very shallow. Adding disassembly and parallel alignment constraints produces more resilient connectors (Figure 33). When using only disassembly constraint, the final stress (59.1) obtained was similar to the baseline case (59.3), while adding parallel alignment constraint to make at least 20% of the connection surface (in each connector) parallel to the disassembly direction increased stress to 70.8.
Optimized: shape, as presented in Figure 35, where our boundary conditions are also shown. We use a similar optimization setup as for tent: first, run stress minimization with a volume constraint, then minimize volume with a stress constraint, and including fixed assembly, disassembly and parallel alignment constraints. Notice in Figure 35 we were able to reduce stress by almost 10 times, while reducing volume in 12.6% compared to the initial dodecahedron. The emerging twisted shape for the bars are optimal for load support.
We also investigated the effect of optimizing only the connectors (balls) the dodecahedron example. Again, we first minimize stress and then, in the second step, minimize the volume. We were able to reduce stress by around 28% to 172, and reduce the volume by 22%. The result can be seen in Figure 36.
Truss. The initial shape and the optimized result are shown in Figure 37. Similar to the dodecahedron, we optimized for strength (keeping initial volume) and then for volume (keeping stress close to its minimum level). We were able to reduce stress by more than half and volume by 20%. Notice that the optimized connectors are smaller than original ones while keeping deep connections. This result can be obtained thanks to the disassembly and parallel alignment energies. Termination and convergence. To understand the convergence behavior of the method, we run our optimization for a fixed number of iterations which we aimed to set high enough for obtaining maximal possible improvement. Table 2 lists several indicators of convergence termination: number of actual iterations performed before stopping, relative gradient norm, and per-step objective reduction, expressed as a percentage of the initial objective value.
In the cases when the optimization reached the maximal number of steps (10 outer iterations for 2D and multi-material examples and 50 for 3D examples), the result is close to a local minimum, or the maximal number of iterations is insufficient. In most cases, relative gradient norm is 5% or significantly less of the original, indicating that the result is close to local minimum. In one case, (Bridge) additional iterations also decrease the relative norm below 5%.
In several cases (hook and two bench variations) the optimization stops for geometric reasons (in some cases, topology changes are needed to achieve an optimum); in this case, the optimization cannot reach a local minimum by additional iterations, although the gradient reduction is also high in most cases.
Examining how the objective changes with iteration suggests that values 𝑚 = 2 and 𝜖 𝑟 = 1% are adequate for all examples (considering the criterion discussed in Section 5 for our stress objective), except truss-like structures, where values 𝑚 = 5 and 𝜖 𝑟 = 0.5% are more suitable due to higher oscillations in the objective early in the optimization. We emphasize that many other convergence criteria can be used; determining an optimal one is not our focus.
We have described a framework for shape optimization with contact and friction, using a smoothed, penalty-based model, and applied it to optimization of assemblies in which parts are held together by contact and friction forces. We demonstrated that in our framework a number of functionals can be optimized reliably, producing significant improvements e.g., in stress concentrations or volume, while maintaining various types of constraints. We validated the results with computational experiments in two and three dimensions and with qualitative experiments using fabricated objects. The experiments confirmed significant improvements in strength. In addition, since we support nonlinear solves, our model can be easily extended to a variety of functionals and constraints. The code for this project will be open-sourced.
Limitations. Our algorithm may (1) not reach a minimum or (2) find a local minimum that does not improve the objective much. In the first case, the optimization stalls: no progress is possible for one of the following reasons: (a) nonlinear solve does not converge (e.g, if the system has no static equilibrium or bad-quality elements not allowing to achieve residual tolerance in max number of iterations); (b) inversions of elements lead to a severe line search step restriction, as mesh quality improvement is not guaranteed to succeed; (c) meshing with tolerance using fTetwild may produce volume meshes for pieces with non-manifold surfaces, while surfacepreserving remeshing does not improve quality enough. In case (2), the algorithm may reach a local minimum without significant decrease of the objective value, while other minimal with lower values exist. Case (1) was the most observed in our experiments, occurring in our hook and bench examples due to narrowing regions during optimization; note that we still obtain a significant stress reduction.
In addition, we use linearized elasticity model that does not account well for large deformations. As the contact and friction terms in our equations are already non-linear and we use a nonlinear solver, extension to nonlinear elasticity is likely to be straightforward. While the functionals we have introduced in this work perform quite well for a number of tasks, there are restriction on the type of connections they can capture. E.g., snap connections are not handled easily with these constraints only. Expanding the range of objectives is an important direction for future work. The method can be considerably sped up but using parameterized geometric models, or model reduction, as well as more efficient solvers.
A FORMULAS FOR THE ADJOINT EQUATION TERMS.
where the first term is related to the elasticity equation and the second to contact and friction terms. In turn, 𝑆 𝑐 𝑚 = 𝑆 𝐶 + 𝑆 𝑆 + 𝑆 𝐶𝐹 + 𝑆 𝑆𝐹 , the sum of terms for different types of contact forces. For these terms we have the following expressions, where we consider uniform surface load 𝑔 = 𝐺 |Γ 𝑁 | :
We have validated the performance of our system on several simple assemblies fabricated with 3D printing and laser-cutting. These are shown in For the lever and hook example as shown in Figures 17,21, the fabrications are done with 3D printing (Ultimaker 3), using PLA material. For these experiments, we used Ultimaker Cura for slicing. We used fine printing setting, with 0.1mm layer height and infill density equal to 100%. The general print speed was also altered to 40mm/s and reduced to 12𝑚𝑚/𝑠 (18𝑚𝑚/𝑠) when printing the inner (outer) walls of the object, in order to have a more precise printing of connections. In addition, to guarantee that objects connect, we offset the polygon shapes by 0.05mm, reducing them to account for the printing thickness.
For the bridge and tower examples (see Figures 18,23), the experimented objects were obtained through laser-cutting 1/8" (3.175mm) acrylic boards, using Epilog Mini cutter. Notice that laser-cutting removes material, so we perturbed the models to account for that, offsetting the polygonal shapes by 0.05mm.
For our 3D examples, we also 3D printed our coat rack optimized model, shown in Figure 30. We used the fast setting on Cura software, with 0.2mm layer height, and reduced the speed to 50mm/s.
Even though we assembled with wood dowels, which can carry imprecision in the measures, the pieces still fit tightly.
In different steps of our algorithm, during the simulation and also shape derivative computation, we need to perform the computation of integrals. For simpler cases using linear elasticity, this is not a challenge and we can simply use Gaussian quadrature, for example. However, when adding penalization due to normal contact and friction, the choice of how to perform these computations becomes very important.
As an example,consider the discretized normal contact force. Here, we assume we are using quadratic FEM elements. Our discretized contact function per contact element Γ 𝑒 𝑆 is: where 𝛼 is a constant, 𝑛 is constant for each element, 𝑝 (𝑥) maps each point 𝑥 into barycentric coordinates of the element, 𝜑 corresponds to the quadratic baricentric basis and 𝑢, the displacement, is also quadratic. ℎ 𝜂 is the same as previously defined:
Notice that 𝑢 (𝑝) • 𝑛 defines the region of the piecewise function. To find the boundary between each of these regions, we would need to solve 𝑢 (𝑝) •𝑛 = ±𝜂. Notice that this corresponds to a quadratic 2D equation and the solution corresponds to a conic section: parabola, hyperbole, ellipse or circle. We have 2 equations, so our original domain is partitioned by a maximum of 2 conic sections. Since the two equations differ only by 𝜂, we expect that the regions boundary are going to have very similar shapes. See an example in Figure 38: Here, 𝑢 𝑗 corresponds to the discrete displacement at each of our 5 nodes of our triangle element. Using Matlab's integral2 function (with absolute tolerance set at 10 -10 ), we can compute 𝑗 ′𝐷 𝑁 ,𝛼 to obtain -0.329516. Using (order 10) Gaussian quadrature, we obtain -0.438527. These errors can accumulate and generate problems in convergence for our simulation.
We implemented then an adaptive approach for computing integrals in our framework. First, we define the scheme to be used by the regions of our triangle. We evenly sample the triangle using 𝑚 points, identifying the regions where these points are contained. In our case, we used 𝑚 = 10. If not all 𝑚 points are from the same region (easily verified by computing 𝑢 (𝑝) • 𝑛 value), we split the triangle into 4 pieces of same area and repeat this process recursively, until a maximum depth is reached. See an example on the right side of Figure 38 In our case, we identified that depth equal to 5 guarantees good results. Inside each of the final triangles, we use (order 10) Gaussian quadrature to compute the integral. Using this method, for the same example as above, we obtained result -0.329560, which is accurate up to 4th digit after the decimal point when compared to the numerical result of the integral2 function.
We use triangle meshes in 2D and tetrahedral meshes in 3D for our shapes.
For 2D, we define shape outlines using Adobe Illustrator, encoding boundary conditions using color; shared boundaries between parts are represented as single polylines. We then transform the geometry in the resulting SVG file into the input format for Triangle [Shewchuk 1996], removing the duplicate edges defining the contact area between different pieces. After obtaining the mesh for the interior, we split the triangle mesh into part meshes. At the end of each inner loop in Algorithm 1, we extract the polygonal boundary of the current solution and remesh it using Triangle again.
In 3D, the input to the meshing process is a collection of surface meshes of individual parts. We assume that in the contact zone, the distance between meshes is small, (in our specific examples we use boolean operations to produce surface for parts, with connectors obtained, e.g., by adding a cylinder to one part, and subtracting it for another, which produces close surfaces. After that, we generate the volumetric mesh using fTetWild [Hu et al. 2020] with surface meshes from the previous step as input. fTetWild has a threshold 𝜖 𝑚 within which close parallel triangles of a surface mesh are merged; we set it to a sufficiently large value to obtain a single contact surface. For most cases, we used the default value of 1e-3 for 𝜖 𝑚 , but increased the value to adjust to different gaps coming from the surface mesh generation.
As in 2D, we apply remeshing to keep a good mesh quality after a sequence of optimization steps. In 3D, remeshing is done in different ways: (a) preserving mesh surface, which is done with CGAL's package for tetrahedral remeshing [Tournois et al. 2021]; or (b) remeshing the boundary as well, by reapplying fTetWild. Choosing between this two options is also a parameter in our framework. We used option (a) for most of our scenarios, but for truss and dodecahedron example, boundary remeshing was used at every 5 remeshing operations. If remeshing the boundary with fTetWild fails, e.g., due to generation of pieces with non-manifold meshes, option (a) is used instead.
ACM Trans. Graph., Vol. 40, No. 6, Article 269. Publication date: December 2021.
ACM Trans. Graph., Vol. 40, No. 6, Article 269. Publication date: December 2021. Optimizing Contact-based Assemblies • 269:7
ACM Trans. Graph., Vol. 40, No. 6, Article 269. Publication date: December 2021. Optimizing Contact-based Assemblies • 269:11
ACM Trans. Graph., Vol. 40, No. 6, Article 269. Publication date: December 2021. Optimizing Contact-based Assemblies • 269:13
ACM Trans. Graph., Vol. 40, No. 6, Article Publication date: December 2021.
ACM Trans. Graph., Vol. 40, No. 6, Article Publication date: December 2021. Optimizing Contact-based Assemblies • 269:17