An algorithm for discrete booleans with applications to finite element modeling of complex systems

Finite Elements in Analysis and Design 68 (2013) 10–27

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

An algorithm for discrete booleans with applications to finite element modeling of complex systems B. Kaan Karamete a,n, Saikat Dey b, Eric L. Mestreau a, Romain Aubry a, Felipe A. Bulat-Jara a a b

Sotera Defense Solutions at US Naval Research Laboratory, 1501 Farm Credit Drive, McLean, VA 22102, USA US Naval Research Laboratory, Code 7130, 4555 Overlook Ave SW, Washington, DC 20375, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 May 2012 Received in revised form 14 January 2013 Accepted 16 January 2013 Available online 11 February 2013

In this paper we describe a robust algorithm for three-dimensional boolean operations between boundary representation objects whose geometry is given by discrete (faceted) data. The algorithm presents a new approach for computing the intersection graph which is critical for robustness. It uses elementary computational-geometry operations such as, facet–segment intersection, point containment in simplices and edge recovery in a plane, to produce high-level boolean operations including union, intersection, difference as well as the imprint of the boundary of one object onto another. We also demonstrate the extension and application of the algorithm to mesh-based volumes. We show the robustness and efficacy of our algorithm by employing it to model complex three-dimensional finite element mesh models such as a complete ship where some of the model components are defined in a CAD-based system while others come from legacy mesh-based facetized representations. Use of our algorithm has enabled automation of modeling of very complex configurations reducing the turnaround time for analyses-ready numerical representations from several months to hours or less. & 2013 Elsevier B.V. All rights reserved.

Keywords: Boolean Triangulation Discrete Mesh Sizing sources Remeshing

1. Introduction Engineering analyses and design of any complex system requires one to get an analyzable representation of the system and its environment. This usually involves, at a minimum, the following aspects:

1. A mathematically rigorous representation of the domain geometry usually described as a boundary-representation (BREP) based geometric model which provides an explicit representation of the topological entities that comprise the boundary of an object along with shape, 2. A discrete representation of the domain geometry usually called the mesh model, and 3. Information needed for the analysis/design that are not in the geometric or the mesh model usually referred to as attributes which typically include material properties and the like. A geometry-centric approach where the mesh model and analysis attributes are associated with, and dependent on the geometry model, offers several advantages that include:

n

Corresponding author. Tel.: þ1 5182696226. E-mail address: [email protected] (B. Kaan Karamete).

0168-874X/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.finel.2013.01.003

1. Robust mesh generation and adaptation due to mathematically rigorous definitions of mesh validity based on concepts of consistency of the mesh model with respect to the geometric model. If the mesh model is not dependent on the geometry model, robust definitions of mesh validity lack and ensuring a valid mesh becomes more expensive. 2. Automated ability to update the information needed by the solvers when geometry is changed. This is not possible in approaches that are mesh-centric where analysis information is attached to the mesh model without any reference to the geometry model associativity. The most productive and straightforward approach to realize a geometry-based analysis and design setup is to use a geometry system1 that provides accurate boundary-representation-based geometry as well as fully automated algorithms for boolean operations involving complex three-dimensional shapes. This provides for the ability to construct the geometric representation of very complex systems (such as an aircraft or a ship) by combining bottom-up and top-down approaches. The geometry-based approach is tailor-made for workflows where geometry may be created from scratch or sourced entirely from a CAD-based system. However, in many cases the description of geometry of the entire system to be analyzed may not exist 1 Such as a computer-aided-design (CAD) system or geometry kernels that underly a CAD system.

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

11

in a CAD system or is in a form that makes it very tedious and time consuming to convert it to a CAD system. One example of such a case is when a new system, like a ship, is to be analyzed for operational effectiveness where the main geometry of the hull and rest of the super structure may come from a CAD-based representation, but many of its component systems, such as a pump, may be reused from older designs for which the only available representation is a simple discrete description of the geometry in the form of a mesh or facet-based data. Another example of interest to us involves modeling urban environments where the ground geometry may be in discrete data coming from a LIDAR scan which needs to be combined with models of vehicles or buildings that may be coming from CAD-based descriptions. In such situations, one can still make the geometry-centric approach work by producing a boundaryrepresentation geometry for those objects that are available only in discrete data form. This requires one to derive a topological representation of the boundary of the object such that the given discrete data represents a consistent discretization of the derived boundary description. It also requires one to assign proper shape to the resulting boundary representation. One can achieve this by either, 1. Reconstructing the geometry of the discretely defined object entirely in a CAD-system or kernel. 2. Reconstructing only the topology of boundary representation and utilizing the shape information derived from the discrete (mesh) data. One way to realize the first approach is to use an interpolation to define the shape of the new object in the geometry system. This incurs errors in the shape and is hard to automate for completely unstructured or irregular data. An alternative is to simply take each mesh entity in the mesh-based description and recreate it as-is in a CAD-based system [1]. This can be entirely automated, but has two main issues; first, it can end up creating too many entities in the geometric model and second, it can produce many small topological features on the boundary of the model that then degrade the quality of the resulting mesh. For these reasons, the second approach is an attractive alternative from the standpoint of providing quick turnaround time for users when faced with complex model data in discrete form as well as fitting into the geometry-(BREP)-based approach in a seamless manner as depicted in Fig. 1. One challenge for the second approach is a robust ability to do boolean operations on the boundary representation objects described with discrete (mesh) data. The format of the input surfaces is either triangular or quadrilateral mesh faces. Our main focus in this paper is implementing a robust algorithm to perform discrete boolean operations between two BREP objects. In particular, the following boolean operations are considered: Union

Given two BREPS, A and B, the result of a union is the BREP C ¼ A [ B that contains the set of points that belong to the closure2 of either A or B. Our implementation supports, both the regularized (manifold) and non-regularized (non-manifold) cases as depicted in Fig. 2(b, c). Difference Given two BREPS, A and B, the result of a difference is the BREP C ¼ A\B that contains the set of points that belong to the closure A but not interior to B as depicted in Fig. 2(d). Intersection Given two BREPS, A and B, the result of an intersection is the BREP C ¼ A \ B that contains the set of points that belong to the closure of both A and B as depicted in Fig. 2(e).

2 Closure denotes the set of points strictly in the interior as well as on the boundary of an object.

Fig. 1. Data flow models depicting the CAD-based and discrete-data based pipelines. The downstream workflows/applications and algorithms work seamlessly irrespective of the data source.

In addition, we also implement the imprint operation which modifies only the boundary topology of one BREP for those portions where it intersects with the boundary of another BREP as depicted in Fig. 2(f). For a more complete description of pointset boolean operations used in geometric modeling, one may refer to [2]. We also show examples employing these algorithms to realize complex modeling workflows including surface and volume meshing of the resulting model, as well as, management of any attribute(s) assigned to the original objects. The rest of the paper is organized as follows: Section 2 defines the topological hierarchy that describes our boundary representation scheme as well as the automated building of boundary topology of the input discrete geometry objects (depicted by the box (A) in Fig. 1). Section 3 gives an outline of the algorithm including the main steps involved; computing the intersection of the boundary of the input objects is described in Section 4; Section 5 presents the determination and topological recovery of the coherent boundary triangulations on the input objects that account for the boundary intersections; Section 6 describes the actual evaluation of the discrete-boolean operations utilizing the boundary intersection information. Section 7 describes the remeshing of the resulting boundary representation based on local mesh modification operations starting from the input discretization. In addition, it describes updating the boundary representation based on the final mesh to reflect a consistent association of the final mesh with the resulting boundary representation of the geometry. It also handles the propagation and update of attribute information associated with the original input objects. Section 8 describes the ability to handle input descriptions that include volume entities and mixed dimension discretizations. Section 9 gives results of the application of the presented algorithm to several examples involving finite element modeling of complex systems and Section 10 concludes with a summary of future work.

2. Preliminaries 2.1. Boundary representation building Since most input discrete data do not originate from a geometry model with a consistent representation of its boundary-topology, we process the discrete geometry input to construct a boundary-representation of a geometric model that corresponds

12

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

Fig. 2. Point-set boolean operations: (a) input BREPS with results for (b) regularized union, (c) non-regularized union, (d) difference, (e) intersection and (f) imprint. Note that the result in (c) has three regions and is non-manifold when used with three-dimensional solid objects.

to the given discrete data. This process is called the discrete geometry building (DGB) and relies on topological boundary feature extraction based on angle thresholds. The process of building the boundary representation of a geometric model for the input data also creates the classification of the mesh entities (vertices, edges, faces and regions) on appropriate entities of the resulting geometric model. For example, a mesh of a cube would result in eight geometrical vertices, twelve geometrical edges and six geometrical faces. This is more than just surface walking and grouping mesh entities because the complete onelevel-up-and-down inter-entity adjacencies are also computed. While not the focus of this paper, this is critical in two ways. First, it enables us to have a consistent definition of validity when performing local mesh modification operations [3,4], and second, it enables us to handle attribute specification based on a boundary representation, a key part of our technical approach. Our approach relies on using the topological data model described in [5,6] and is depicted in Fig. 4. 2.2. Intersection graphs The discrete boolean operations rely on a critical topological concept called the intersection-graph [7]. It is a graph consisting of a chain of edges connected by shared vertices. It is a discrete representation of the intersection-set of the boundary of the two input objects. Fig. 3 depicts intersection graph examples; note that in some special cases the graph may have just a vertex and no (null) edges. The main algorithm for discrete booleans will be explained based on a simple case consisting of two intersecting cylindrical objects as depicted in Fig. 5. The discrete data (mesh) for the two objects is used by the discrete geometry builder algorithm to create the topology of the boundary representation for each object. The resulting boundary representations along with its support mesh entities are illustrated in Fig. 6. The statistics of the boundary topology for the boundary representations are listed in Fig. 6(b, c); the number of geometric model entities created are dependent on the angle thresholds used and, in general, the number of entities in the geometric model is less than those in the input discrete representation. The main steps of the algorithm are discussed next.

Fig. 3. The intersection graph for facets from two BREPS A (thin solid lines) and B (broken lines) is defined by two vertex-connected edge-chains fv1 ,v2 ,v3 g (thick lines) and fv4 ,v4 g (no edge).

common with respect to each other so that various boolean operations could be performed by using simple point-set logic. This requires computation of the intersection graph between the boundaries of the two objects. One key requirement is that the intersection graph must be unique and identical for both input objects. This is needed for two reasons. First, it ensures the proper detection of common volumes between the input BREPs; second, it enables the merging of the input (intersecting) BREPs such that the result is topologically consistent. This is because edges on both BREPs corresponding to the same end vertex intersection locations can be merged as long as the intersection graph is uniquely matching. Merging of edge pairs results in one connected, conformal BREP out of the two intersecting BREPs. Therefore, the first goal is to design an efficient algorithm to create matching intersection graphs on both BREPs. To this end, there are two major approaches:

3. Main steps

a. Explicit creation of intersection-graph ring-topologies (incremental approach) [8,9], b. Implicit creation of intersection-graph ring-topologies (global approach) [10,11].

Discrete boolean operations require the classification of portions of boundaries of the two input objects as inside, outside and

Forming matching edge-pairs requires edge and face splits on each input BREP. Use of one of the two approaches above is

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

13

Fig. 4. Data Model: Geometric entities have one-level-up-and-down adjacencies. A geometrical model may be composed of a set of BREPs. A number of mesh models can be associated with the same geometry model. Mesh entities can be defined on geometrical entities based on the classification rule, e.g., a mesh face can be classified on a geometry face or a geometry region but not on a geometry edge. It is possible to create many different mesh views based on mesh inter-entity adjacency relations. Moreover, partial adjacencies are also possible by the use of scopes, e.g., the definition of a mesh face from a set of mesh vertices classified on a geometrical face will also implicitly create mesh edges if mesh face to mesh edge adjacencies are also requested if the scope of mesh edges are also set at least at GFace level; if the scope of mesh edges are set at GEdge level, then partial adjacency from mesh faces to mesh edges are created, i.e., only edges classified on geometric edges will appear in mesh faces’ mesh edge inter-adjacency list.

Fig. 5. The input is two inter penetrating BREPS. (a) The transparent top view, (b) Hidden line perspective view of the two cylinder BREPs.

required to compute the intersection graph where the splits are to be made. The approach in [8,9] favors the deduction of the intersection graph chains in an explicit manner in which the intersections are performed starting from a seed facet–segment in a walking pattern by visiting the facets and segments always in the immediate neighborhood of previously computed intersection set. This is an incremental process that eventually forms the connected chain(s). The order in which the facets and segments (walk) are visited is also the order in which the split operations are performed. This ensures that the splits create matching edgepair chains for the subsequent edge merging operation. However,

as the authors state, the explicit formation of the intersection graph rings is prone to algorithmic and numerical issues, particularly, for when the mesh/facet sizes for the input BREPs differ significantly. A representative case is illustrated in Fig. 7 when one large mesh face/facet intersects many facets with the other BREP and the resulting intersection points may be part(s) of two or more rings (chains) of the to-be-extracted intersection graph. In the incremental approach, one has to be very careful to keep the parallel-history of walking on each of these potential chain-paths to decipher the connected rings of the common intersection graph. We have developed a new and novel implicit approach to computing the matching intersection graphs that addressed limitations with the explicit approach and is more robust. Implicit formation of the intersection graph requires compiling of relational information between intersecting facets of the two BREPs using unordered facet–segment intersections and localizations. The face–face intersections are computed similar to [10] and [11]. Specific steps of the overall discrete boolean algorithm, using the implicit method for matching intersection graph, are itemized below.

1. Find intersections between BREP1 edges and BREP2 faces. 2. Find intersections between BREP2 edges and BREP1 faces. 3. Localize intersection points on intersected edges and faces.

14

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

4. Create a relational table between faces of BREP1 and BREP2 based on intersection data. 5. Resolve the relational table for each face to create sets of intersection point pairs (intersection graph edges). 6. Triangulate each face in such a way that the entire intersection graph edge set per face is preserved.

BREP1 Face2 BREP2 Face5 ge5 Ed

7. Make sure that all intersection graph face edges are recovered. 8. Edges on both sides must be conformal, create an edge and vertex map on both sides due to (7). This is the conformal intersection graph of BREP1 and BREP2. 9. Categorize faces on BREP1 that are outside and inside of BREP2. 10. Categorize faces on BREP2 that are outside and inside of BREP1. 11. Remove respective categorized faces based on the specified boolean operation. 12. Merge edges of BREP1 with the edges of BREP2 specified in the edge map. Merge Brep1 with Brep2. 13. Remove all geometrical classifications. 14. Rebuild the discrete geometric topology based on angle thresholds, i.e., re-classify the mesh. 15. Transfer of user attributes to new geometrical topologies. 16. Locally remesh to improve quality around the intersection graph. Steps 1–5 will be explained in Section 4, steps 6–8 in Section 5, steps 9–12 in Section 6, and steps 13–16 in Section 7.

4. Computing intersection between discrete boundary representations To create the intersection graph between two input objects, BREP1 and BREP2, it is necessary to compute the facet–segment intersections between

Fig. 6. Discrete geometry builder (DGB) builds the association between geometric topologies and mesh entities. (a) Discrete geometry faces, namely, Face 2 of BREP1 and Face 5 of BREP2 are intersecting. (b) The discrete geometry topological tree with one-level-up-and-down inter-topology adjacencies are shown. Note that Edge 5 is used twice by Face 5 as it is a seam. (c) The number of different level mesh entities are classified against the discrete topology, e.g., 4 mesh vertices are associated with 4 geometric vertices, a total of 92 mesh edges are classified on 6 geometrical edges, and 924 mesh faces are classified on 6 geometrical faces, 3 from each BREP.

i. the edges of BREP1 and the faces of BREP2, and ii. the edges of BREP2 and the faces of BREP1, consecutively as seen in Fig. 10. To compute above intersections efficiently, we build two special octrees [12] using an adaptive subdivision scheme based on the number of facets stored inside each terminal octant. Firstly, we detect intersected geometry faces between the two

Fig. 7. The conceptual layout of two intersecting facet sets on two penetrating BREPs. Robustness issues with creating matching intersection graph using the explicit approach: (a) The parallel formation of two distinct chains, that are tied to the same facet on one BREP, requires complex algorithm to keep the intersection history execution simultaneously, (b) the ambiguity to distinguish between one or two chain formations due to the same facet being involved in intersections, (c) the split on one side of the facet being intersected will result in the modification of the intersection set on the other side which would require the upper chain to be recomputed.

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

15

BREPs using fast bounding-box intersections. Note that geometry faces are clustered set of mesh faces based on some geometrical predicates such as facet normal deviations and analysis attributes. Secondly, the facets (mesh faces) classified on intersected geometry faces of BREP1 are inserted into the first octree. A second octree is generated identically by repeating the same insertion process for BREP2 faces. These two octrees store potentially intersecting sets of mesh faces based on bounding boxes. The intersection of an edge of BREP1 against the faces of BREP2 in item (i) above, is computed in the following manner:

a. Identify potentially intersecting BREP2 faces using the bounding box of BREP1 edge on the tree generated using BREP2 faces, then b. Compute true (‘‘go-through’’) intersection(s) using facet–segment intersections for the faces identified by the boundingbox test. Facet–segment intersection algorithm has been studied by many researchers [8,9,13–15] and briefly summarized here in Fig. 8. Having found the intersection point between the plane of the facet and the edge segment, it is necessary to check if the point of intersection is contained within the bounded triangular face entity. A barycentric face containment check is used for this purpose where the barycentric areas are computed form each edge to the point of intersection. If any of these sub-barycentric areas is less than zero then the intersection point is outside of the face. Otherwise, the intersection point corresponding to the facet entity and the sub-entity where the intersection point lies on the other BREP is localized and recorded (see Fig. 9). All possible facet–segment intersections are computed between BREP1 faces and BREP2 edges and vice versa as depicted in items (i) and (ii) above. Constructing this raw intersection record between the two BREPs is illustrated in Fig. 10. This raw intersection record is not readily useful since splits at these locations would subsequently destroy the entities in the records. Additionally, splits independently performed on each BREP do not guarantee that the results will have matching conformal edge chains between intersection point pairs. This is because there is no implicit information on how these intersection points connect

Fig. 9. Entity localization and containment check for the intersection point p inside the face f for splits. First the distance from p to corner vertices is checked against the tolerance E to decide if p is localized on one of the face vertices vi. Then the distance from p to each edge ei is computed using the barycentric areas Ai. If Ai is less than zero then p is outside of f. If only one Ai is less than the edge relative E tolerance then p is localized on ei. Otherwise, p is deemed to be localized on f.

with each other in this raw form. Therefore, a different form of data structure is required such that the splits would not be destructive to the information in the records and reflect how the intersection points connect with each other to create matching conformal edge chains on both BREPs. To this end, a relational data structure based on faces on both BREPs is constructed. This is described in Section 4.1. 4.1. Face-to-face relational multi-map generation In the raw intersection data, each intersection point belonging to a BREP is associated with a face of the same BREP and a localized mesh sub-entity of the other BREP (Fig. 11a, b). The subentity field of this raw intersection record is replaced by as many records as there are adjacent faces to the localized sub-entity. The inter-entity adjacencies are used to multiply the data records so that all possible face to face relations over each intersection point could be generated as depicted in Fig. 11(c). Next step involves extraction of intersection points that share the same face pairs from the modified records. This is achieved by inverting the table in Fig. 11(c). The result is a multi-keyed look-up table (multimap) from paired faces to each intersection point and depicted in Fig. 11(d). 4.2. Resolving of relational face maps

Fig. 8. Computing the ‘‘go-through’’ intersection location p between a triangular plane v0 ,v1 ,v2 and a segment p0 ,p1 . If the plane–segment intersection condition is satisfied, i.e., the segment end vertices lie on opposite sides of the plane, the scalar parameter t is computed by inserting the line equation into the plane equation. The coordinates of intersection point p can then be computed via the line equation.

A more useful form to deduce intersection graph edges is to relate a set of intersection point pairs within each face. To this end, intersection point pairs are composed from identical face pairs appearing in the multi-map. For each face of the face pair in the multi-map, the algorithm searches through the multi-map to pick points that have the same second face. For instance, the search for face F3 in Fig. 11(d), starts searching from the first record to extract points that have the identical second face to F3 through the depth of multi-map. This search algorithm first finds

16

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

Fig. 10. Illustration of constructing the intersection record for booleans on a small portion of the inter-penetrating B-reps, BREP1 and BREP2. (a) Two intersecting facets of BREP1 and BREP2. (b) Edges of BREP2, E52 and E50 intersecting with the facets of BREP1, F1 and F2, respectively. (c) Edges of BREP1, E67 intersecting with the facet of BREP2, F3. (d) Intersection vertices of BREP1, 1, 2, 3 and BREP2, 10 , 20 , 30 . (e) The raw intersection record on BREP1 and BREP2 are tabulated as the pairs of intersection points and mesh entities that intersections occur on both BREPs.

Fig. 11. Face intersection map (FIM) generation algorithm: (a) The facet–segment intersections between the two BREPs. (b) The result of intersections: a set of facet entity on BREP1, intersection sub-entity on BREP2 for intersection points 1, 2 and facet entity on BREP2 with intersected sub-entity on BREP1 for point 3. (c) Intersected subentities are related to the faces using entity-to-face adjacency, e.g., E67 is adjacent to F1 and F2. (d) face-to-face pairwise relational intersection multi-map is generated; for each intersection point, corresponding face pairs are extracted using the inverse of previous form. (e) From the multi-map, a set of intersection point pairs is deduced for each face. (f) This FIM constitutes the intersection graph edge data for each mesh face with coinciding intersection nodes between 110 ,220 and 330 on BREP1 and BREP2, respectively.

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

two intersection points, (1, 3), corresponding to F1 as identical face pair component. Similarly, the search further continues to extract another two points, (3, 2), corresponding to F2 as identical face pair to F3. The resulting form is simply called Face Intersection Map (FIM) and depicted in Fig. 11(e). The visual representation of FIM is illustrated on two BREP1 faces F1 and F2 and one BREP2 face F3 as a set of connected edges in Fig. 11(f). 4.3. Degeneracy and singularity resolution The face–segment intersection tests, whose equations are pictorially given in Fig. 8 and 9 may result in a situation that the segment is just touching on the face or going through one of its vertices. Our goal is to create matching entities on both BREPs for these cases. The latter one is addressed trivially that the entity at which the intersection happens is found to be one of the facet vertices by the containment checks, and inserted into FIM accordingly. The former case, in which the segment touches on the face, results in one vertex that is not connected to any intersection graph edge. The intersection graph in this case is the vertex itself. This is reflected and tagged in FIM specially. Consequently, the split of the face and the split of the intersecting segment edge at the vertex location recovers the intersection graph vertex on both BREPs identically. In this section, the facet–segment intersections are discussed for ‘‘go-through’’ cases. However, if the facet and the segment lie on the same plane, i.e., the plane of the facet, then extra measures have to be taken. We have experimented the resolution proposed by Aftosmis [10] by using ‘‘virtual perturbation algorithms’’, but concluded to address the issue more concretely as an improvement to our algorithm in a future article.

5. Triangulating faces Face-intersection-map (FIM) is constructed per face, i.e., it is the imprint of the one BREP onto the other one over each face. The task at hand is to insert the vertex-pair sets of FIM into each face and ensure the edges per vertex pair are recovered after the splits.

1

This is the piecewise definition of intersection graph between the two BREPs. The special case of only one vertex touching a mesh face requires the recovery of the vertex which is trivially satisfied by the split operation. For all other cases, the insertion of the two vertices of the vertex-pair does not guarantee that there would be an edge between the vertices as depicted in Fig. 13. The strategy of sorting the vertex-pairs such that the vertex pairs forms an ordered chain increases the chances of having the edges between each vertex pair after the splits. The split of vertices can happen on an existing edge or a face. The original face in which the intersection map vertex pairs are tabulated can be defined as the master face to indicate that the following vertex splits and edge recovery manipulations would be contained within the closure of this master face. The split operation of the very first vertex in the chain requires the master face to be deleted from the mesh data structure. However, retaining this face has algorithmic merits in limiting the operations within the convex boundary of the master face until all vertex splits are performed. The number of new faces after each split operation will not be many. Therefore, the decision of finding if an edge or face split has to be performed at the vertex location are carried by containment checks on all of the newly created faces inside the master face. Again, barycentric face containment check is used to localize the entity for splits. First, the distance from the intersection point, i.e., the vertex coordinate of one the vertex-pairs to each face’s corner vertex is checked against a tolerance to decide if vertex is localized on one of the face vertices. Then the distance from the vertex to each face-edge is computed using the barycentric areas. If any of these areas is less than zero then the vertex is outside of the face which is generally not the situation for the master face due to the nature of the FIM construction process. However, it may happen after the first split as there will be more faces due to the prior splits inside the master face. In such a case, we move to the next adjacent face of the edge whose area coordinate is negative, and start localizing the vertex within the next adjacent face [16]. If only one barycentric area is less than the edge relative tolerance, then the vertex is localized on the corresponding face edge. Otherwise, the vertex is deemed to be localized on the face topology itself. The implementation details of the containment check and

F3

F2

2

3

17

3

3’ 2’ 1’

F1

1

1

3

2

2

3 3

1’

3’

2’

Fig. 12. Triangulate each face independently in FIM splitting the face at each of the intersection point pairs: (a) The intersection graph edges on BREP1 faces F1 and F2 and BREP2 face F3. (b) Face and edge splits are performed within F1, F2 and F3 on nodes 110 ,220 and 330 . (c) Splits should create matching edges on BREP1 and BREP2, this is ensured by edge recovery within each face based on FIM.

18

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

localization of the split entity is very important and has huge impact on the robustness of the entire algorithm as this operation is computationally inexact. The localization algorithm for splits and face containment check are summarized and illustrated in Fig. 9. Having localized the split entity for the vertex in the set of vertex pairs which are already sorted in a loop sense, the next procedure is to actually split the entity at the vertex coordinate and update the set of faces due to the split operation for the repeat of the localization process to the next vertex in the looping pairs (see Fig. 12). 5.1. Edge recovery It is possible after the splits that some edges for inserted vertex pairs are missing. The ordering of the splits lowers the likelihood of missing edges. However, there is no simple split order that will guarantee having the edges between each vertex pair. Therefore, a more sophisticated algorithm incorporating explicit edge-recovery is needed to ensure an edge exist between each vertex pair in FIM. The edge recovery procedure on planar mesh faces is a guaranteed operation using random edge swaps [17]. The edge recovery algorithm is explained extensively in [18]. The recovery of missing edges is not a limiting issue due to the fact that all splits happen on the planar geometry of the original master face. A possible scenario in which there are two edges

missing is shown in Fig. 13. The edge recovery algorithm is applied after the splits for each face appearing in FIM.

5.2. Creation of conformal edge vertex pairs There is a corresponding vertex for each intersection point on BREP1 and BREP2. The edges recovered by splits on the face boundaries are simply the pairwise connections between these vertices specified in FIM. The algorithm forming face-to-face relations in terms of vertex pairs in FIM ensures having matching edges. Therefore, it is a straightforward process to match edges on both sides by using the vertex look-up table on the intersection records. The merging of two BREPs is possible at this stage due to the fact that both BREPs have conformal one-to-one matching edges along their intersection graph. These two matching intersection graphs or imprints of BREP1 onto BREP2 and vice versa is shown in Fig. 14. Boolean operations may require retention and/or deletion of certain portions of BREP boundaries. The BREPs are kept separate to make categorization of faces as needed by a specific boolean operation. Therefore, merging or stitching of BREPs into one BREP is delayed until after categorization and deletion of faces as specified by each boolean operation. The categorization of faces to carry out a specific boolean operation will be discussed next in Section 6.

Fig. 13. The edge recovery between each vertex pair in FIM. (a) The set of vertex pairs reflecting the signature of the intersection graph edges in a mesh face appearing in FIM. (b) After the splits, two edges are missing depicted as dashed lines. (c) The missing edges are recovered by the random swaps inside the planar boundary of the master face.

Fig. 14. The intersection graph edges on both BREPs are conformal and one-to-one matching. (a) The imprint of BREP2 onto BREP1 in terms of a set of mesh edges specified in FIM for each face. (b) The imprint of BREP1 onto BREP2 in terms of a set of mesh edges specified in FIM for each face.

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

19

6. Realizing boolean operations Due to the faceted nature of the input object boundaries, the boolean operations can be realized using set-operation rules on facet-sets of each BREP. The compilation of boundary portions of the two BREPs as inside, outside and common against each other is required to apply these rules. Therefore, the boundary qualification as to which portions are contained or common against the other BREP require categorization of BREP faces. There are a total of four different face types to reflect all possible combinations for two containment states, i.e., inside and outside for two BREPs;

   

BREP1 BREP1 BREP2 BREP2

faces faces faces faces

inside of BREP2 outside of BREP2 outside of BREP1 inside of BREP1

6.1. Categorization of faces Faces are categorized by tagging them with four integer values corresponding to the four distinct types itemized above. The categorization is a simple surface-walking algorithm over edgeconnected faces starting from one of the intersection-graph edges on any BREP. The initiating intersection graph edge shown in Fig. 15 is the edge connecting vertices V2 and V3. Note that there is also a pair of corresponding vertices on BREP2 for V2 and V3 and a corresponding edge since BREPs are not yet merged. Face F4 on BREP2 has a consistent orientation with respect to the region it

Fig. 16. Table for boolean operations and corresponding face tags. The deletion of faces tagged as shown in the second column as combinations of 10, 11, 20, 21 from the merged BREP results in the configurations corresponding to various boolean operations between BREP1 and BREP2 such as mesh implant, union, difference, intersect, imprint and non-regularized union operation as depicted in the first column.

is enclosing. Therefore, if the convention is such that the normal of the face is oriented outwards of the region it is enclosing then the box product of its oriented vertices and the opposite vertex of one of the adjacent faces on BREP1, say F2, will indicate if F2 is located outside or inside of BREP2 as depicted in Fig. 15. F2 is then tagged accordingly, for example, a positive box product value depicts that it is a BREP1 face outside of BREP2. The search continues with the tagged face of F2 and walks to the adjacent faces through its edges. The adjacent faces are tagged with the same tag until another intersection edge is crossed. The search then switches the tag to the opposite of that of the walked face. For instance, F3 is tagged as inside of BREP2 when crossed from F2 whose tag is outside of BREP2. The search process continues on both BREPs until all faces are tagged. 6.2. Employing boolean operation: deletion of faces There is a direct relationship between the face tags and a specific boolean operation. For instance, the boolean operation to union both BREPs would require for removing BREP1 faces inside BREP2 and BREP2 faces inside BREP1. This is equivalent of saying faces with corresponding tags, i.e., 11 and 21 should be deleted for the boolean union operation. The table showing which set of face tags should be deleted for a particular boolean operation is illustrated in Fig. 16. Note that for a non-regularized union, no faces are deleted. 6.3. Merging of the input BREPs

Fig. 15. Faces of BREP1 and BREP2 are categorized in four classes to carry out the boolean operations. (a) A cut-out mesh section which has the intersection graph edge passing through V2 and V3, three faces F1, F2 and F3 from BREP1 and three faces F4, F5, and F6 from BREP2. (b) The categorization is carried out by tagging the faces. To tag F2 as a face of BREP1 outside of BREP2, one needs to check the sign of the box product for the volume of the fictitious tetrahedron between F4 of BREP2 and the vertex V4 of BREP1 to be positive with the assumption that the BREP faces are oriented such that the face normals are showing outside of the region(s) of the BREP.

The final step of the discrete boolean algorithm is to stitch the two BREPs along their matching intersection graph edges. Both BREPs are already evaluated and modified by deleting the required faces for a specific boolean operation as explained in Section 6.2. The mesh data structure is modified to change the ownership of all mesh entities from BREP1 and BREP2 to a new BREP id. Also the inter-entity adjacencies are modified according to the topological vertex and edge merging on the mesh data structure. In the low level topological vertex and edge merging operation implemented in the mesh data structure, the mesh edges of BREP1 on the intersection graph is replaced by the corresponding BREP2 edges using the edge map. The final results are illustrated for the two inter-penetrating cylinder BREPs in Fig. 17 for every boolean operation. 7. Remeshing The mesh quality of the BREP resulting from the boolean is usually degraded due to the edge and face splits of the boolean

20

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

Fig. 17. Boolean results on two intersecting cylinder BREPs: (a) The two inter-penetrating BREPs. The results of boolean operations: (b) Imprint (non-regular union), (c) Union, (d) Subtract BREP2 from BREP1, (e) Subtract BREP1 from BREP2, (f) Intersect.

algorithm. Splits cause high vertex valence (vertex to edge connection number) and needle like sliver mesh face formations around the intersection graph as shown in Fig. 14. There have been several strategies proposed in improving the mesh quality after the booleans [13,15]. Explicit creation of intersection graph ring topologies (incremental approach) is advocated by these studies. Redistribution of the nodes around the intersection graph rings to create a uniform sized mesh was common in both of these approaches. Meshing strategy was to remove the faces locally around the intersection graph and fill the gap by an advancing front algorithm. Although this is a very fast way to generate quality meshes, the size of the gap is a concern since the geometry needs to be redefined in the gap area. If there is a huge discrepancy in the geometry of boundary edges defining the gap area, respecting the input geometry becomes highly susceptible. To eliminate such concerns, we have devised a remeshing approach using local mesh modifications which ensures recovery of the input discrete geometry as there is no gap formation. The mesh is modified incrementally by edge splits/swaps and collapses and vertex relocations that are defined over existing mesh topology by respecting the geometry of the mesh. Our remeshing strategy is similar to the works of [19] and [20] regarding the use of local edge splits/swaps and collapses. However, the former employs these operations over the parametrized 2D domain, and the latter uses an advancing front algorithm in 3D.

7.1. Rebuilding discrete geometry After remeshing, the discrete geometry topology is no longer a representative of geometry of the mesh; portions of the mesh may have been deleted (subtraction), carved out (intersections), or merged (union). Therefore, mesh needs to be reclassified by running the discrete geometry building DGB process (see Section 1) one more time to account for the change in the mesh topology impacted by the boolean operation. The mesh is disassociated from its discrete model first and then using the same angle thresholds and geometrical predicates, the DGB process is re-run again over the booleaned mesh topology. The new geometrical

Fig. 18. Copying the attributes over to the intersection booleaned BREP: (a) The frustum and sphere BREPs have user attributes A1, A2 and A3 defined on three geometry faces. (b) Intersection of these two BREPs preserves the attributes by reattaching them onto the newly created geometrical entities as a result of the DGB process performed after the boolean operation.

entities are defined within the context of the BREP of the boolean as shown in Fig. 18(b). 7.2. Attribute management The input meshes may have already had the user attributes associated with them. These include analysis-relevant data such as material properties, and boundary conditions. All such attributes are associated with the geometric model and not directly on the mesh by design. This is a key towards having a geometrybased analysis representation [5,6]. Therefore it is critical to ensure that any attribution data is not lost due to the booleans. As a result, before running the DGB process, a temporary look-up storage, a map for each user attribute corresponding to a set of classified mesh entities, is generated. After the DGB process, the data is attached back onto the new geometrical entities associated with the mesh entities defined in the temporary map. This is similar to the attribution management approaches implemented by CAD-systems [21–24] to handle topology modifying operations

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

such as the booleans. In such kernel implementations, the choices of transferring an attribute onto the new topologies, are often enumerated as options to the attribution itself- whether to keep or delete attributes from older component topologies. In other words, attributes are defined with their transfer option information so that they would be handled by the geometry modification operations accordingly. In our implementation, we have opted to transfer all attributes from both BREPs to the new BREP topologies as they are provided. A sample case depicting the transfer of attributes for the boolean intersection operation is shown in Fig. 18. 7.3. Assigning mesh sizes The edge and face splits of the boolean algorithm alter the sizing field of the input mesh. To improve the mesh quality, any remeshing strategy would need to have smoothly varying sizing field. We have opted to have the sizing implemented independent from meshing due to the following reasons:

 User has complete control of the sizing process.  There is no heuristics implicit to the meshers in terms of mesh sizes.

 Unbiased comparison of mesh quality can be done using different meshers.

 Experimenting on different sizing proxy implementations can be done efficiently without having to modify any meshing code. There have been many studies on mesh sizing approaches in the literature [25–31]. In many of these studies, background grids are used to encode a size distribution in space. They support a size distribution which is interpolated to provide the size value often by means of call-backs from the meshers. In the pioneering works of [32] and [26], boundary spacings are used to interpolate the mesh sizes. These are computed for each boundary vertex by averaging the lengths of the adjacent boundary edges. Edge meshing (sampling) followed by the triangulation of the domain using only the boundary vertices, enables computation of sizes by barycentric interpolation within each triangle (tetrahedron in 3D).

21

However, smoothness of the sizing field is not guaranteed and the sampling of the interior surface shape variations is also harder to accommodate with this formulation. Sampling of the geometry, driven by the user specified sizing parameters and storing these sampled sizes in an octree is another background grid alternative [33,34,27]. The smoothness of the sizing field is satisfied by forcing the octree refinement with 2:1 ratio everywhere in the tree. A problem with this approach is the inability to resolve sizes different than factors of two. There were also remedies to tackle this problem by scaling the entire octree. However, this extra refinement of the octree to satisfy 2:1 ratio and accurate size interpolation becomes another concern in terms of its large memory footprint that needs to be stored during the entire meshing cycle. Querying the size at a spatial location by the mesher has to be a very fast procedure as this is needed a lot during any meshing process. Therefore, the speed of querying sizes becomes another concern and a significant factor in the overall meshing time. Some researchers use auxiliary data structures such as uniform bins to aid the binary search time spent in the octree [20], and some use modified octree structures akin to meshes with octant adjacencies to retrieve data more quickly than logarithmic rate of the binary search [33]. None of these approaches are perfect but there is no solution without some limitation either. The advantages and drawbacks of various mesh sizing approaches have been explained in detail by the survey papers [28,29]. We have adopted a sizing implementation using sources [27,30,31]. In our implementation, user mesh sizing parameters, namely, global size, local sizes attached on geometrical entities, curvature refinement over curved surfaces, minimum number of edges on each geometrical edge with allowed minimum size, are all translated into generation of point, segment and triangular sources with constant and linearly varying strengths. The final mesh size is calculated from the contributions of all the sources. The minimum size prevails. The domain is divided into regular bins and each bin stores only a subset of the sources that will be effective to the size computation inside each bin’s bounding box to speed up the process as depicted in Fig. 19. This source based sizing implementation will be referred to as sizing-proxy in the next Section 7.4.

Fig. 19. Calculating mesh sizes using sources: (a) Linear and higher order spherical, segmented and triangular sources. (b) The closest point C is computed on each source Srci to the point P where the mesh size is requested. The mesh size from Srci is then calculated using the formula Sizei ¼ S þ Dnðf 1Þ with the geometrical progression factor of f. (c) The final mesh size at point P is calculated from the contributions of all the sources. The minimum size prevails. The domain is divided into regular bins and each bin stores only a subset of the sources that will be effective to the size computation inside its bounding box to speed up the process of computing sizes.

22

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

7.4. Local mesh modifications Remeshing algorithm is implemented using local mesh modifications in 3D space. Remeshing may be limited to the immediate zone of the intersection graph of the booleans, thereby not affecting the mesh outside of the impact zone of the edge/face splits. This approach is efficient due to its local nature, however, figuring out how far the remeshing will be performed moving out from the intersection graph rings is not always easy. Therefore, remeshing algorithm can have its input either in terms of a set of mesh face layers around the intersection graph or the entire mesh of the BREP. Regardless of this choice, details of local mesh modification operations used in surface remeshing algorithm will be the same. A source based mesh-sizing proxy provides a smooth sizing field for the remesher as explained in Section 7.3. Remesher then queries the sizing proxy for the target size at the mid location of each mesh edge. An edge split or collapse operation is performed depending on the comparison of the edge’s length versus the target size from the proxy as shown in Fig. 20. The goal is to adapt

the actual local mesh size to the target size of the proxy’s smooth sizing field. The edge splits are followed by the swaps of the opposite edges to the split vertex. The swap of the opposite edge is only realized if it minimizes the maximum face angles compared to the angles before the swap as shown in Fig. 21. Minimizing the maximum angles or maximizing the minimum

P Pnew

Ai

Ai

Ci

Ci

Pnew= (Ai Ci) / (Ai) Fig. 22. Local mesh quality optimization by vertex smoothing: (a) The location of point P is moved to Pnew due to smoothing of the vertex at P using either constrained Laplacian or CVT schemes. (b) Laplacian scheme: The new location Pnew is computed by area averaging the centroids of its surrounding mesh faces. (c) CVT scheme: The new location Pnew is computed by nearest neighbor area averaging of the center of circumcenters of its surrounding faces.

Fig. 20. Remeshing using local mesh modifications: (a) Mesh size at point P is queried from the sizing proxy at the mid of each mesh edge whose length is L. (b) If the length of the edge L is less than the required size then edge is collapsed. (c) If the edge length L is more than twice the required size, then edge is split at its mid location. (d) Opposite edges to the newly split vertex is checked if an edge swap would improve mesh quality. (e) The split vertex is smoothed locally using either the CVT or the constrained Laplacian schemes.

Fig. 21. Local mesh quality optimization by edge swaps: (a) The opposite edges to each vertex is checked if a swap of the edge would improve surrounding mesh face quality by maximizing the minimum internal angles aj to bj . (c) After the swap of the edge, new opposite edges to the vertex is pushed back into the list of edges to check for subsequent swap operations.

Fig. 23. The results of remeshing of boolean unioned of two penetrating cylinder BREPs. (a) Mesh of the boolean-unioned BREP. (b) Local remeshing applied around the intersection graph of unioning operation to improve mesh quality. (c) Remeshing applied globally to the BREP using the intersection graph local mesh size as segment sources. (d) Remeshing using uniform mesh size over the unioned BREP.

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

angles is the weak Delaunay-property argument for surface meshes [35,36]. Satisfying this weak Delaunay property is used to accomplish the main goal of improving the mesh quality. Finally, a vertex smoothing scheme is used to pre-condition the mesh topology for further swaps and splits. The idea in vertex smoothing is to relocate the vertex such that the worst quality face surrounding the vertex is improved. To this end, the new location is computed either by constrained Laplacian or centroidal voronoi tessellation (CVT) [37] schemes as shown in Fig. 22. The mesh of the BREP resulting from a boolean can be improved in different ways:

 Remesh in the vicinity of the intersection graph across a number of mesh face layers.

 Remesh the entire BREP using the mesh spacing of the intersection graph dictated as segment sources.

 Mesh the entire BREP with a uniform size or with any user sizing input using either remesher or any other available surface mesher.

The intersection graph mesh spacing is expected to be finer than the original size of the mesh. This local fine spacing can be respected by translating the fine sizes as segment sources with constant strength computed as the equal division of the total geometry edge length by the number of mesh edges due to splits. This natural way of generating segment sources in sizing proxy enables the remeshing algorithm to improve quality while maintaining smooth transition from the fine resolution of the intersection graph zone. The results from various remeshing options are depicted in Fig. 23.

23

8. Volume meshes The algorithm for discrete booleans can be extended to cover volume meshes. One or both BREPs can have volume meshes. Hence, a set of faces appearing in FIM may be adjacent to mesh regions (tetrahedral or hexahedral). These adjacent mesh regions need to be conformally modified when the face and edge splits are performed over the faces of FIM on geometry surfaces. Otherwise, the mesh becomes invalid and disconnected from its volume elements. Therefore, the main task is to match intersection graph edges while maintaining valid volume meshes. A similar problem has already been studied [18] for a different purpose of matching cavity boundaries adjacent to incomplete volume meshes. In the present case, we are trying to match the boundary of the mesh region to one or more of its faces appearing in FIM subject to edge and face splits. However, it is not necessary to modify the boolean algorithm until after the creation of matching intersection graph edges explained in Section 4. The boundary matching algorithm is defined as Simple Convex Triangulation (SCT) and depicted in Fig. 24. In SCT, there is an ordered manner of splitting of mesh region topologies by allowing the creation of non-conformal mesh entities temporarily. An example mesh face listed in FIM is split at the intersection locations of BREP1 and BREP2 as shown in Fig. 24(c). Convex triangulation is defined as an atomic operation connecting corner vertices of a topology (face/region) to a vertex created at its centroid. First, we need to make faces of the hexahedron conformal by applying convex triangulation for its faces. For example, one edge of the FIM face, is also adjacent to another interior side face of the hexahedron. The side interior face is ‘‘convex-triangulated’’ by connecting its corner vertices and the intersection vertex on the edge to the face centroid as depicted in Fig. 24(e). Next step is to apply the ‘‘convex-triangulation’’ to the region of the hexahedron itself by connecting all the corner and convex triangulation vertices of its faces to its centroid as depicted in Fig. 24(f). Finally, we have accomplished creation of conformal mesh region(s) with one or more of its face boundaries on the intersection graph. The rest of the boolean algorithm is identical. The only caveat is that when deleting the faces we also need to categorize the volume mesh regions and delete them in the same manner as depicted in Fig. 16.

9. Complex geometric operations and application to engineering design and analysis Boolean operations are indispensable asset in a geometric modeling framework. It is not possible to create complex shapes without the ability to apply boolean operations over primitive shapes. The basic boolean algorithm is designed to operate on two

Fig. 24. Simple-Convex-Triangulation (SCT) steps: (a) A face adjacent to a mesh region in intersection set. (b) Intersection edges. (c) Split of the face that keeps intersection edges. (d) Hexahedral mesh region adjacent to the face. (e) Convex triangulation of the faces of the mesh region constrained by edge splits connecting to face centers. (f) Faces being triangulated are connected to regions’s center to ensure a valid mesh.

Fig. 25. The result of discrete boolean unite applied to equally spaced tori around a circle. (a) Eight (8) disjoint tori, (b) Result of non-regular unioning followed by remeshing.

24

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

Fig. 26. Non-regular unioning of three BREPs. (a) Three penetrating BREPs: a sphere, a torus and a corner sheet. (b) Internal details after non-regular boolean operation. (c) Remeshing of the booleaned BREP using boundary edges as segment sources. (d) Using a frontal mesher with the same sizing.

Fig. 27. Boolean algorithm applied on BREPs with different mesh scales. (a) Three cylinder BREPs intersecting a box BREP with almost comparable sizing. (b) The result of boolean unite of all BREPs. (c) Bottom face of the box mesh after boolean. (d) Bottom face of the box mesh remeshed using the boundary segments as sources. (e) Three cylinder BREPs with same mesh size and a box BREP with a much coarser mesh size. (f) The result of unioning of all BREPs. (g) The bottom face of the cube after boolean operation. (d) Meshing the bottom face using a frontal approach with a uniform size and segment sources created around holes.

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

BREPs in this study. Nonetheless, it is straightforward to deal with models with multiple BREPs. Complex geometries can be created by the successive applications of booleans over more than two BREPs as depicted in Fig. 25. The resultant BREP at the end of each boolean operation, can be exploded into disjoint parts in terms of new BREPs to make the next boolean operation possible. The final stitched BREP from eight disjoint BREPs, in this example, is created by applying the boolean operation seven times between two BREPs consecutively, one BREP from the previous boolean and the BREP of an original torus as shown in Fig. 25(b). Booleans can also be applied on sheet or solid BREPs. The mixture of sheet and solid BREPs is also possible as depicted in Fig. 26 in which the non-regularized union between a solid sphere, a solid torus and a rectangular sheet-elbow BREP is realized. The mesh quality due to the splits at the intersection of these three (3) BREPs is vastly improved by remeshing depicted in Fig. 26(c). Boundary edge spacing along the intersection graph is distributed evenly in terms of segment sources to provide smooth sizing for the remesher. The same sizing proxy is used for a frontal meshing strategy for quality comparisons without having to change any meshing code as shown in Fig. 26(d). The effect of different input mesh scales on the boolean algorithm is depicted in Fig. 27 in which three slender cylinder BREPs are boolean unioned with a cube BREP of two significantly different mesh sizes. The algorithm is robust in dealing with different mesh scales in the input. Similar to previous tests, the intersection graph of the boolean result is used to create evenly distributed boundary spacing in terms of segment sources for remeshing. However, the choice is completely optional whether to create boundary spacing in the sizing proxy or remesh without it. The differences in mesh quality between the two different meshers can be compared as shown in Fig. 27(d) and (h). A realistic example of an application of interest to the military is depicted in Fig. 28 in which a component BREP representing a gun-system is implanted over the deck of a ship BREP. The guncomponent is moved to an appropriate location over the deck as

25

shown in Fig. 28(b). The two BREPs are then implanted and remeshed as depicted in Fig. 28(b) and in a close up in Fig. 28(e). Optionally, the implant can also be performed without impacting the mesh outside of the boolean intersection zone. The remeshing option that only modifies and improves the mesh around a few layers in the immediate vicinity of the intersection zone of the

Fig. 29. Input data courtesy of Naval Air Systems Command, Pax-River: Three air vehicle part BREPs are stitched to each other using boolean unite. (a) A conceptual engine, wing and aircraft body BREPs gathered into the same geometry model. (b) The boolean algorithm unite applied over three BREPs consecutively. (d) The detail of the mesh on the aircraft body shows the result of boolean algorithm.

Fig. 28. Input data courtesy of Naval Surface Warfare Center, Carderock Division: The implant of gun component mesh into the discrete ship mesh: (a) The individual component and ship BREPs. (b) The gun BREP is moved to the location that will be installed onto the ship deck. (c) The two BREPs are implanted. (d) The close-up mesh after implanting. (e) The close-up mesh after remeshing of the stitched BREPs.

26

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

Fig. 30. An automobile BREP is stitched with the floor of the wind tunnel mesh. (a) Automobile is moved to the middle of the floor. (b) The boolean algorithm unite applied to stitch the automobile tires over the floor of the wind tunnel. (c) The imprint of the flooring over the automobile tires. (d) Remeshing applied to improve mesh quality using the boundary spacing as segment sources for smooth mesh sizing around the imprint of the tires with the floor. (e) Zoom-out of the mesh around the imprint zone.

two BREPs is obviously the least intrusive in terms of preserving the input mesh topology as shown in Fig. 28(d). Another real-world example in aviation research is depicted in Fig. 29 demonstrating two boolean union operations for three (3) BREPs between a gas turbine engine with its pylon to a wing followed by another union with the aircraft body. The detail of the mesh and the result of the union between the aircraft body and the wing section is shown in Fig. 29(c). Similar to the aviation example, an automobile BREP is stitched to the floor of the wind tunnel. The remeshing greatly improves the mesh quality as depicted in Fig. 30(d, e).

10. Closing remarks The ability to apply various boolean operations to create complex shapes from simpler forms is a valuable asset and a key component for any modeling suite for engineering analysis and design. In this paper, we have devised and demonstrated a robust algorithm to perform discrete booleans on solid or sheet BREPs with differing input mesh scales. The algorithm is also capable of preserving user attributes by associating them with the newly created booleaned topologies efficiently. It is implemented in Cþþ as part of a geometry, mesh and attribution modeling platform called ‘‘Capstone’’. We point our several future enhancements to the process of discrete boolean algorithm, including improved remeshing using local re-parametrization, more flexible attribute management instead of migrating all attributes, and rebuilding the discrete geometry locally. Remeshing with local parametric projections may help improving the mesh quality for two reasons: more accurate projections to the original BREP topologies, and, a more precise computation for CVT which is not a well-defined scheme over 3D surface meshes. More options in attribute migration from operand BREPs over to the resultant BREP is also required similar to CAD operations which will increase the user control over the attribute transfer mechanism. Globally rebuilding the geometry

topology is a safe and robust solution to create new BREP topologies, however, it is not the most efficient resolution. A more efficient approach could be to run DGB locally only within the closure of the original topologies impacted by the boolean operation. Finally, an important and necessary improvement is a volume remeshing algorithm using local modification to increase the mesh quality after SCT for volume and mixed dimension meshes.

Acknowledgments This work was funded by the DoD High Performance Computing and Modernization Program (HPCMP) as part of the Computational Research & Engineering Acquisition Tools and Environment (CREATE) Meshing and Geometry (MG) Project. References [1] S. Owen, D.R. White, Mesh-based geometry: a systematic approach to constructing geometry from a finite element mesh, in: Tenth International Meshing Roundtable, 2001, pp. 83–96. [2] M. Mantlya, An Introduction to Solid Modeling, Computer Science Press, 1988. [3] W.J. Schroeder, M.S. Shephard, On rigorous conditions for automatically generated finite element meshes, in: J. Turner, J. Pegna, M. Wozny (Eds.), Product Modeling for Computer-Aided Design and Manufacturing, 1991, pp. 267–281. [4] S. Dey, M.S. Shephard, M.K. Georges, Elimination of the adverse effects of small model features by the local modification of automatically generated meshes, Eng. Comput. 13 (1997) 134–152. [5] B.K. Karamete, Design and implementation of a new topological data structure: double link topology structure (DLTS), in: Eleventh U.S. National Congress on Computational Mechanics, Minnesota, MN, 2011. [6] M.W. Beall, M.S. Shephard, A general topology-based mesh data structure, Int. J. Numer. Methods Eng. 40 (19) (1997) 1573–1596. [7] J. Corney, 3d Modeling Using the Acis Kernel and Toolkit, 1st edn., John Wiley, 1997. [8] S.H. Lo, W.X. Wang, A fast robust algorithm for the intersection of triangulated surfaces, Eng. Comput. 20 (2004) 11–21. [9] S.H. Lo, W.X. Wang, Finite element mesh generation over intersecting curved surfaces by tracing of neighbours, Eng. Comput. 41 (2005) 351–370.

B. Kaan Karamete et al. / Finite Elements in Analysis and Design 68 (2013) 10–27

[10] M. Aftosmis, M.J. Berger, J.E. Melton, Robust and efficient cartesian mesh generation for component-based geometry, AIAA J. 36-6 (1998) 952–960. [11] O. Devilers, P. Guigue, Faster Triangle–Triangle Intersection Tests, Technical Report. 4488, INRIA, 2002. [12] H. Samet, Design and Analysis of Spatial Data Structures, Addison Wesley, 1990. [13] A.A. Shostko, R. Lohner, W.C. Sandberg, Surface triangulation over intersecting geometries, Int. J. Numer. Methods Eng. 44 (1999) 1359–1376. [14] K. Guo, L.C. Zhang, C. Wang, S.H. Huang, Boolean operations of STL models based on loop detection, Int. J. Adv. Manuf. Technol. 33 (2007) 627–633. [15] S.H. Lo, Automatic mesh generation over intersecting surfaces, Int. J. Numer. Methods Eng. 38 (1995) 943–954. [16] B.K. Karamete, R.V. Garimella, M.S. Shephard, Recovery of an arbitrary edge on an existing surface mesh using local mesh modifications, Int. J. Numer. Methods Eng. 50 (2001) 1389–1409. [17] P. George, F. Hecht, E. Saltel, Automatic mesh generator with specified boundaries, Comput. Methods Appl. Mech. Eng. 92 (1991) 269–288. [18] B.K. Karamete, M.W. Beal, M.S. Shephard, Triangulation of arbitrary polyhedra to support automatic mesh generators, Int. J. Numer. Methods Eng. 49 (2000) 167–191. [19] J. Remacle, C. Geuzaine, G. Compe re, E. Marchandise, High quality surface remeshing using harmonic maps, Int. J. Numer. Methods Eng. 83 (4) (2010) 403–425. [20] R. Aubry, G. Houzeaux, M. Va´zquez, A surface remeshing approach, Int. J. Numer. Methods Eng. 85 (12) (2011) 1475–1498. [21] C. Spatial, 3D InterOp Data Exchange, URL /http://www.spatial.com/3d-inter operabilityS, 2011. [22] S.P. Siemens, Parasolid interoperability, URL /http://www.plm.automation. siemens.comS, 2011. [23] C. PTC, Interoperability for Pro/ENGINEER, URL /http://www.ptc.comS, 2011. [24] SMLIB, Data Translators, URL /http://www.smlib.comS, 2011.

27

[25] N. Weatherill, M. Marchant, O. Hassan, D. Marcum, Grid adaptation using a distribution of sources applied to inviscid compressible flow simulations, Int. J. Numer. Methods Fluids 19 (1994) 739–764. [26] J. Peraire, J. Peiro, K. Morgan, Adaptive remeshing for three-dimensional compressible flow computations, J. Comput. Phys. 103 (1992) 269–285. ¨ [27] R. Lohner, Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods, 2nd edn., Wiley, 2008. [28] S.J. Owen, S. Saigal, Neighborhood-based element sizing control for finite element surface meshing, in: Sixth International Meshing Roundtable, 1997, pp. 143–154. [29] S. Owen, S. Saigal, Surface mesh sizing control, Int. J. Numer. Methods Eng. 47 (2000) 497–511. [30] A. Cunha, S. Canann, S. Saigal, Automatic boundary sizing for 2D and 3D meshes, in: AMD Trends in Unstructured Mesh Generation, ASME 220, 1997, pp. 65–72. [31] W. Quadros, V. Vyas, M. Brewer, S. Owen, K. Shimada, A computational framework for generating sizing function in assembly meshing, in: Fourteenth International Meshing Roundtable, 2005, pp. 55–73. [32] N. Weatherill, O. Hassan, Efficient three dimensional Delaunay triangulation with automatic point creation and imposed boundary constraints, Int. J. Numer. Methods Eng. 37 (1994) 2005–2039. [33] M. Shephard, M. Georges, Automatic three-dimensional mesh generation by the finite octree technique, Int. J. Numer. Methods Eng. 32 (1991) 709–749. ¨ [34] R. Aubry, K. Karamete, E. Mestreau, F. Bulat-Jara, S. Dey, R. Lohner, Geodesicsbased surface reparametrization, in: 50th AIAA, pp. 1258–2012. [35] S. Sloan, A fast algorithm for generating constrained Delaunay triangulations, Comput. Struct. 47 (3) (1993) 441–450. [36] H. de Cougny, Parallel Unstructured Distributed Three Dimensional Mesh Generation, Ph.D. Thesis, Mechanical Engineering, Aeronautical Engineering and Mechanics, Rensselear Polytechnic Institute, Scientific Computation Research Center, Troy NY, 1998. [37] Q. Du, M. Gunzburger, Grid generation and optimization based on centroidal voronoi tessellations, Appl. Math. Comput. 133 (2002) 591–607.