- Email: [email protected]
Generative Design Conversion to Editable and Watertight Boundary Representation
Computer-Aided Design 115 (2019) 194–205
Contents lists available at ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Generative Design Conversion to Editable and Watertight Boundary Representation✩ , ✩✩ ∗
Martin Marinov , Marco Amagliani, Tristan Barback, Jean Flower, Stephen Barley, Suguru Furuta, Peter Charrot, Iain Henley, Nanda Santhanam, G. Thomas Finnigan, Siavash Meshkat, Justin Hallet, Maciej Sapun, Pawel Wolski Autodesk, Inc., United Kingdom, USA, Australia, Poland
article
info
Article history: Received 24 April 2019 Accepted 4 May 2019 Keywords: Generative design Topology optimization Mesh segmentation Mesh parameterization Mesh to B-Rep conversion Solid models
a b s t r a c t We present the first, to our knowledge, fully automatic method for the conversion of a general generative design to a watertight B-Rep composed of the exact residual geometry of the input solids and editable T-NURCC surfaces. The design can be synthesized by any combination of both additive and subtractive in-synthesis processes, and constrained by any number of keep-in and -out solids. Our method requires only the input generative setup solids, rudimentary data from the solve, and a boundary mesh of the optimized design. Leveraging the generative solve data, we augment the boundary mesh with ‘‘incidence’’ attributes linking it to the input solids, and partition it into organic and incident regions. The organic regions are parameterized and approximated with T-NURCCs. Then, their boundaries are ‘‘pulled’’ into the adjacent input solids to construct clean intersection curves. Finally, the organic surfaces are combined with the input solids to compose a watertight solid B-Rep of the generative design. Without such an automatic conversion, users seeking to exploit generative designs in downstream workflows (e.g. setting up assembly constraints, simulation analyses, performing parametric edits) have to manually convert the design to a B-Rep, slowing down the path to manufacturing and increasing the product time to market. As a result of our work, engineers can now quickly obtain useful CAD models of any number of generative designs computed by a variety of algorithms, settings and iterations. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction Generative design is a modern approach to synthesizing design alternatives by exploring the entire design space to achieve objectives such as minimal structural compliance, given a set of functional requirements and constraints based on geometry, material specifications, and manufacturing considerations. A precursor to this technology is topology optimization [1], a niche tool employed by a small number of expert users using specialized simulation software and know-how along with a significant budget. The availability of cloud as an ‘‘infinite computing’’ platform ✩ This paper has been recommended for acceptance by Pierre Alliez, Yong-Jin Liu & Xin Li. ✩✩ One or more of the authors of this paper have disclosed potential or pertinent conflicts of interest, which may include receipt of payment, either direct or indirect, institutional support, or association with an entity in the biomedical field which may be perceived to have potential conflict of interest with this work. For full disclosure statements refer to https://doi.org/10.1016/j. cad.2019.05.016. ∗ Corresponding author. E-mail address: [email protected] (M. Marinov). https://doi.org/10.1016/j.cad.2019.05.016 0010-4485/© 2019 Elsevier Ltd. All rights reserved.
and novel manufacturing methods like additive manufacturing have opened up the possibility of using generative technology as a new approach to design. This capability is starting to become available in CAD packages such as [2], enabling users to leverage generative designs in diverse areas including aerospace, automotive, consumer products and architecture. The generative design geometry constraints are defined as solid B-Reps and are placed in a R3 optimization domain to avoid interference (keep-out) and specify attachment joints (keep-in) with other parts in a CAD assembly. However, many generative solvers, e.g., [3–5], do not operate directly on the exact geometry constraint B-Reps. Instead, B-Reps are sampled and replaced with volumetric representations such as level sets or tetra-/hexahedral meshes, which are significantly more convenient and efficient for physical simulations and material synthesis. Inevitably, this conversion introduces a loss of surface precision due to the limitations of the target representation (Fig. 1). While this loss is negligible in the context of the generative computations, in general it prevents the output of a generative design that accurately represents the input solid surface geometry away from the organic material synthesized by the solver.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
195
B-Rep by leveraging the generative solver input and representation to:
• embed the exact input solid boundary surfaces where the design coincides with the input,
• approximate everywhere else the design boundary with globally smooth, editable ‘‘organic’’ surfaces, and
• join all surfaces to form a generative design output B-Rep. Our method is applicable to generative designs that represent arbitrary combinations of material additions and subtractions, and leverages the solver geometry representation. 2. Related work
Fig. 1. A fragment of a mesh extracted from the level set representation of a generative design. The yellow segment represents the overlaid semi-transparent input solid. The accuracy loss is especially apparent near sharp edges and holes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Subsequently, when integrating the design within the complete model assembly, these surface inaccuracies can cause issues such as part interferences and manufacturing deficiencies (Fig. 3). To resolve these inaccuracies, our fully automatic method (Fig. 2) converts a generative design into an editable, watertight
General B-Rep surface reconstruction from meshes and point clouds is a highly prominent and long-standing research topic with particularly important applications in the reverse engineering of existing physical objects [6–9]. Reviewing this work is beyond the scope of our paper, so we refer to several recent comprehensive surveys: [10–12]. More relevant to our method is research work which focuses on partitioning and approximating a mesh with piece-wise smooth surfaces [13–19] and networks of G1 continuous surface patches [20,21]. However, none of these techniques are known to automatically produce results fulfilling the requirement to represent accurately the residual input solid surfaces in the generative design. Approaches approximating the entire mesh with a globally smooth subdivision surface [22–29] appear at first sight directly applicable to our problem, especially if equipped with a global
Fig. 2. Our generative design conversion in brief. Top: Input solids (left) and 0-isosurface mesh (right), partitioned by our method into incident and organic parts. Bottom: Editable T-NURCC surfaces constructed for the organic part (left) and output watertight B-Rep composed from the input and organic surfaces (right). Note that before composition, T-spline patches are merged where possible to produce the simpler patch layout on the right.
196
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Fig. 3. Left: a rod (in blue), designed to fit through a couple of input solids (in yellow). Middle: a general surface reconstruction violates the generative constraints causing interferences (in red). Right: The rod fits properly when assembled with our generative B-Rep, which reuses the input geometry. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
quad parameterization algorithm such as [30–37]. However, the shortcomings of such unrestricted global approximation schemes for generative design models are easily apparent (Fig. 3). It is tempting to ‘‘fix’’ such defects by applying Boolean operations to remove interferences (subtract keep-out solids) and add missing volume (unite keep-in solids). In practice though, Booleans on nearly overlapping surfaces are very unreliable, which leads to outright failures and B-Reps that are problematic for downstream modeling (Fig. 4).
In the CAD industry, some commercial products, e.g. [38,39], offer capabilities to build a B-Rep model using an input mesh as a geometry reference. However, even for moderately complex shapes, these products require significant manual effort to meaningfully align the smooth surfaces with the processed mesh geometry. For the complex, feature-rich shapes found in most generative designs, such a manually-intensive approach is extremely time-consuming and often cost-prohibitive. Instead, our method is fully automatic and does not rely on user guidance. While specialized solvers such as [40] can process input solid B-Reps directly and output a consistent B-Rep instead of a mesh, the shapes they produce are significantly less varied (e.g. lattices) and only applicable for specific physical problems. In contrast, our approach is applicable to generative solvers optimizing for arbitrary physical criteria, producing a variety of unique and useful shapes. Another alternative to our method is the direct embedding of generative design meshes in the output B-Rep instead of smooth, precise surfaces [41,42]. While this approach can be used for some downstream editing and processing, most modeling operations available in the current generation of B-Rep kernels are not applicable to mesh geometry, thus significantly limiting the usefulness of B-Reps with mesh embedding. 3. Overview 3.1. Input data The primary input to our method consists of the generative solver input solid B-Reps {Si } and output mesh M representing the generative design G. To link M to {Si }, we can also process the intermediate solver representation (level set or volumetric mesh). Alternatively, the solver representation accuracy error, e.g., the level set voxel grid edge length, can be supplied instead. {Si } designate sub-spaces of the optimization domain. We use the standard notation ∂ Si and Sio to denote the solid boundary surface and interior. {Si } can be of the following types:
• Keep-in: Si is completely contained, i.e., Si ⊂ G. • Keep-out: Sio is outside G, ∂ Si can overlap. • Seed: Si is fully or partially inside G.
Fig. 4. Top: An organic surface approximating globally the generative mesh in Fig. 2 violates the keep-in and keep-out constraints as shown in Fig. 3. Bottom: Attempts to repair the violations by Boolean operations fail for the yellow keepin solid and produce poor outcomes everywhere else. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
We unite overlapping solids of the same type, so without loss of generality, Si ∩ Sj = ∅ for every i, j of the same type. We denote the union of all input solids (regardless of type) as S = ⋃ Si . M = (V, E, T) is a closed 2-manifold mesh approximating G, where V = {Vi }, E = {Ei }, T = {Ti } are correspondingly the mesh vertices, edges and triangles. If a non-triangular mesh is provided, we triangulate any quads and general polygons.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
197
3.2. Workflow The main stages (Fig. 2) and algorithm steps of our method are listed in order of execution below: 1. Partition the generative design mesh into two regions: incident and organic. The former corresponds to the input solid surfaces and the latter to the generative synthesis surface: (a) build an incidence relation I from M to the input B-Reps {Si }, (b) use I to determine the elements of M corresponding to each Si . 2. Approximate the organic region with globally smooth editable ‘‘organic’’ surfaces {Oi }: (a) relax the boundary between the organic and incident regions, (b) compute cross- and size-field-guided quad parameterization of the organic region, (c) construct the smooth surfaces {Oi } approximating the organic region. 3. Compose a watertight B-Rep by combining the input solids and the constructed organic surfaces: (a) pull the boundaries of {Oi } in the interior of S, (b) construct contact curves between {Oi } and {Si }, (c) compose a watertight B-Rep from {Oi } and {Si }. 4. Mesh partition The first stage of our algorithm partitions M into a pair of disjoint triangle subsets T = C ⊔ O, where C represents the input solid B-Reps in M, and O represents the boundary surface of the ‘‘organic’’ material synthesized by the solver. We use ⊔ to denote the union of disjoint sets. To compute this partition, our method builds an ‘‘incidence’’ relation I : M → P ({Si }) using the solver geometry representation, where P ({Si }) is the power set of {Si }. I can map a mesh element to zero, one or multiple input solids. 4.1. Incidence data We start by constructing incidence data linking the generative mesh to the input solids. The incidences are stored as attribute(s) for each mesh vertex that is considered ‘‘a representative’’ of an input solid surface. A vertex can be incident with multiple input solids. Notably, our method does not require a representative vertex to lie precisely on the incident input surface(s), and our algorithms are designed to perform robustly when the vertex positions are displaced within the generative solver B-Rep discretization accuracy, e.g., level set grid edge size. The incidence data construction varies on how the solver represents the input solids internally. We propose methods for incidence extraction for two common solver representations, level set and volumetric mesh. We also describe a fallback distancebased algorithm that can be applied to any generative output mesh if an incidence distance is provided. 4.1.1. Level set incidences Level sets are one of the most popular contemporary representations for generative design, used both in research work and in the CAD industry [4]. M is extracted from the 0-isosurface ¯ using a contouring algorithm such as [43–45], where G ¯ is ∂G
Fig. 5. Top: 2D level set incidences for three input solids {A, B, C }: keep-in, seed and keep-out. Organic material added {O} and seed solid material subtracted {E } by the solver as shown on the grid points. Bottom: Grid vertex incidences are accumulated at mesh vertices extracted in each grid cell. Vertices extracted from purely additive organic cells do not have incidence attributes.
the level set representation of the generative design. Each mesh vertex position is computed by combining the corner positions of the level set grid cell containing the vertex. Correspondingly, we construct the mesh vertex incidences as the union of the grid cell corner incidences as follows (Fig. 5): 1. Construct the incidence set I(gk ) = {Si |gk ∈ Si } for each grid point gk using the standard ray intersection count test [46]. For efficiency, only grid points participating in the mesh extraction are tested, i.e., the corner points of the narrow ¯. band of cells intersecting ∂ G ¯ ∧ I(gk ) = {Si } for a seed Si , the volume around gk 2. If gk ∈ /G ¯ , so set I(gk ) = {Ei }, where Ei ⊆ Si has been eroded from G is the eroded subset of Si . 3. For each Vi extracted in a grid cell Gj with corners {gj,k }7k=0 , ⋃7 construct the incidence set I(Vi ) = k=0 I(gj,k ). 4.1.2. Volumetric mesh incidences Volumetric meshes, e.g., tetra- and hexahedral meshes, are another commonly used generative design representation, especially popular in topology optimization applications [47,48]. ¯ is straightIncidence data extraction from the volumetric mesh G ¯: forward since M ≡ ∂ G Let K ⊆ S be the union of the input keep-in and seed solids. ¯ is initialized by piece-wise-linear sampling of K, i.e., G ¯ ≊ K. G
198
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Hence, initially M ≊ ∂ K. During initialization, simply set ∀Vi the incidence attribute I(Vi ) = {Sj } where ∂ Sj is the input solid boundary sampled to obtain the position of Vi . The attributes ¯ , and are recovered persist during the generative optimization of G ¯ from the final M ≡ ∂ G for all vertices that have not been either eroded or added by the solver. 4.1.3. Distance incidences The generative design input to our algorithm may not be based on volumetric optimization, or the solver geometry representation may be unavailable or proprietary. In such cases, we can construct incidence data if we are given a distance d > 0 representing the accuracy of the underlying solver representation. Incidences are constructed as follows: ∀Vi and ∀Sj , compute Dj (Vi ) where Dj (·) is the signed distance function of ∂ Sj , Dj (p) ≤ 0 for p ∈ Sj . Set I(Vi ) = {Sj |Dj (V { i ) ≤ }d}. If I(Vi ) = {Sj } ∧ Dj (Vi ) < −d for a seed Sj , set I(Vi ) = Sj , Ei . For efficiency when |V| and |{Si }| are large, we can precompute a bounding hierarchy of S for proximity searches [49–51]. 4.2. Partitioning We now extend the incidence data to the mesh triangles as follows: If all three vertices of a triangle Tl = {Vi , Vj , Vk } are incident with some solid(s), construct the triangle incidence I(Tl ) = I(Vi ) ∪ I(Vj ) ∪ I(Vk ). Otherwise, set I(Tl ) = ∅: if a triangle has at least one organic vertex, the entire triangle is considered organic, even if another vertex is incident with a solid. From the triangle incidences we construct the mesh subsets corresponding to each input solid Si ↔ Ci = {Tj |Si ∈ I(Tj )} (see inset below). Notably, Ci ∩ Cj ̸ = ∅ if ∂ Si ∩ ∂ Sj ̸ = ∅, and Ci = ∅ if Si is consumed by organic material. We now construct the incident ⋃ C= Ci and organic O = T \ C parts to complete the partition T = C ⊔ O. 4.3. Organic components A pair of triangles Ti and Tj are edge-connected if they share an edge Ek⨆= Ti ∩ Tj . We refer to the elements {Oi } of the smallest partition Oi = O in disjoint, edge-path-connected subsets as the organic components of O. This partition is trivially computed by ‘‘flood-fill’’ through edge connections.
5.1. Boundary vertices Two edges Ei and Ej are connected( if they ) share a vertex Vk = Ei ∩Ej . We refer to Vk as boundary if ∃ Ei , Ej : Ei ∈ B∧Ej ∈ B. Most boundary vertices connect a single pair of boundary edges and are known as 2-manifold. However, occasionally we encounter partitions where a vertex Vk is connecting multiple pairs of boundary edges, i.e., Vk is a non-manifold boundary vertex. Such vertices introduce complications in our downstream processing that must be resolved by special case handling. In practice, it is more convenient to repair such configurations by slightly altering the partitioning as shown on the inset below: for all non-manifold vertices {Vk } we reassign {Tj |Tj ∩ {Vk } ̸= ∅ ∧ Tj ∈ C} to O. {Vk } are no longer boundary, but now ⨆ we have to recompute the organic component partition Oi = O, since the reassigned triangles introduce new edge connections in O which may connect previously disconnected subsets. We iterate this process until all boundary vertices become 2-manifold and typically, it takes only a few iterations to complete. Convergence is ensured for pathological input data as well, since the entire mesh is eventually reassigned to the organic region.
5.2. Boundary curves Let us consider ⨆ the elements {Bi } of the smallest organic boundary partition Bi = B in path-connected, disjoint subsets. We compute {Bi } by tracing the vertex connections using an arbitrary unprocessed boundary edge as a seed. The algorithm completes the current subset Bi when closing the traced edge sequence. This construction shows that {Bi } represent closed, piece-wise-linear curves, and we refer to them as the boundary curves of the partition T = C ⊔ O.
5. Organic boundary
5.3. Boundary relaxation
Let B be the subset of edges separating O and C. We refer to B as the organic boundary and to its elements as the organic boundary edges, or just boundary edges: this shorthand notation is not ambiguous since M is closed by definition.
We further precondition the next stage of our algorithm by relaxing the shape of the boundary curves {Bi } to reduce jaggedness and length differences between adjacent boundary edges. This step is essential to obtain high quality, simpler quad-maps,
Fig. 6. Our boundary relaxation leads to high-quality contact geometry between the organic surfaces and the input solids. Left: input mesh zoomed in on a high-curvature organic boundary segment. Middle: the boundary segment after relaxation. Right: the organic surface and the input solid.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
which in turn lead to improved editability of the output organic surfaces. To avoid shape shrinkage caused by smoothing algorithms such as [52,53], our curve relaxation is restricted to the surface of M and works as follows: 1. Adaptively remesh a small channel of triangles around B to reduce adjacent boundary edge length differences. Our channel simply includes all triangles that are in the tworing neighborhood of the boundary vertices, but alternative definitions are also applicable. We remesh by applying a curvature-adaptive algorithm [54] and refine the adjacent mesh triangles to conform to the changes. 2. Smooth B by applying the surface-restricted curve smoothing algorithm [55], which iteratively decreases angular and length deviations on the boundary curve by shifting the boundary vertex positions towards their averages. While step (1) is not strictly required, we found that it leads to much better relaxation outcomes, mostly because the improved mesh sampling significantly reduces the adjacent boundary edge length differences (Fig. 6). 6. Organic surfaces We construct an editable smooth surface Oi for each organic component Oi . As the construction is identical ∀Oi , in this section the component index i is omitted to simplify the notation. Our representation of choice is T-NURCCs [56] due to the following useful properties:
199
curvature-aligned cross-field computation algorithms, e.g., [62], are also applicable. To estimate curvature tensors on O, we apply the method in [63], but other mesh curvature estimation methods are applicable as well, e.g., [64–66]. The computed cross-field consists of a pair of orthogonal directions ∀Ti ∈ O representing the desired quad edge alignment and contains a small number of singularities assigned to vertices where the field flow bifurcates. Since singularities break the quad edge flow, affecting the editability of O, we add regularity constraints { ⏐ on all boundary vertices } with a ‘‘flat’’ interior surface angle Vi ⏐ π/2 < ∢(Vi ) < 3π/2 (see inset). Due to the boundary relaxation step in Section 5.3, vertices where this condition is violated are very few and far between. We further propagate the boundary curve tangent directions over a small channel along B, replacing the estimated principal curvature directions. This reinforces the boundary structure by aligning the cross-field to the boundary directions within the channel. The channel size is set automatically to match the expected local control mesh edge length. This setup leads to a simple, as regular as possible, quad structure along the control mesh boundary (Fig. 8), allowing improved outcomes when constructing contact curves between O and {Si }.
• The surface is defined by a quad control mesh, a very popular tool for shape editing employed by industry designers.
• Local refinement [57] allows for error reduction by adding control points only where needed [27,28].
• General CAD algorithms, e.g., intersection, offset, can access the geometry as an atlas of C 2 continuous quad patches. • The surface can be easily converted to a G1 continuous NURBS network for export to legacy CAD formats. To obtain a quad partition of the organic shape suitable for defining a meaningfully aligned control mesh, we compute the Integer Grid Map (IGM) φ : O → Ω ⊂ R2 [34,37]. This map has some essential advantages compared to alternative mesh parameterization techniques [58,59]:
• The quad mesh edges can be aligned along the shape principal curvatures by computing a smooth cross-field [32].
• Shapes containing features of vastly varying sizes can be represented efficiently by specifying an adaptive size-field.
• No need to specify a cut network for the high genus organic forms often present in generative designs (Fig. 7). Instead, only a small number of cone singularities are required [31], which we obtain from the cross-field. • Map singularities are mapped to integer points, and the image boundary φ : B → ∂ Ω consists of integer isolines, which allows us to extract a quad partitioning of Ω . 6.1. Cross-field Following the established MIQ method [32,60], we compute a globally smooth cross-field approximating the principal curvature directions on O (away from isotropic regions) to specify the desired edge orientation of the quad mesh. As shown in [61], quad edges sampled along principal curvature directions closely match human-designed meshes, provide optimal approximation, naturally represent symmetries, and allow meaningful shape deformations when used as a tool for surface editing. Alternative
Fig. 7. High genus (88 and 54) generative B-Rep conversions with our method.
200
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Fig. 8. Our boundary relaxation (left) allows regularity constraints almost on the entire organic boundary without causing map distortion. This leads to very regular and easily editable organic surface boundaries (right).
6.2. Size-field Generative shapes are often composed of features of vastly different sizes. To avoid wasteful control mesh complexity on large scale features and improve editability, we specify an anisotropic size-field along the cross-field directions. The desired edge lengths (di,u , di,v ) ∀Ti ∈ O are set from Equation (6) in [61]:
Fig. 9. Top: a generative design, intentionally specified with low sampling accuracy, causes a large visible gap between the boundary of the organic surface and the contact input solid. Middle: our pull algorithm applied on the organic surface boundary. Bottom: the intersection curve is homeomorphic with the boundary which allows us to construct a watertight B-Rep.
√ ( ) 2 di,t (κi,t ) = 2 ϵt − ϵt , |κi,t |
7. Boundary contact
t = {u, v},
where (κi,u , κi,v ) are the estimated curvatures along the crossfield directions on Ti . (ϵu , ϵv ) are automatically chosen to retain the smallest shape features by ensuring ϵt < 2/maxi κi,t . 6.3. Approximation Equipped with cross- and size-fields, we now compute the IGM φ , extract a quad partitioning of Ω and construct the organic surface O. Depending on the application requirements, different strategies can be used to approximate the organic geometry: 1. If simpler, easier to edit organic surfaces are preferable, we compute the control mesh vertex positions by solving the linear min L2 (O, O) problem without bounding the maximum deviation. In this case O can be represented as a Catmull–Clark surface [67] as well. 2. If a tight fit of O is sought, we can solve the non-linear bound problem L∞ (O, O) < δ for a tolerance δ > 0 to output an adaptively refined T-NURCC surface by applying a similar approach to [27]. We set δ to match to the √ generative design accuracy, e.g., for level set designs δ = 3× the voxel grid edge length. 3. As a compromise between surface complexity and accuracy, we can compute an adaptively refined T-NURCC surface which approximates tightly only the organic boundary by solving L∞ (B, ∂ O) < δ for the boundary control points and then min L2 (O, O) for the interior. This ensures that the next step in our algorithm does not change significantly the organic shape near the boundary when constructing contact curves with the incident solids.
⋃ Typically, the constructed organic surface boundary ∂ O = ∂ Oi is within a small distance (≊ the generative solver accuracy) away from the input solid constraints S. For the composed output B-Rep to be watertight, it is necessary that ∂ O ⊂ S. Hence, we apply a small ‘‘pull’’ modification to each organic surface boundary ∂ Oi to move it in the interior of S (Fig. 9). Then, } corresponding B-Rep contact intersection curves {we construct the Ci,j = Oi ∩ ∂ Sj to ensure the successful outcome of our subsequent watertight B-Rep composition. As the algorithm applies identically ∀Ci,j , we omit the index pair i, j. 7.1. Boundary correspondence From the mesh partitioning computed in Section 4.2 we determine the solid S incident with B. If there are multiple adjacent solids (of different types) incident with the same boundary curve, we consider S to be their union. By construction, such solids have overlapping boundaries and therefore their union defines a valid, path-connected solid. Then, from the IGM we find the corresponding organic surface boundary φ : B → ∂ O. 7.2. Boundary pull Let D(·) be the signed distance to ∂ S, D(p) ≤ 0 for p ∈ S. Our algorithm ‘‘pulls’’ ∂ O so that D(∂ O) = maxp∈∂ O D(p) ≤ η, where η < 0 is a small negative safety tolerance adapted to ensure that the intersection C = O ∩ ∂ S and ∂ O are homeomorphic curves (see Section 7.3). D(∂ O) ≤ η < 0 implies ∂ O ⊂ S o , therefore the contact between S and O is watertight. ∑ Using the parameterization c(t) = i ci Ni (t) of ∂ O, where {Ni } are the NURBS basis functions along ∂ O and {ci } are the corresponding boundary control points, we can evaluate D(∂ O)
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
maxp∈∂ O D(p) to {ci }:
=
=
maxt D(c(t)) and ∇ D(∂ O) with respect
( ) ∂ maxt D(c(t)) ∇ D(∂ O) = . . . , ,... . ∂ ci We now apply the gradient descent method to move {ci } along −∇ D(∂ O) until D (∂ O) ≤ η. To efficiently implement the evaluations of D(c(t)), we use a hierarchical space parameterization around S as in Section 4.1.3. We also speed-up the line search maxt D(c(t)) by maintaining a dense uniform sample grid { {t}j } over the domain of c(t) caching the basis function values Ni (tj ) . Loosely approximating the partial derivatives ∂ D(∂ O)/∂ ci with ⃗ i at the projection point qi ∈ ∂ S of ci works the unit normal n quite well and reduces the computation time significantly. If qi happens to lie on a sharp edge or a corner vertex, we define ni as the average of the surface normals incident with qi . The pull process typically converges quickly in a few iterations. Occasionally, we encounter cases where the algorithm fails to pull small segments of the boundary inside S within a certain iteration number (e.g. 100). To correct such segments, we bisect them by local insertion [57] of a new boundary vertex and rerun the pull algorithm. We iterate boundary pull and refinement until convergence: in most observed instances of this issue, only a few new vertices and refinement steps were needed. 7.3. Safety tolerance In practice, due to accuracy limitations of the underlying geometry representations, surface–surface intersection implementations may fail to compute a correct intersection result if the boundary of one of the surfaces is nearly coincident with the other surface. Therefore, if ∂ O is nearly coincident with ∂ S, the computed intersection curve C may not be homeomorphic to ∂ O. To handle such situations, our pull algorithm initially sets η = −1× the accuracy tolerance used by the underlying intersection algorithm. Once the pull is complete, we compute C = O ∩ ∂ S and check if it is homeomorphic with ∂ O. If this is not the case, we ‘‘double up’’ the safety tolerance η → 2η and resume the pull. Applying this strategy iteratively allows our method to work around such accuracy limitations quite robustly. As a performance optimization, it may be preferable to initialize |η| with a larger initial tolerance if it is known that the underlying surface accuracy is limited. 7.4. Tolerant modeling Some B-Rep modeling kernels, e.g., [68], support ‘‘tolerant’’ edges and vertices representing contact curves and points for non- or partially-intersecting surfaces [69]. If tolerant modeling capabilities are available, a distance µ > 0 representing the largest supported tolerance bounding the contact surfaces gap can be specified. We can then relax the safety tolerance η and set it to the positive value η = µ > 0 to pull ∂ O within µ from S. Hence, there is no need to compute the intersection C = O ∩ ∂ S, as we can simply prescribe C = ∂ O = O ∩ Sµ , where Sµ = {p|D(p) ≤ µ}. 7.5. Organic shape change The incidence construction in Section 4.1 implies that B is contained within a distance d > 0 from S, where d ≈ the generative solver accuracy. Consequently, we have D (∂ O) ≤ 2d initially if fitting strategies (2) or (3) in Section 6.3 are applied. Therefore, the shape change introduced by a ‘‘pull’’ of ∼2d is
201
comparable with the design accuracy tolerance and does not affect the validity of the generative result. If approach (1) in Section 6.3 is used instead, the initial D (∂ O) is unbounded. However, in practice we observed that ∂ O is sufficiently close to S for the introduced changes to be acceptable. B-Rep modeling computations are many orders of magnitude more accurate than the generative representation sampling tolerance. Hence, |η| ≪ d and the shape change caused by the additional safety tolerance pull D (∂ O) ≤ η < 0 is negligible. 8. Composition The final stage of our method composes the watertight B-Rep representing the generative design by combining the organic surfaces {Oi } and the input solids {Si }. Contact curves for ∀ {∂ Oi } are now constructed, so it only remains to determine what volume must be assigned to the output generative B-Rep. So far, we did not disambiguate in our algorithms between keep-in and keep-out input solids. We now define the subsets K ∪ L ∪ P = S where K, L and P are respectively the keep-in, keep-out and seed input solids. 8.1. Surface Boolean
S and {Oi } can be combined using existing well-established Boolean operations that function on solids and orientated sheets [70–72] as follows: 1. Combine the ⨆surfaces {Oi } into a single (disjoint) orientated sheet O = Oi . 2. Perform a Boolean unite G = S∪O . Because by construction ∂ O ⊂ S, the operation forms a valid solid G. 3. Subtract the keep-out solids to yield the desired output B-rep G \ L. 8.2. Cellular composition An alternative algorithm, that we apply in our method, uses cellular modeling techniques, which are well established in the art [73,74]. This has several advantages:
• Boolean operations on oriented sheets and partial solids have always proved confusing and counter-intuitive both for users and developers. • The cellular approaches yield an intermediate cell decomposition of the generative design, which can be used to customize the result or diagnose failures due to bad input. 8.2.1. Cell creation Let U ⊃ S ∪ {Oi } be a ‘‘universe’’ void-less solid, e.g., an expanded bounding box or sphere. We partition }U into path{ ⋃ connected solid cells Ui | Ui = U ∧ Uio ∩ Ujo = ∅ bounded by the organic surfaces {Oi }, the solid surfaces {∂ Si }, and the universe surface ∂ U, i.e., {∂ Ui } = {Oi } ∪ {∂ Si } ∪ ∂ U. By construction, there is a single cell U0 containing the universe boundary ∂ U ⊂ ∂ U0 . Discard U0 and the remaining cells satisfy {∂ Ui } = {Oi } ∪ {∂ Si }. 8.2.2. Cell classification We determine now if each Ui is either inside or outside of our final generative B-Rep. If Ui ⊆ Sj for some Sj ⊆ K, Ui is classified as inside, otherwise as outside. { } For all other cells, we inspect the cell boundary ∂ U . Let Fi,j be the i { ⋃ } surface face partition of ∂ Ui , Fi,j | Fi,j = ∂ Ui ∧ Fio,j ∩ Fio,k = ∅ , where each face Fi,j coincides with a path-connected subset of a single surface Ok or ∂ Sk :
202
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Fig. 10. B-Rep conversions of generative designs. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
1. If ∀Fi,j , Fi,j ⊂ {∂ Sk }, Ui is classified outside if it is a void of S, or inside otherwise, since then Ui ⊆ Sj ∈ P. 2. If ∃Fi,j ⊆ Ok , and Ok has the same orientation as Fi,j , then the cell is provisionally classified inside. 3. If ∃Fi,j ⊆ Ok , and Ok has the opposite orientation of Fi,j , then the cell is provisionally classified outside. 4. If both (2) and (3) are true, then the cell is classified inside and a diagnostic indicating an invalid input is reported. 5. Otherwise the provisional classification is adopted.
the input mesh is shown after being partitioned into incident
We compose the output generative solid B-Rep by uniting all inside cells {Ui }.
tive design synthesis which can take many hours. For models
9. Results We include five more generative designs conversions output by our method (Figs. 10 and 11, right column). In all illustrations,
(yellow) and organic (white) regions. Input keep-in (transparent, red tint) and keep-out (transparent, blue tint) bodies are overlaid on the organic T-NURCC surfaces (white). In all results, T-spline patches have been merged where possible to produce a simpler patch layout for the output B-Rep. All conversions shown in the paper are computed in less than 10 min, a small fraction of the time required for the generawith very dense organic regions, the map and approximation computations can become expensive. In such cases, error-bound decimation [75] of the organic mesh away from the organic boundaries can yield much better performance at a negligible accuracy loss.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
203
Fig. 11. B-Rep conversions of generative designs (continued). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
9.1. Robustness
Ω can become quite distorted near concave segments of ∂ O.
Our method is tested extensively and applied in a production environment where hundreds of newly synthesized generative designs are converted to CAD B-Rep models every week. The failures we have observed and investigated so far are caused by one or more of the following factors:
Such maps are difficult to process by the underlying geometry modeling kernel when computing the intersection ∂ S ∩ O and can lead to failed contact curve computations. Reducing the IGM distortion in these cases is an open problem. Map optimization techniques such as [76] offer a promising direction, but are very computationally intensive for complex shapes.
Insufficient solver resolution. Generative design constraint setups occasionally contain very small gaps that are inadvertently left between adjacent input solids, usually of different types, e.g., keep-in or and keep-out. In practice, solvers ignore such gaps and ‘‘unite’’ the solids in their representation. This is especially evident in level set representations where it is often impractical to reduce the voxel grid size to capture tiny gaps. If ∂ O intersects such gaps, our pull algorithm (Section 7.2) cannot converge as the small parts of ∂ O inside the gaps cannot be pulled in So . A useful approach to resolve this issue is to detect tiny gaps in the generative setup and expand the B-Rep solids to fill the gaps. Organic surface map distortion. Due to the limitations of the applied parameterization technique (Section 6), the IGM φ : O →
Organic surface defects. Our pull algorithm can cause small surface self-intersections near high-curvature regions of ∂ O. If a surface self-intersection overlaps with the contact curve C = ∂ S ∩ O, the intersection algorithm can fail to compute C . Surface fairing energy terms [77] can be added to the optimization in Section 7.2 to mitigate this problem, but we have not implemented this as it occurs very infrequently. 9.2. Limitations A limitation of our algorithm is the inability to recognize G1 joints between organic and incident regions. As a post-processing step, the user can edit the organic surface boundary to make a joint G1 , or they can construct a blend feature (available in most CAD packages) to replace an intersection curve. However,
204
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
such manual interaction can be costly for complex models where the organic boundary consists of hundreds of loops. Hence, an automatic detection of such shape characteristics would be very valuable. 10. Conclusion We propose a fully automatic method for generative design to B-Rep conversion capable of processing extremely complex and intricate shapes, implemented and tested over a wide range of data. We hope it can be a driver for growth of the popularity and adoption of generative design models. Engineers and designers looking to incorporate generative designs in their workflows are no longer constrained by the inaccuracy and detail loss of the solver output: instead they can work with an editable, watertight B-Rep model in a CAD system of choice without the additional, often insurmountable cost of manual conversion. Acknowledgments We would like to thank the anonymous reviewers for their helpful feedback, Max Lyon, David Bommes and Leif Kobbelt for their support with integrating IGM, Andreas Bastian, Bryce Heventhal, Karl Willis and Michael Smell for providing generative design test data, and our colleagues in the Autodesk modeling components, simulation, Fusion 360, platform and research groups for supporting the development of our algorithms. This work is funded by Autodesk, Inc, United States. References [1] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optim 1989;1(4):193–202. [2] Autodesk, Fusion 360. [3] MSC, Nastran - multidisciplinary structural analysis. [4] Autodesk Research, Project dreamcatcher. [5] PTC, Frustum generate. [6] Hoppe H, DeRose T, Duchamp T, McDonald J, Stuetzle W. Surface reconstruction from unorganized points. In: Proc. of ACM SIGGRAPH 92; 1992, p. 71–8. [7] Bajaj CL, Bernardini F, Xu G. Reconstructing surfaces and functions on surfaces from unorganized three-dimensional data. Algorithmica (New York) 1997;19(1):243–61. [8] Kruth JP, Kerstens A. Reverse engineering modelling of free-form surfaces from point clouds subject to boundary conditions. J Mater Process Technol 1998;76(1–3):120–7. [9] Benkő P, Martin RR, Várady T. Algorithms for reverse engineering boundary representation models. Comput Aided Des 2001;33(11):839–51. [10] Berger M, Tagliasacchi A, Seversky LM, Alliez P, Guennebaud G, Levine JA, et al. A survey of surface reconstruction from point clouds. Comput Graph Forum 2017;36(1):301–29. [11] Geng Z, Bidanda B. Review of reverse engineering systems - current state of the art. Virtual Phys Prototyp 2017;12(2):161–72. [12] Buonamici F, Carfagni M, Furferi R, Governi L, Lapini A, Volpe Y. Reverse engineering modeling methods and tools: a survey. Comput Aided Des Appl 2018;15(3):443–64. [13] Mangan A, Whitaker R. Partitioning 3D surface meshes using watershed segmentation. IEEE Trans Vis Comput Graphics 1999;5(4):308–21. [14] Garland M, Willmott A, Heckbert PS. Hierarchical face clustering on polygonal surfaces. In: 2001 ACM symposium on interactive 3D graphics; 2001, p. 49–58. [15] Cohen-Steiner D, Alliez P, Desbrun M, Cohen-Steiner D, Alliez P, Desbrun M. Variational shape approximation. In: Proc. of ACM SIGGRAPH 04, vol. 23; 2004, p. 905–14. [16] Wu J, Kobbelt L. Structure recovery via hybrid variational surface approximation. Comput Graph Forum 2005;24(3):277–84. [17] Attene M, Falcidieno B, Spagnuolo M. Hierarchical mesh segmentation based on fitting primitives. Vis Comput 2006;22(3):181–93. [18] Yan D-M, Liu Y, Wang W. Quadric surface extraction by variational shape approximation. In: Proc. of GMP 06; 2006, p. 73–6. [19] Yan D-M, Wang W, Liu Y, Yang Z. Variational mesh segmentation via quadric surface fitting. Comput Aided Des 2012;44(11):1072–82. [20] Eck M, Hoppe H. Automatic reconstruction of B-spline surfaces of arbitrary topological type. In: Proc. of ACM SIGGRAPH 96; 1996, p. 325–34.
[21] Rouhani M, Sappa AD, Boyer E. Implicit B-spline surface reconstruction. IEEE Trans Image Process 2015;24(1):22–32. [22] Hoppe H, DeRose T, Duchamp T, Halstead M, Jin H, McDonald J, et al. Piecewise smooth surface reconstruction. In: Proc. of ACM SIGGRAPH 94; 1994, p. 295–302. [23] Krishnamurthy V, Levoy M. Fitting smooth surfaces to dense polygon meshes. In: Proc. of ACM SIGGRAPH 96; 1996, p. 313–24. [24] Litke N, Levin A, Schröder P. Fitting subdivision surfaces. In: IEEE visualization 2001; 2001, p. 319–24. [25] Ma W, Ma X, Tso S-K, Pan Z. A direct approach for subdivision surface fitting from a dense triangle mesh. Comput Aided Des 2004;36(6):525–36. [26] Marinov M, Kobbelt L. Optimization methods for scattered data approximation with subdivision surfaces. Graph Models 2005;67(5):452–73. [27] Li W-C, Ray N, Lévy B. Automatic and interactive mesh to T-spline conversion. In: Proc. of Eurographics SGP’06; 2006, p. 191–200. [28] Wang Y, Zheng J. Adaptive T-spline surface approximation of triangular meshes. In: Proc. of ICICS’07; 2007, p. 1–5. [29] Wu X, Zheng J, Cai Y, Li H. Variational reconstruction using subdivision surfaces with continuous sharpness control. Comput Vis Media 2017;3(3):217–28. [30] Gu X, Yau S-T. Global conformal surface parameterization. In: Proc. of eurographics SGP’03; 2003, p. 127–37. [31] Ray N, Li WC, Lévy B, Sheffer A, Alliez P. Periodic global parameterization. ACM Trans Graph 2006;25:1460–85. [32] Bommes D, Zimmer H, Kobbelt L. Mixed-integer quadrangulation. In: Proc. of ACM SIGGRAPH 09; 2009, p. 77:1–10. [33] Zhang M, Huang J, Liu X, Bao H. A wave-based anisotropic quadrangulation method. In: Proc. of ACM SIGGRAPH 10; 2010, p. 118:1–8. [34] Bommes D, Campen M, Ebke H-C, Alliez P, Kobbelt L. Integer-grid maps for reliable quad meshing. ACM Trans Graph 2013;32(4):1. [35] Bommes D, Lévy B, Pietroni N, Puppo E, Silva C, Tarini M, et al. Quad-mesh generation and processing: A survey. Comput Graph Forum 2013;32(6):51–76. [36] Myles A, Pietroni N, Zorin D. Robust field-aligned global parametrization. ACM Trans Graph 2014;33(4):1–14. [37] Campen M, Bommes D, Kobbelt L. Quantized global parametrization. ACM Trans Graph 2015;34(6):1–12. [38] 3D Systems, Geomagic design X scan-to-CAD solid model software. [39] Autodesk, PowerShape. [40] LimitState, LimitState:3D. [41] Siemens, Convergent modeling. [42] Lai LM, Yuen MM. Blending of mesh objects to parametric surface. Comput Graph 2015;46(C):283–93. [43] Lorensen WE, Cline HE. Marching cubes: A high resolution 3d surface construction algorithm. In: Proc. of ACM SIGGRAPH 87; 1987, p. 163–9. [44] Kobbelt LP, Botsch M, Schwanecke U, Seidel H-P. Feature sensitive surface extraction from volume data. In: Proc. of ACM SIGGRAPH 01; 2001, p. 57–66. [45] Ju T, Losasso F, Schaefer S, Warren J. Dual contouring of hermite data. In: Proc. of ACM SIGGRAPH 02; 2002, p. 339–46. [46] Chen M, Townsend T. Efficient and consistent algorithms for determining the containment of points in polygons and polyhedra. In: EG 1987-technical papers; 1987. [47] Brackett D, Ashcroft I, Hague R. Topology optimization for additive manufacturing. In: Proceedings of the solid freeform fabrication symposium, vol. 1, S. Austin, TX; 2011, p. 348–62. [48] Chen X, Peng D, Gao S. SVM-based topological optimization of tetrahedral meshes. In: Proceedings of the 21st international meshing roundtable; 2013, p. 211–24. [49] Siemens, Collision detection manager (D-cubed CDM). [50] Ehmann SA. Accurate and fast proximity queries between polyhedra using convex surface decomposition. Comput Graph Forum 2001;20(3). Ming C. Lin.. [51] Alliez P, Tayeb S, Wormser C. 3D fast intersection and distance computation. In: CGAL user and reference manual. 4th ed. 2018. [52] Chaikin GM. An algorithm for high speed curve generation. Comput Graph Image Process 1974;20(3):346–9. [53] Taubin G. Curve and surface smoothing without shrinkage. In: Proc. of ICCV 95; 1995, p. 852–7. [54] Alliez P, Ucelli G, Gotsman C, Attene M. Recent advances in remeshing of surfaces. Technical report, AIM@SHAPE network of excellence, 2005. [55] Schmidt RM. Mesh boundary smoothing. 2013. [56] Sederberg TW, Zheng J, Bakenov A, Nasri A. T-splines and T-NURCCs. In: Proc. of ACM SIGGRAPH 03; 2003, p. 477–84. [57] Sederberg TW, Cardon DL, Finnigan GT, North NS, Zheng J, Lyche T. T-spline simplification and local refinement. In: Proc. of ACM SIGGRAPH 04; 2004, p. 276. [58] Floater MS, Hormann K. Surface parameterization: a tutorial and survey. In: Dodgson NA, Floater MS, Sabin MA, editors. Advances in multiresolution for geometric modelling. 2005, p. 157–86.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205 [59] Sheffer A, Praun E, Rose K. Mesh parameterization methods and their applications. Found Trends Comput Graph Vis 2006;2(2):105–71. [60] Bommes D, Zimmer H, Kobbelt L. System and method for generating quadrangulations. 2010. [61] Alliez P, Cohen-Steiner D, Devillers O, Lévy B, Desbrun M. Anisotropic polygonal remeshing. In: Proc. of ACM SIGGRAPH 03; 2003, p. 485–93. [62] Knöppel F, Crane K, Pinkall U, Schröder P. Globally optimal direction fields. ACM Trans Graph 2013;32(4):1. [63] Cohen-Steiner D, Morvan J-M. Restricted delaunay triangulations and normal cycle. In: 19th annu. ACM sympos. comput. geom.; 2003, p. 237–46. [64] Meyer M, Desbrun M, Schröder P, Barr AH. Discrete differential-geometry operators for triangulated 2-manifolds. 2003, p. 35–57. [65] Cazals F, Pouget M. Estimating differential quantities using polynomial fitting of osculating jets. Comput Aided Geom Design 2005;22(2):121–46. [66] Grinspun E, Gingold Y, Reisman J, Zorin D. Computing discrete shape operators on general meshes. Comput Graph Forum 2006;25(3):547–56. [67] Catmull E, Clark J. Recursively generated B-spline surfaces on arbitrary topological meshes. Comput Aided Des 1978;10(6):350–5. [68] Siemens Parasolid, Parasolid.
205
[69] Jackson DJ. Boundary representation modelling with local tolerances. In: Proc. of ACM SMA 95; 1995, p. 247–54. [70] Nef W. Beiträge zur Theorie der Polyeder : mit Anwendungen in der Computergraphik. 1978. [71] Granados M, Hachenberger P, Hert S, Kettner L, Mehlhorn K, Seel M. Boolean operations on 3D selective Nef complexes: Data structure, algorithms, and implementation. In: Di Battista G, Zwick U, editors. Algorithms - ESA 2003. 2003, p. 654–66. [72] Hachenberger P, Kettner L. 3D boolean operations on nef polyhedra. In: CGAL user and reference manual. 4th ed. 2018. [73] Arnon D. A cellular decomposition algorithm for semialgebraic sets. In: Ng EW, editor. Symbolic and algebraic computation. 1979, p. 301–15. [74] Damiand G. Linear cell complex. In: CGAL user and reference manual. 4th ed. 2018. [75] Garland M, Heckbert P. Surface simplification using quadric error metrics. In: Proc. of ACM SIGGRAPH 97; 1997, p. 209–16. [76] Campen M, Kobbelt L. Quad layout embedding via aligned parameterization. Comput Graph Forum 2014;33(8):69–81. [77] Halstead M, Kass M, DeRose T. Efficient, fair interpolation using Catmull–Clark surfaces. In: Proc. of SIGGRAPH 93; 1993, p. 35–44.
Contents lists available at ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Generative Design Conversion to Editable and Watertight Boundary Representation✩ , ✩✩ ∗
Martin Marinov , Marco Amagliani, Tristan Barback, Jean Flower, Stephen Barley, Suguru Furuta, Peter Charrot, Iain Henley, Nanda Santhanam, G. Thomas Finnigan, Siavash Meshkat, Justin Hallet, Maciej Sapun, Pawel Wolski Autodesk, Inc., United Kingdom, USA, Australia, Poland
article
info
Article history: Received 24 April 2019 Accepted 4 May 2019 Keywords: Generative design Topology optimization Mesh segmentation Mesh parameterization Mesh to B-Rep conversion Solid models
a b s t r a c t We present the first, to our knowledge, fully automatic method for the conversion of a general generative design to a watertight B-Rep composed of the exact residual geometry of the input solids and editable T-NURCC surfaces. The design can be synthesized by any combination of both additive and subtractive in-synthesis processes, and constrained by any number of keep-in and -out solids. Our method requires only the input generative setup solids, rudimentary data from the solve, and a boundary mesh of the optimized design. Leveraging the generative solve data, we augment the boundary mesh with ‘‘incidence’’ attributes linking it to the input solids, and partition it into organic and incident regions. The organic regions are parameterized and approximated with T-NURCCs. Then, their boundaries are ‘‘pulled’’ into the adjacent input solids to construct clean intersection curves. Finally, the organic surfaces are combined with the input solids to compose a watertight solid B-Rep of the generative design. Without such an automatic conversion, users seeking to exploit generative designs in downstream workflows (e.g. setting up assembly constraints, simulation analyses, performing parametric edits) have to manually convert the design to a B-Rep, slowing down the path to manufacturing and increasing the product time to market. As a result of our work, engineers can now quickly obtain useful CAD models of any number of generative designs computed by a variety of algorithms, settings and iterations. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction Generative design is a modern approach to synthesizing design alternatives by exploring the entire design space to achieve objectives such as minimal structural compliance, given a set of functional requirements and constraints based on geometry, material specifications, and manufacturing considerations. A precursor to this technology is topology optimization [1], a niche tool employed by a small number of expert users using specialized simulation software and know-how along with a significant budget. The availability of cloud as an ‘‘infinite computing’’ platform ✩ This paper has been recommended for acceptance by Pierre Alliez, Yong-Jin Liu & Xin Li. ✩✩ One or more of the authors of this paper have disclosed potential or pertinent conflicts of interest, which may include receipt of payment, either direct or indirect, institutional support, or association with an entity in the biomedical field which may be perceived to have potential conflict of interest with this work. For full disclosure statements refer to https://doi.org/10.1016/j. cad.2019.05.016. ∗ Corresponding author. E-mail address: [email protected] (M. Marinov). https://doi.org/10.1016/j.cad.2019.05.016 0010-4485/© 2019 Elsevier Ltd. All rights reserved.
and novel manufacturing methods like additive manufacturing have opened up the possibility of using generative technology as a new approach to design. This capability is starting to become available in CAD packages such as [2], enabling users to leverage generative designs in diverse areas including aerospace, automotive, consumer products and architecture. The generative design geometry constraints are defined as solid B-Reps and are placed in a R3 optimization domain to avoid interference (keep-out) and specify attachment joints (keep-in) with other parts in a CAD assembly. However, many generative solvers, e.g., [3–5], do not operate directly on the exact geometry constraint B-Reps. Instead, B-Reps are sampled and replaced with volumetric representations such as level sets or tetra-/hexahedral meshes, which are significantly more convenient and efficient for physical simulations and material synthesis. Inevitably, this conversion introduces a loss of surface precision due to the limitations of the target representation (Fig. 1). While this loss is negligible in the context of the generative computations, in general it prevents the output of a generative design that accurately represents the input solid surface geometry away from the organic material synthesized by the solver.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
195
B-Rep by leveraging the generative solver input and representation to:
• embed the exact input solid boundary surfaces where the design coincides with the input,
• approximate everywhere else the design boundary with globally smooth, editable ‘‘organic’’ surfaces, and
• join all surfaces to form a generative design output B-Rep. Our method is applicable to generative designs that represent arbitrary combinations of material additions and subtractions, and leverages the solver geometry representation. 2. Related work
Fig. 1. A fragment of a mesh extracted from the level set representation of a generative design. The yellow segment represents the overlaid semi-transparent input solid. The accuracy loss is especially apparent near sharp edges and holes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Subsequently, when integrating the design within the complete model assembly, these surface inaccuracies can cause issues such as part interferences and manufacturing deficiencies (Fig. 3). To resolve these inaccuracies, our fully automatic method (Fig. 2) converts a generative design into an editable, watertight
General B-Rep surface reconstruction from meshes and point clouds is a highly prominent and long-standing research topic with particularly important applications in the reverse engineering of existing physical objects [6–9]. Reviewing this work is beyond the scope of our paper, so we refer to several recent comprehensive surveys: [10–12]. More relevant to our method is research work which focuses on partitioning and approximating a mesh with piece-wise smooth surfaces [13–19] and networks of G1 continuous surface patches [20,21]. However, none of these techniques are known to automatically produce results fulfilling the requirement to represent accurately the residual input solid surfaces in the generative design. Approaches approximating the entire mesh with a globally smooth subdivision surface [22–29] appear at first sight directly applicable to our problem, especially if equipped with a global
Fig. 2. Our generative design conversion in brief. Top: Input solids (left) and 0-isosurface mesh (right), partitioned by our method into incident and organic parts. Bottom: Editable T-NURCC surfaces constructed for the organic part (left) and output watertight B-Rep composed from the input and organic surfaces (right). Note that before composition, T-spline patches are merged where possible to produce the simpler patch layout on the right.
196
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Fig. 3. Left: a rod (in blue), designed to fit through a couple of input solids (in yellow). Middle: a general surface reconstruction violates the generative constraints causing interferences (in red). Right: The rod fits properly when assembled with our generative B-Rep, which reuses the input geometry. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
quad parameterization algorithm such as [30–37]. However, the shortcomings of such unrestricted global approximation schemes for generative design models are easily apparent (Fig. 3). It is tempting to ‘‘fix’’ such defects by applying Boolean operations to remove interferences (subtract keep-out solids) and add missing volume (unite keep-in solids). In practice though, Booleans on nearly overlapping surfaces are very unreliable, which leads to outright failures and B-Reps that are problematic for downstream modeling (Fig. 4).
In the CAD industry, some commercial products, e.g. [38,39], offer capabilities to build a B-Rep model using an input mesh as a geometry reference. However, even for moderately complex shapes, these products require significant manual effort to meaningfully align the smooth surfaces with the processed mesh geometry. For the complex, feature-rich shapes found in most generative designs, such a manually-intensive approach is extremely time-consuming and often cost-prohibitive. Instead, our method is fully automatic and does not rely on user guidance. While specialized solvers such as [40] can process input solid B-Reps directly and output a consistent B-Rep instead of a mesh, the shapes they produce are significantly less varied (e.g. lattices) and only applicable for specific physical problems. In contrast, our approach is applicable to generative solvers optimizing for arbitrary physical criteria, producing a variety of unique and useful shapes. Another alternative to our method is the direct embedding of generative design meshes in the output B-Rep instead of smooth, precise surfaces [41,42]. While this approach can be used for some downstream editing and processing, most modeling operations available in the current generation of B-Rep kernels are not applicable to mesh geometry, thus significantly limiting the usefulness of B-Reps with mesh embedding. 3. Overview 3.1. Input data The primary input to our method consists of the generative solver input solid B-Reps {Si } and output mesh M representing the generative design G. To link M to {Si }, we can also process the intermediate solver representation (level set or volumetric mesh). Alternatively, the solver representation accuracy error, e.g., the level set voxel grid edge length, can be supplied instead. {Si } designate sub-spaces of the optimization domain. We use the standard notation ∂ Si and Sio to denote the solid boundary surface and interior. {Si } can be of the following types:
• Keep-in: Si is completely contained, i.e., Si ⊂ G. • Keep-out: Sio is outside G, ∂ Si can overlap. • Seed: Si is fully or partially inside G.
Fig. 4. Top: An organic surface approximating globally the generative mesh in Fig. 2 violates the keep-in and keep-out constraints as shown in Fig. 3. Bottom: Attempts to repair the violations by Boolean operations fail for the yellow keepin solid and produce poor outcomes everywhere else. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
We unite overlapping solids of the same type, so without loss of generality, Si ∩ Sj = ∅ for every i, j of the same type. We denote the union of all input solids (regardless of type) as S = ⋃ Si . M = (V, E, T) is a closed 2-manifold mesh approximating G, where V = {Vi }, E = {Ei }, T = {Ti } are correspondingly the mesh vertices, edges and triangles. If a non-triangular mesh is provided, we triangulate any quads and general polygons.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
197
3.2. Workflow The main stages (Fig. 2) and algorithm steps of our method are listed in order of execution below: 1. Partition the generative design mesh into two regions: incident and organic. The former corresponds to the input solid surfaces and the latter to the generative synthesis surface: (a) build an incidence relation I from M to the input B-Reps {Si }, (b) use I to determine the elements of M corresponding to each Si . 2. Approximate the organic region with globally smooth editable ‘‘organic’’ surfaces {Oi }: (a) relax the boundary between the organic and incident regions, (b) compute cross- and size-field-guided quad parameterization of the organic region, (c) construct the smooth surfaces {Oi } approximating the organic region. 3. Compose a watertight B-Rep by combining the input solids and the constructed organic surfaces: (a) pull the boundaries of {Oi } in the interior of S, (b) construct contact curves between {Oi } and {Si }, (c) compose a watertight B-Rep from {Oi } and {Si }. 4. Mesh partition The first stage of our algorithm partitions M into a pair of disjoint triangle subsets T = C ⊔ O, where C represents the input solid B-Reps in M, and O represents the boundary surface of the ‘‘organic’’ material synthesized by the solver. We use ⊔ to denote the union of disjoint sets. To compute this partition, our method builds an ‘‘incidence’’ relation I : M → P ({Si }) using the solver geometry representation, where P ({Si }) is the power set of {Si }. I can map a mesh element to zero, one or multiple input solids. 4.1. Incidence data We start by constructing incidence data linking the generative mesh to the input solids. The incidences are stored as attribute(s) for each mesh vertex that is considered ‘‘a representative’’ of an input solid surface. A vertex can be incident with multiple input solids. Notably, our method does not require a representative vertex to lie precisely on the incident input surface(s), and our algorithms are designed to perform robustly when the vertex positions are displaced within the generative solver B-Rep discretization accuracy, e.g., level set grid edge size. The incidence data construction varies on how the solver represents the input solids internally. We propose methods for incidence extraction for two common solver representations, level set and volumetric mesh. We also describe a fallback distancebased algorithm that can be applied to any generative output mesh if an incidence distance is provided. 4.1.1. Level set incidences Level sets are one of the most popular contemporary representations for generative design, used both in research work and in the CAD industry [4]. M is extracted from the 0-isosurface ¯ using a contouring algorithm such as [43–45], where G ¯ is ∂G
Fig. 5. Top: 2D level set incidences for three input solids {A, B, C }: keep-in, seed and keep-out. Organic material added {O} and seed solid material subtracted {E } by the solver as shown on the grid points. Bottom: Grid vertex incidences are accumulated at mesh vertices extracted in each grid cell. Vertices extracted from purely additive organic cells do not have incidence attributes.
the level set representation of the generative design. Each mesh vertex position is computed by combining the corner positions of the level set grid cell containing the vertex. Correspondingly, we construct the mesh vertex incidences as the union of the grid cell corner incidences as follows (Fig. 5): 1. Construct the incidence set I(gk ) = {Si |gk ∈ Si } for each grid point gk using the standard ray intersection count test [46]. For efficiency, only grid points participating in the mesh extraction are tested, i.e., the corner points of the narrow ¯. band of cells intersecting ∂ G ¯ ∧ I(gk ) = {Si } for a seed Si , the volume around gk 2. If gk ∈ /G ¯ , so set I(gk ) = {Ei }, where Ei ⊆ Si has been eroded from G is the eroded subset of Si . 3. For each Vi extracted in a grid cell Gj with corners {gj,k }7k=0 , ⋃7 construct the incidence set I(Vi ) = k=0 I(gj,k ). 4.1.2. Volumetric mesh incidences Volumetric meshes, e.g., tetra- and hexahedral meshes, are another commonly used generative design representation, especially popular in topology optimization applications [47,48]. ¯ is straightIncidence data extraction from the volumetric mesh G ¯: forward since M ≡ ∂ G Let K ⊆ S be the union of the input keep-in and seed solids. ¯ is initialized by piece-wise-linear sampling of K, i.e., G ¯ ≊ K. G
198
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Hence, initially M ≊ ∂ K. During initialization, simply set ∀Vi the incidence attribute I(Vi ) = {Sj } where ∂ Sj is the input solid boundary sampled to obtain the position of Vi . The attributes ¯ , and are recovered persist during the generative optimization of G ¯ from the final M ≡ ∂ G for all vertices that have not been either eroded or added by the solver. 4.1.3. Distance incidences The generative design input to our algorithm may not be based on volumetric optimization, or the solver geometry representation may be unavailable or proprietary. In such cases, we can construct incidence data if we are given a distance d > 0 representing the accuracy of the underlying solver representation. Incidences are constructed as follows: ∀Vi and ∀Sj , compute Dj (Vi ) where Dj (·) is the signed distance function of ∂ Sj , Dj (p) ≤ 0 for p ∈ Sj . Set I(Vi ) = {Sj |Dj (V { i ) ≤ }d}. If I(Vi ) = {Sj } ∧ Dj (Vi ) < −d for a seed Sj , set I(Vi ) = Sj , Ei . For efficiency when |V| and |{Si }| are large, we can precompute a bounding hierarchy of S for proximity searches [49–51]. 4.2. Partitioning We now extend the incidence data to the mesh triangles as follows: If all three vertices of a triangle Tl = {Vi , Vj , Vk } are incident with some solid(s), construct the triangle incidence I(Tl ) = I(Vi ) ∪ I(Vj ) ∪ I(Vk ). Otherwise, set I(Tl ) = ∅: if a triangle has at least one organic vertex, the entire triangle is considered organic, even if another vertex is incident with a solid. From the triangle incidences we construct the mesh subsets corresponding to each input solid Si ↔ Ci = {Tj |Si ∈ I(Tj )} (see inset below). Notably, Ci ∩ Cj ̸ = ∅ if ∂ Si ∩ ∂ Sj ̸ = ∅, and Ci = ∅ if Si is consumed by organic material. We now construct the incident ⋃ C= Ci and organic O = T \ C parts to complete the partition T = C ⊔ O. 4.3. Organic components A pair of triangles Ti and Tj are edge-connected if they share an edge Ek⨆= Ti ∩ Tj . We refer to the elements {Oi } of the smallest partition Oi = O in disjoint, edge-path-connected subsets as the organic components of O. This partition is trivially computed by ‘‘flood-fill’’ through edge connections.
5.1. Boundary vertices Two edges Ei and Ej are connected( if they ) share a vertex Vk = Ei ∩Ej . We refer to Vk as boundary if ∃ Ei , Ej : Ei ∈ B∧Ej ∈ B. Most boundary vertices connect a single pair of boundary edges and are known as 2-manifold. However, occasionally we encounter partitions where a vertex Vk is connecting multiple pairs of boundary edges, i.e., Vk is a non-manifold boundary vertex. Such vertices introduce complications in our downstream processing that must be resolved by special case handling. In practice, it is more convenient to repair such configurations by slightly altering the partitioning as shown on the inset below: for all non-manifold vertices {Vk } we reassign {Tj |Tj ∩ {Vk } ̸= ∅ ∧ Tj ∈ C} to O. {Vk } are no longer boundary, but now ⨆ we have to recompute the organic component partition Oi = O, since the reassigned triangles introduce new edge connections in O which may connect previously disconnected subsets. We iterate this process until all boundary vertices become 2-manifold and typically, it takes only a few iterations to complete. Convergence is ensured for pathological input data as well, since the entire mesh is eventually reassigned to the organic region.
5.2. Boundary curves Let us consider ⨆ the elements {Bi } of the smallest organic boundary partition Bi = B in path-connected, disjoint subsets. We compute {Bi } by tracing the vertex connections using an arbitrary unprocessed boundary edge as a seed. The algorithm completes the current subset Bi when closing the traced edge sequence. This construction shows that {Bi } represent closed, piece-wise-linear curves, and we refer to them as the boundary curves of the partition T = C ⊔ O.
5. Organic boundary
5.3. Boundary relaxation
Let B be the subset of edges separating O and C. We refer to B as the organic boundary and to its elements as the organic boundary edges, or just boundary edges: this shorthand notation is not ambiguous since M is closed by definition.
We further precondition the next stage of our algorithm by relaxing the shape of the boundary curves {Bi } to reduce jaggedness and length differences between adjacent boundary edges. This step is essential to obtain high quality, simpler quad-maps,
Fig. 6. Our boundary relaxation leads to high-quality contact geometry between the organic surfaces and the input solids. Left: input mesh zoomed in on a high-curvature organic boundary segment. Middle: the boundary segment after relaxation. Right: the organic surface and the input solid.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
which in turn lead to improved editability of the output organic surfaces. To avoid shape shrinkage caused by smoothing algorithms such as [52,53], our curve relaxation is restricted to the surface of M and works as follows: 1. Adaptively remesh a small channel of triangles around B to reduce adjacent boundary edge length differences. Our channel simply includes all triangles that are in the tworing neighborhood of the boundary vertices, but alternative definitions are also applicable. We remesh by applying a curvature-adaptive algorithm [54] and refine the adjacent mesh triangles to conform to the changes. 2. Smooth B by applying the surface-restricted curve smoothing algorithm [55], which iteratively decreases angular and length deviations on the boundary curve by shifting the boundary vertex positions towards their averages. While step (1) is not strictly required, we found that it leads to much better relaxation outcomes, mostly because the improved mesh sampling significantly reduces the adjacent boundary edge length differences (Fig. 6). 6. Organic surfaces We construct an editable smooth surface Oi for each organic component Oi . As the construction is identical ∀Oi , in this section the component index i is omitted to simplify the notation. Our representation of choice is T-NURCCs [56] due to the following useful properties:
199
curvature-aligned cross-field computation algorithms, e.g., [62], are also applicable. To estimate curvature tensors on O, we apply the method in [63], but other mesh curvature estimation methods are applicable as well, e.g., [64–66]. The computed cross-field consists of a pair of orthogonal directions ∀Ti ∈ O representing the desired quad edge alignment and contains a small number of singularities assigned to vertices where the field flow bifurcates. Since singularities break the quad edge flow, affecting the editability of O, we add regularity constraints { ⏐ on all boundary vertices } with a ‘‘flat’’ interior surface angle Vi ⏐ π/2 < ∢(Vi ) < 3π/2 (see inset). Due to the boundary relaxation step in Section 5.3, vertices where this condition is violated are very few and far between. We further propagate the boundary curve tangent directions over a small channel along B, replacing the estimated principal curvature directions. This reinforces the boundary structure by aligning the cross-field to the boundary directions within the channel. The channel size is set automatically to match the expected local control mesh edge length. This setup leads to a simple, as regular as possible, quad structure along the control mesh boundary (Fig. 8), allowing improved outcomes when constructing contact curves between O and {Si }.
• The surface is defined by a quad control mesh, a very popular tool for shape editing employed by industry designers.
• Local refinement [57] allows for error reduction by adding control points only where needed [27,28].
• General CAD algorithms, e.g., intersection, offset, can access the geometry as an atlas of C 2 continuous quad patches. • The surface can be easily converted to a G1 continuous NURBS network for export to legacy CAD formats. To obtain a quad partition of the organic shape suitable for defining a meaningfully aligned control mesh, we compute the Integer Grid Map (IGM) φ : O → Ω ⊂ R2 [34,37]. This map has some essential advantages compared to alternative mesh parameterization techniques [58,59]:
• The quad mesh edges can be aligned along the shape principal curvatures by computing a smooth cross-field [32].
• Shapes containing features of vastly varying sizes can be represented efficiently by specifying an adaptive size-field.
• No need to specify a cut network for the high genus organic forms often present in generative designs (Fig. 7). Instead, only a small number of cone singularities are required [31], which we obtain from the cross-field. • Map singularities are mapped to integer points, and the image boundary φ : B → ∂ Ω consists of integer isolines, which allows us to extract a quad partitioning of Ω . 6.1. Cross-field Following the established MIQ method [32,60], we compute a globally smooth cross-field approximating the principal curvature directions on O (away from isotropic regions) to specify the desired edge orientation of the quad mesh. As shown in [61], quad edges sampled along principal curvature directions closely match human-designed meshes, provide optimal approximation, naturally represent symmetries, and allow meaningful shape deformations when used as a tool for surface editing. Alternative
Fig. 7. High genus (88 and 54) generative B-Rep conversions with our method.
200
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Fig. 8. Our boundary relaxation (left) allows regularity constraints almost on the entire organic boundary without causing map distortion. This leads to very regular and easily editable organic surface boundaries (right).
6.2. Size-field Generative shapes are often composed of features of vastly different sizes. To avoid wasteful control mesh complexity on large scale features and improve editability, we specify an anisotropic size-field along the cross-field directions. The desired edge lengths (di,u , di,v ) ∀Ti ∈ O are set from Equation (6) in [61]:
Fig. 9. Top: a generative design, intentionally specified with low sampling accuracy, causes a large visible gap between the boundary of the organic surface and the contact input solid. Middle: our pull algorithm applied on the organic surface boundary. Bottom: the intersection curve is homeomorphic with the boundary which allows us to construct a watertight B-Rep.
√ ( ) 2 di,t (κi,t ) = 2 ϵt − ϵt , |κi,t |
7. Boundary contact
t = {u, v},
where (κi,u , κi,v ) are the estimated curvatures along the crossfield directions on Ti . (ϵu , ϵv ) are automatically chosen to retain the smallest shape features by ensuring ϵt < 2/maxi κi,t . 6.3. Approximation Equipped with cross- and size-fields, we now compute the IGM φ , extract a quad partitioning of Ω and construct the organic surface O. Depending on the application requirements, different strategies can be used to approximate the organic geometry: 1. If simpler, easier to edit organic surfaces are preferable, we compute the control mesh vertex positions by solving the linear min L2 (O, O) problem without bounding the maximum deviation. In this case O can be represented as a Catmull–Clark surface [67] as well. 2. If a tight fit of O is sought, we can solve the non-linear bound problem L∞ (O, O) < δ for a tolerance δ > 0 to output an adaptively refined T-NURCC surface by applying a similar approach to [27]. We set δ to match to the √ generative design accuracy, e.g., for level set designs δ = 3× the voxel grid edge length. 3. As a compromise between surface complexity and accuracy, we can compute an adaptively refined T-NURCC surface which approximates tightly only the organic boundary by solving L∞ (B, ∂ O) < δ for the boundary control points and then min L2 (O, O) for the interior. This ensures that the next step in our algorithm does not change significantly the organic shape near the boundary when constructing contact curves with the incident solids.
⋃ Typically, the constructed organic surface boundary ∂ O = ∂ Oi is within a small distance (≊ the generative solver accuracy) away from the input solid constraints S. For the composed output B-Rep to be watertight, it is necessary that ∂ O ⊂ S. Hence, we apply a small ‘‘pull’’ modification to each organic surface boundary ∂ Oi to move it in the interior of S (Fig. 9). Then, } corresponding B-Rep contact intersection curves {we construct the Ci,j = Oi ∩ ∂ Sj to ensure the successful outcome of our subsequent watertight B-Rep composition. As the algorithm applies identically ∀Ci,j , we omit the index pair i, j. 7.1. Boundary correspondence From the mesh partitioning computed in Section 4.2 we determine the solid S incident with B. If there are multiple adjacent solids (of different types) incident with the same boundary curve, we consider S to be their union. By construction, such solids have overlapping boundaries and therefore their union defines a valid, path-connected solid. Then, from the IGM we find the corresponding organic surface boundary φ : B → ∂ O. 7.2. Boundary pull Let D(·) be the signed distance to ∂ S, D(p) ≤ 0 for p ∈ S. Our algorithm ‘‘pulls’’ ∂ O so that D(∂ O) = maxp∈∂ O D(p) ≤ η, where η < 0 is a small negative safety tolerance adapted to ensure that the intersection C = O ∩ ∂ S and ∂ O are homeomorphic curves (see Section 7.3). D(∂ O) ≤ η < 0 implies ∂ O ⊂ S o , therefore the contact between S and O is watertight. ∑ Using the parameterization c(t) = i ci Ni (t) of ∂ O, where {Ni } are the NURBS basis functions along ∂ O and {ci } are the corresponding boundary control points, we can evaluate D(∂ O)
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
maxp∈∂ O D(p) to {ci }:
=
=
maxt D(c(t)) and ∇ D(∂ O) with respect
( ) ∂ maxt D(c(t)) ∇ D(∂ O) = . . . , ,... . ∂ ci We now apply the gradient descent method to move {ci } along −∇ D(∂ O) until D (∂ O) ≤ η. To efficiently implement the evaluations of D(c(t)), we use a hierarchical space parameterization around S as in Section 4.1.3. We also speed-up the line search maxt D(c(t)) by maintaining a dense uniform sample grid { {t}j } over the domain of c(t) caching the basis function values Ni (tj ) . Loosely approximating the partial derivatives ∂ D(∂ O)/∂ ci with ⃗ i at the projection point qi ∈ ∂ S of ci works the unit normal n quite well and reduces the computation time significantly. If qi happens to lie on a sharp edge or a corner vertex, we define ni as the average of the surface normals incident with qi . The pull process typically converges quickly in a few iterations. Occasionally, we encounter cases where the algorithm fails to pull small segments of the boundary inside S within a certain iteration number (e.g. 100). To correct such segments, we bisect them by local insertion [57] of a new boundary vertex and rerun the pull algorithm. We iterate boundary pull and refinement until convergence: in most observed instances of this issue, only a few new vertices and refinement steps were needed. 7.3. Safety tolerance In practice, due to accuracy limitations of the underlying geometry representations, surface–surface intersection implementations may fail to compute a correct intersection result if the boundary of one of the surfaces is nearly coincident with the other surface. Therefore, if ∂ O is nearly coincident with ∂ S, the computed intersection curve C may not be homeomorphic to ∂ O. To handle such situations, our pull algorithm initially sets η = −1× the accuracy tolerance used by the underlying intersection algorithm. Once the pull is complete, we compute C = O ∩ ∂ S and check if it is homeomorphic with ∂ O. If this is not the case, we ‘‘double up’’ the safety tolerance η → 2η and resume the pull. Applying this strategy iteratively allows our method to work around such accuracy limitations quite robustly. As a performance optimization, it may be preferable to initialize |η| with a larger initial tolerance if it is known that the underlying surface accuracy is limited. 7.4. Tolerant modeling Some B-Rep modeling kernels, e.g., [68], support ‘‘tolerant’’ edges and vertices representing contact curves and points for non- or partially-intersecting surfaces [69]. If tolerant modeling capabilities are available, a distance µ > 0 representing the largest supported tolerance bounding the contact surfaces gap can be specified. We can then relax the safety tolerance η and set it to the positive value η = µ > 0 to pull ∂ O within µ from S. Hence, there is no need to compute the intersection C = O ∩ ∂ S, as we can simply prescribe C = ∂ O = O ∩ Sµ , where Sµ = {p|D(p) ≤ µ}. 7.5. Organic shape change The incidence construction in Section 4.1 implies that B is contained within a distance d > 0 from S, where d ≈ the generative solver accuracy. Consequently, we have D (∂ O) ≤ 2d initially if fitting strategies (2) or (3) in Section 6.3 are applied. Therefore, the shape change introduced by a ‘‘pull’’ of ∼2d is
201
comparable with the design accuracy tolerance and does not affect the validity of the generative result. If approach (1) in Section 6.3 is used instead, the initial D (∂ O) is unbounded. However, in practice we observed that ∂ O is sufficiently close to S for the introduced changes to be acceptable. B-Rep modeling computations are many orders of magnitude more accurate than the generative representation sampling tolerance. Hence, |η| ≪ d and the shape change caused by the additional safety tolerance pull D (∂ O) ≤ η < 0 is negligible. 8. Composition The final stage of our method composes the watertight B-Rep representing the generative design by combining the organic surfaces {Oi } and the input solids {Si }. Contact curves for ∀ {∂ Oi } are now constructed, so it only remains to determine what volume must be assigned to the output generative B-Rep. So far, we did not disambiguate in our algorithms between keep-in and keep-out input solids. We now define the subsets K ∪ L ∪ P = S where K, L and P are respectively the keep-in, keep-out and seed input solids. 8.1. Surface Boolean
S and {Oi } can be combined using existing well-established Boolean operations that function on solids and orientated sheets [70–72] as follows: 1. Combine the ⨆surfaces {Oi } into a single (disjoint) orientated sheet O = Oi . 2. Perform a Boolean unite G = S∪O . Because by construction ∂ O ⊂ S, the operation forms a valid solid G. 3. Subtract the keep-out solids to yield the desired output B-rep G \ L. 8.2. Cellular composition An alternative algorithm, that we apply in our method, uses cellular modeling techniques, which are well established in the art [73,74]. This has several advantages:
• Boolean operations on oriented sheets and partial solids have always proved confusing and counter-intuitive both for users and developers. • The cellular approaches yield an intermediate cell decomposition of the generative design, which can be used to customize the result or diagnose failures due to bad input. 8.2.1. Cell creation Let U ⊃ S ∪ {Oi } be a ‘‘universe’’ void-less solid, e.g., an expanded bounding box or sphere. We partition }U into path{ ⋃ connected solid cells Ui | Ui = U ∧ Uio ∩ Ujo = ∅ bounded by the organic surfaces {Oi }, the solid surfaces {∂ Si }, and the universe surface ∂ U, i.e., {∂ Ui } = {Oi } ∪ {∂ Si } ∪ ∂ U. By construction, there is a single cell U0 containing the universe boundary ∂ U ⊂ ∂ U0 . Discard U0 and the remaining cells satisfy {∂ Ui } = {Oi } ∪ {∂ Si }. 8.2.2. Cell classification We determine now if each Ui is either inside or outside of our final generative B-Rep. If Ui ⊆ Sj for some Sj ⊆ K, Ui is classified as inside, otherwise as outside. { } For all other cells, we inspect the cell boundary ∂ U . Let Fi,j be the i { ⋃ } surface face partition of ∂ Ui , Fi,j | Fi,j = ∂ Ui ∧ Fio,j ∩ Fio,k = ∅ , where each face Fi,j coincides with a path-connected subset of a single surface Ok or ∂ Sk :
202
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
Fig. 10. B-Rep conversions of generative designs. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
1. If ∀Fi,j , Fi,j ⊂ {∂ Sk }, Ui is classified outside if it is a void of S, or inside otherwise, since then Ui ⊆ Sj ∈ P. 2. If ∃Fi,j ⊆ Ok , and Ok has the same orientation as Fi,j , then the cell is provisionally classified inside. 3. If ∃Fi,j ⊆ Ok , and Ok has the opposite orientation of Fi,j , then the cell is provisionally classified outside. 4. If both (2) and (3) are true, then the cell is classified inside and a diagnostic indicating an invalid input is reported. 5. Otherwise the provisional classification is adopted.
the input mesh is shown after being partitioned into incident
We compose the output generative solid B-Rep by uniting all inside cells {Ui }.
tive design synthesis which can take many hours. For models
9. Results We include five more generative designs conversions output by our method (Figs. 10 and 11, right column). In all illustrations,
(yellow) and organic (white) regions. Input keep-in (transparent, red tint) and keep-out (transparent, blue tint) bodies are overlaid on the organic T-NURCC surfaces (white). In all results, T-spline patches have been merged where possible to produce a simpler patch layout for the output B-Rep. All conversions shown in the paper are computed in less than 10 min, a small fraction of the time required for the generawith very dense organic regions, the map and approximation computations can become expensive. In such cases, error-bound decimation [75] of the organic mesh away from the organic boundaries can yield much better performance at a negligible accuracy loss.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
203
Fig. 11. B-Rep conversions of generative designs (continued). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
9.1. Robustness
Ω can become quite distorted near concave segments of ∂ O.
Our method is tested extensively and applied in a production environment where hundreds of newly synthesized generative designs are converted to CAD B-Rep models every week. The failures we have observed and investigated so far are caused by one or more of the following factors:
Such maps are difficult to process by the underlying geometry modeling kernel when computing the intersection ∂ S ∩ O and can lead to failed contact curve computations. Reducing the IGM distortion in these cases is an open problem. Map optimization techniques such as [76] offer a promising direction, but are very computationally intensive for complex shapes.
Insufficient solver resolution. Generative design constraint setups occasionally contain very small gaps that are inadvertently left between adjacent input solids, usually of different types, e.g., keep-in or and keep-out. In practice, solvers ignore such gaps and ‘‘unite’’ the solids in their representation. This is especially evident in level set representations where it is often impractical to reduce the voxel grid size to capture tiny gaps. If ∂ O intersects such gaps, our pull algorithm (Section 7.2) cannot converge as the small parts of ∂ O inside the gaps cannot be pulled in So . A useful approach to resolve this issue is to detect tiny gaps in the generative setup and expand the B-Rep solids to fill the gaps. Organic surface map distortion. Due to the limitations of the applied parameterization technique (Section 6), the IGM φ : O →
Organic surface defects. Our pull algorithm can cause small surface self-intersections near high-curvature regions of ∂ O. If a surface self-intersection overlaps with the contact curve C = ∂ S ∩ O, the intersection algorithm can fail to compute C . Surface fairing energy terms [77] can be added to the optimization in Section 7.2 to mitigate this problem, but we have not implemented this as it occurs very infrequently. 9.2. Limitations A limitation of our algorithm is the inability to recognize G1 joints between organic and incident regions. As a post-processing step, the user can edit the organic surface boundary to make a joint G1 , or they can construct a blend feature (available in most CAD packages) to replace an intersection curve. However,
204
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205
such manual interaction can be costly for complex models where the organic boundary consists of hundreds of loops. Hence, an automatic detection of such shape characteristics would be very valuable. 10. Conclusion We propose a fully automatic method for generative design to B-Rep conversion capable of processing extremely complex and intricate shapes, implemented and tested over a wide range of data. We hope it can be a driver for growth of the popularity and adoption of generative design models. Engineers and designers looking to incorporate generative designs in their workflows are no longer constrained by the inaccuracy and detail loss of the solver output: instead they can work with an editable, watertight B-Rep model in a CAD system of choice without the additional, often insurmountable cost of manual conversion. Acknowledgments We would like to thank the anonymous reviewers for their helpful feedback, Max Lyon, David Bommes and Leif Kobbelt for their support with integrating IGM, Andreas Bastian, Bryce Heventhal, Karl Willis and Michael Smell for providing generative design test data, and our colleagues in the Autodesk modeling components, simulation, Fusion 360, platform and research groups for supporting the development of our algorithms. This work is funded by Autodesk, Inc, United States. References [1] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optim 1989;1(4):193–202. [2] Autodesk, Fusion 360. [3] MSC, Nastran - multidisciplinary structural analysis. [4] Autodesk Research, Project dreamcatcher. [5] PTC, Frustum generate. [6] Hoppe H, DeRose T, Duchamp T, McDonald J, Stuetzle W. Surface reconstruction from unorganized points. In: Proc. of ACM SIGGRAPH 92; 1992, p. 71–8. [7] Bajaj CL, Bernardini F, Xu G. Reconstructing surfaces and functions on surfaces from unorganized three-dimensional data. Algorithmica (New York) 1997;19(1):243–61. [8] Kruth JP, Kerstens A. Reverse engineering modelling of free-form surfaces from point clouds subject to boundary conditions. J Mater Process Technol 1998;76(1–3):120–7. [9] Benkő P, Martin RR, Várady T. Algorithms for reverse engineering boundary representation models. Comput Aided Des 2001;33(11):839–51. [10] Berger M, Tagliasacchi A, Seversky LM, Alliez P, Guennebaud G, Levine JA, et al. A survey of surface reconstruction from point clouds. Comput Graph Forum 2017;36(1):301–29. [11] Geng Z, Bidanda B. Review of reverse engineering systems - current state of the art. Virtual Phys Prototyp 2017;12(2):161–72. [12] Buonamici F, Carfagni M, Furferi R, Governi L, Lapini A, Volpe Y. Reverse engineering modeling methods and tools: a survey. Comput Aided Des Appl 2018;15(3):443–64. [13] Mangan A, Whitaker R. Partitioning 3D surface meshes using watershed segmentation. IEEE Trans Vis Comput Graphics 1999;5(4):308–21. [14] Garland M, Willmott A, Heckbert PS. Hierarchical face clustering on polygonal surfaces. In: 2001 ACM symposium on interactive 3D graphics; 2001, p. 49–58. [15] Cohen-Steiner D, Alliez P, Desbrun M, Cohen-Steiner D, Alliez P, Desbrun M. Variational shape approximation. In: Proc. of ACM SIGGRAPH 04, vol. 23; 2004, p. 905–14. [16] Wu J, Kobbelt L. Structure recovery via hybrid variational surface approximation. Comput Graph Forum 2005;24(3):277–84. [17] Attene M, Falcidieno B, Spagnuolo M. Hierarchical mesh segmentation based on fitting primitives. Vis Comput 2006;22(3):181–93. [18] Yan D-M, Liu Y, Wang W. Quadric surface extraction by variational shape approximation. In: Proc. of GMP 06; 2006, p. 73–6. [19] Yan D-M, Wang W, Liu Y, Yang Z. Variational mesh segmentation via quadric surface fitting. Comput Aided Des 2012;44(11):1072–82. [20] Eck M, Hoppe H. Automatic reconstruction of B-spline surfaces of arbitrary topological type. In: Proc. of ACM SIGGRAPH 96; 1996, p. 325–34.
[21] Rouhani M, Sappa AD, Boyer E. Implicit B-spline surface reconstruction. IEEE Trans Image Process 2015;24(1):22–32. [22] Hoppe H, DeRose T, Duchamp T, Halstead M, Jin H, McDonald J, et al. Piecewise smooth surface reconstruction. In: Proc. of ACM SIGGRAPH 94; 1994, p. 295–302. [23] Krishnamurthy V, Levoy M. Fitting smooth surfaces to dense polygon meshes. In: Proc. of ACM SIGGRAPH 96; 1996, p. 313–24. [24] Litke N, Levin A, Schröder P. Fitting subdivision surfaces. In: IEEE visualization 2001; 2001, p. 319–24. [25] Ma W, Ma X, Tso S-K, Pan Z. A direct approach for subdivision surface fitting from a dense triangle mesh. Comput Aided Des 2004;36(6):525–36. [26] Marinov M, Kobbelt L. Optimization methods for scattered data approximation with subdivision surfaces. Graph Models 2005;67(5):452–73. [27] Li W-C, Ray N, Lévy B. Automatic and interactive mesh to T-spline conversion. In: Proc. of Eurographics SGP’06; 2006, p. 191–200. [28] Wang Y, Zheng J. Adaptive T-spline surface approximation of triangular meshes. In: Proc. of ICICS’07; 2007, p. 1–5. [29] Wu X, Zheng J, Cai Y, Li H. Variational reconstruction using subdivision surfaces with continuous sharpness control. Comput Vis Media 2017;3(3):217–28. [30] Gu X, Yau S-T. Global conformal surface parameterization. In: Proc. of eurographics SGP’03; 2003, p. 127–37. [31] Ray N, Li WC, Lévy B, Sheffer A, Alliez P. Periodic global parameterization. ACM Trans Graph 2006;25:1460–85. [32] Bommes D, Zimmer H, Kobbelt L. Mixed-integer quadrangulation. In: Proc. of ACM SIGGRAPH 09; 2009, p. 77:1–10. [33] Zhang M, Huang J, Liu X, Bao H. A wave-based anisotropic quadrangulation method. In: Proc. of ACM SIGGRAPH 10; 2010, p. 118:1–8. [34] Bommes D, Campen M, Ebke H-C, Alliez P, Kobbelt L. Integer-grid maps for reliable quad meshing. ACM Trans Graph 2013;32(4):1. [35] Bommes D, Lévy B, Pietroni N, Puppo E, Silva C, Tarini M, et al. Quad-mesh generation and processing: A survey. Comput Graph Forum 2013;32(6):51–76. [36] Myles A, Pietroni N, Zorin D. Robust field-aligned global parametrization. ACM Trans Graph 2014;33(4):1–14. [37] Campen M, Bommes D, Kobbelt L. Quantized global parametrization. ACM Trans Graph 2015;34(6):1–12. [38] 3D Systems, Geomagic design X scan-to-CAD solid model software. [39] Autodesk, PowerShape. [40] LimitState, LimitState:3D. [41] Siemens, Convergent modeling. [42] Lai LM, Yuen MM. Blending of mesh objects to parametric surface. Comput Graph 2015;46(C):283–93. [43] Lorensen WE, Cline HE. Marching cubes: A high resolution 3d surface construction algorithm. In: Proc. of ACM SIGGRAPH 87; 1987, p. 163–9. [44] Kobbelt LP, Botsch M, Schwanecke U, Seidel H-P. Feature sensitive surface extraction from volume data. In: Proc. of ACM SIGGRAPH 01; 2001, p. 57–66. [45] Ju T, Losasso F, Schaefer S, Warren J. Dual contouring of hermite data. In: Proc. of ACM SIGGRAPH 02; 2002, p. 339–46. [46] Chen M, Townsend T. Efficient and consistent algorithms for determining the containment of points in polygons and polyhedra. In: EG 1987-technical papers; 1987. [47] Brackett D, Ashcroft I, Hague R. Topology optimization for additive manufacturing. In: Proceedings of the solid freeform fabrication symposium, vol. 1, S. Austin, TX; 2011, p. 348–62. [48] Chen X, Peng D, Gao S. SVM-based topological optimization of tetrahedral meshes. In: Proceedings of the 21st international meshing roundtable; 2013, p. 211–24. [49] Siemens, Collision detection manager (D-cubed CDM). [50] Ehmann SA. Accurate and fast proximity queries between polyhedra using convex surface decomposition. Comput Graph Forum 2001;20(3). Ming C. Lin.. [51] Alliez P, Tayeb S, Wormser C. 3D fast intersection and distance computation. In: CGAL user and reference manual. 4th ed. 2018. [52] Chaikin GM. An algorithm for high speed curve generation. Comput Graph Image Process 1974;20(3):346–9. [53] Taubin G. Curve and surface smoothing without shrinkage. In: Proc. of ICCV 95; 1995, p. 852–7. [54] Alliez P, Ucelli G, Gotsman C, Attene M. Recent advances in remeshing of surfaces. Technical report, AIM@SHAPE network of excellence, 2005. [55] Schmidt RM. Mesh boundary smoothing. 2013. [56] Sederberg TW, Zheng J, Bakenov A, Nasri A. T-splines and T-NURCCs. In: Proc. of ACM SIGGRAPH 03; 2003, p. 477–84. [57] Sederberg TW, Cardon DL, Finnigan GT, North NS, Zheng J, Lyche T. T-spline simplification and local refinement. In: Proc. of ACM SIGGRAPH 04; 2004, p. 276. [58] Floater MS, Hormann K. Surface parameterization: a tutorial and survey. In: Dodgson NA, Floater MS, Sabin MA, editors. Advances in multiresolution for geometric modelling. 2005, p. 157–86.
M. Marinov, M. Amagliani, T. Barback et al. / Computer-Aided Design 115 (2019) 194–205 [59] Sheffer A, Praun E, Rose K. Mesh parameterization methods and their applications. Found Trends Comput Graph Vis 2006;2(2):105–71. [60] Bommes D, Zimmer H, Kobbelt L. System and method for generating quadrangulations. 2010. [61] Alliez P, Cohen-Steiner D, Devillers O, Lévy B, Desbrun M. Anisotropic polygonal remeshing. In: Proc. of ACM SIGGRAPH 03; 2003, p. 485–93. [62] Knöppel F, Crane K, Pinkall U, Schröder P. Globally optimal direction fields. ACM Trans Graph 2013;32(4):1. [63] Cohen-Steiner D, Morvan J-M. Restricted delaunay triangulations and normal cycle. In: 19th annu. ACM sympos. comput. geom.; 2003, p. 237–46. [64] Meyer M, Desbrun M, Schröder P, Barr AH. Discrete differential-geometry operators for triangulated 2-manifolds. 2003, p. 35–57. [65] Cazals F, Pouget M. Estimating differential quantities using polynomial fitting of osculating jets. Comput Aided Geom Design 2005;22(2):121–46. [66] Grinspun E, Gingold Y, Reisman J, Zorin D. Computing discrete shape operators on general meshes. Comput Graph Forum 2006;25(3):547–56. [67] Catmull E, Clark J. Recursively generated B-spline surfaces on arbitrary topological meshes. Comput Aided Des 1978;10(6):350–5. [68] Siemens Parasolid, Parasolid.
205
[69] Jackson DJ. Boundary representation modelling with local tolerances. In: Proc. of ACM SMA 95; 1995, p. 247–54. [70] Nef W. Beiträge zur Theorie der Polyeder : mit Anwendungen in der Computergraphik. 1978. [71] Granados M, Hachenberger P, Hert S, Kettner L, Mehlhorn K, Seel M. Boolean operations on 3D selective Nef complexes: Data structure, algorithms, and implementation. In: Di Battista G, Zwick U, editors. Algorithms - ESA 2003. 2003, p. 654–66. [72] Hachenberger P, Kettner L. 3D boolean operations on nef polyhedra. In: CGAL user and reference manual. 4th ed. 2018. [73] Arnon D. A cellular decomposition algorithm for semialgebraic sets. In: Ng EW, editor. Symbolic and algebraic computation. 1979, p. 301–15. [74] Damiand G. Linear cell complex. In: CGAL user and reference manual. 4th ed. 2018. [75] Garland M, Heckbert P. Surface simplification using quadric error metrics. In: Proc. of ACM SIGGRAPH 97; 1997, p. 209–16. [76] Campen M, Kobbelt L. Quad layout embedding via aligned parameterization. Comput Graph Forum 2014;33(8):69–81. [77] Halstead M, Kass M, DeRose T. Efficient, fair interpolation using Catmull–Clark surfaces. In: Proc. of SIGGRAPH 93; 1993, p. 35–44.