The meccano method for automatic 3-d triangulation and volume parametrization of complex solids

Abstract In this paper we review the novel meccano method. We summarize the main stages (subdivision, mapping, optimization) of this automatic tetrahedral mesh generation technique and we concentrate the study to complexgenus-zero solids. In this case, our procedure only requires a surface triangulation of the solid. A crucial consequence of our method is the volume parametrization of the solid to a cube. We construct volume T-meshes for isogeometric analysis by using this result. The efficiency of the proposed technique is shown with several examples. A comparison between the meccano method and standard mesh generation techniques is introduced. Keywords: Tetrahedral mesh generation, adaptive refinement/derefinement, nested meshes, mesh smoothing, mesh untangling, surface and volume parametrization. 1 Introduction Manyauthorshave devotedgreateffort to solvingtheautomaticmeshgenerationproblem in different ways [4, 9, 21, 22, 37, 38, 39, 41], but the 3-D problem is still open [2, 3]. In the past, the main objective has been to achieve high quality adaptive meshes of complex solids with minimal user intervention and low computational cost. At present, it is well known that most mesh generators are based on Delaunay triangulation and advancing front technique. However, problems related to mesh quality, mesh adaptation and mesh conformity with the solid boundary, still remain. We have recently introduced the meccano technique in [6, 7, 33, 34] for constructing adaptive tetrahedral meshes of solids. The method requires a surface triangulation of the solid, a meccano and a tolerance that fixes the desired approximation of the solid surface. The name of the method stems from the fact that the process starts from an outline of the solid, i.e. a meccano composed by connected polyhedral pieces. The method builds a 3-D triangulation of the solid as a deformation of an appropriate tetrahedral mesh of the meccano. The meccano can be made up of different types of pieces (cuboids, pyramids, prisms, etc). A particular case is a meccano consisting only of connected cubes, i.e. a polycube [28, 40, 44]. The main idea of the new mesh generator is to combine an automatic parametrization of surface triangulations [14], a local refinement algorithm for 3-D nested triangulations [26] and a simultaneous untangling and smoothing procedure [10]. Those previous procedures (mapping, refinement, untangling and smoothing) are not in themselves new, but the overall integration is an original contribution. In this paper, we present significant advances in the method. We define an automatic parametrization of a solid surface triangulation to the meccano boundary. For this purpose, we first divide the surface triangulation into patches with the same topological connection as the meccano faces. Then, a discrete mapping from each surface patch to the corresponding meccano face is constructed by using the parameterization of surface triangulations proposed in [14, 15, 16, 17]. The shape-preserving parametrizations, which are planar triangulations on the meccano faces, are the solutions of linear systems based on convex combinations. Specifically, we describe the procedure for a solid whose boundary is a surface of genus 0; i.e. a surface that is homeomorphic to the surface of a sphere. In this case, the meccano is a single cube, and the global mapping is the combination of six patch-mapping. The solution to several compatibility problems on the cube edges will be discussed. The extension to more general solids is possible if the construction of an appropriate meccano is assumed. In the near future, more effort should be made in an automatic construction of the meccano when the genus of the solid surface is greater than zero. Currently, several authors are working on this aspect in the context of polycube-maps, see for example [28, 40, 44]. They are analyzing how to construct a polycube for a generic solid and, simultaneously, how to define a conformal mapping between the polycube boundary and the solid surface. Although surface parametrization has been extensively studied in the literature, only a few works deal with volume parametrization and this problem is still open. A meshless procedure is presented in [27] as one of the first tentative approaches to solve the problem. In addition, [18] gives a simple counterexample to show that convex combination mappings over tetrahedral meshes are not necessarily one-to-one. Our method automatically obtains a volume parametrization of the solid. In this paper we introduce this result for isogeometric analysis [3, 8]. In the following section we present a brief description of the main stages of the method for a generic meccano composed of polyhedral pieces. In Section 3 we analyze the algorithm in the case that the meccano is formed by a simple cube. In Section 4 we show test problems and practical applications that illustrate the efficiency of this strategy. Finally, the conclusions and future research are presented in Section 5. 2 The Algorithm of the Meccano Method Themain stepsofthemeccano tetrahedralmeshgenerationalgorithmaresummarized in this section. A first approach of this method can be found in [33, 6, 7, 34]. The input data of the algorithm are the definition of the solid boundary (for example a surface triangulation) and a given precision (corresponding to the approximation of the solid boundary). The following algorithm describes the mesh generation approach. Meccano tetrahedral mesh generation algorithm 1. Construct a meccano approximation of the 3-D solid formed by polyhedral pieces. 2. Defineanadmissible mappingbetweenthemeccanoboundaryfaces and the solid boundary. 3. Construct a coarse tetrahedral mesh of the meccano. 4. Generatea local refined tetrahedral mesh of the meccano, such that the mapping (according step 2) of the meccano boundary triangulation approximates the solid boundary for a given precision. 5. Move the boundary nodes of the meccano to the solid surface according to the mapping defined in 2. 6. Relocate the inner nodes of the meccano. 7. Optimize the tetrahedral mesh by applying the simultaneous untangling and smoothing procedure. The first step of the procedure is to construct a meccano approximation by connecting different polyhedral pieces. The meccano and the solid must be equivalent from a topological point of view, i.e., their surfaces must have the same genus. Once the meccano is assembled, we have to define an admissible one-to-one mapping between the boundary faces of the meccano and the boundary of the solid. In step 3, the meccano is decomposed into a coarse tetrahedral mesh by an appropriate subdivision of its initial polyhedral pieces. This mesh is locally refined and its boundary nodes are virtually mapped to the solid surface until it is approximated to within a given precision. Then, we construct a mesh of the domain by mapping the boundary nodes from the meccano faces to the true surface and by relocating the inner nodes at a reasonable position. After those two steps, the resulting mesh is generally tangled. Therefore, a simultaneous untangling and smoothing procedure is applied and a valid adaptive tetrahedral mesh of the solid is obtained. 3 The Meccano Method for Genus-Zero Solids In this section, we present the application of the meccano algorithm in the case of the solid surface being genus-zero and the meccano being formed by one cube. We assume a triangulation of the solid surface as a datum. We introduce an automatic parametrization between the surface triangulation of the solid and the cube boundary. To that end, we divide the surface triangulation into six patches, with the same topological connection that cube faces, so that each patch is mapped to a cube face. We note that the meccano method constructs a high quality tetrahedral mesh even when the initial triangulation has poor quality. 3.1 Meccano A simple cube, C, is defined as meccano. We associate a planar graph, GC, to the meccano in the following way: •Each face of the meccano corresponds to a vertex of the graph. •Two vertices of the graph are connected if their corresponding meccano faces share an edge. Figure 1 shows the numbering of cube faces and their connectivities, and Figure 2 represents the corresponding planar graph. The position of the cube is crucial to define an admissible mapping between the cube and solid boundary, as we analyze later. However, its size is less important, because it only affects the efficiency of the mesh optimization step. If the center of the cube is placed inside the solid, the existence of an admissible mapping is ensured. (a) (b) Figure 1: Meccano formed by one cube: (a) notation of nodes and faces of the cube and (b) connectivities of faces 3.2 Mapping from Cube Faces to Solid Surface Patches Once the cube is fixed, we have to determine a mapping between the cube faces and the solid surface triangulation. First, we define the concept of admissible mapping for a cube. Let ΣCbe the boundary of the cube and ΣSthe boundary of the solid, Figure 2: Planar graph GCassociated to the cube given by a surface triangulation TS. We denote by Σi Cthe i-th face of the cube, i.e. ΣC=5 i=0 Σi C. Let Π:Σ C→ΣSbe a piecewise function, such that Π|Σi C=Π iwhere Πi:Σ i C→Πi(Σi C)⊂ΣS. Then, Πis called an admissible mapping if it satisfies: a) Functions {Πi}5 i=0 are compatible on ΣC. That is Πi |Σi C∩Σj C =Π j |Σj C∩Σi C ,∀i, j = 0,...,5, with i=jand Σi C∩Σj C=∅. b) Global mapping Πis continuous and bijective between ΣCand ΣS. Note that admissible mapping is not unique. We define an automatic admissible mapping in the following sections. For this purpose, we first construct a partition of the solid surface triangulation into six patches, maintaining the topology of the graph of the Figure 2, and then we parametrize each patch to a cube face. 3.2.1 Partition of the Solid Surface Triangulation In the following, we call connected subtriangulation a set of triangles of TSwhose interior is a connected set. Given a decomposition of the surface triangulation TSin any set of connected subtriangulations, we can associate a planar graph, GS, to this partition in the following way: •Each subtriangulation corresponds to a vertex of the graph. •Two vertices of the graph are connected if their corresponding subtriangulations have at least one common edge. We say that a solid surface partition and the meccano are compatible if their graphs are isomorphic, GS=GC. In our case, since the solid surface is isomorphic to a sphere, it is clear that a compatible partition exists. We now propose an algorithm to obtain a decomposition of the given solid surface triangulation TSinto six subtriangulations {T i S}5 i=0. We distinguish three steps: a) Subdivision in connected subtriangulations. We construct the Voronoi diagram associated to the centers of the six cube faces. We consider that a triangle F∈T Sbelongs to the i-th Voronoi cell if its barycenter is inside this cell. We generate a partition of TSin maximal connected subtriangulations with this criterion, i.e. two subtriangulations belonging to the same cell can not be connected. We denote as Tij Sthe j-th connected subtriangulation belonging to the i-th Voronoi cell, and niis the total number of subtriangulation in the i-th cell. The number of subtriangulation depends only on the position of the cube center. If this point is placed inside the solid and the surface triangulation is fine enough, there is a subtriangulation in each cell. If any ni=0the process is aborted and the center of the meccano must be modified. For a complex solid the value of niis usuallygreater than 1, therefore, a modification of the partition is necessary to obtain a compatible decomposition of TS. b) Construction of the graph. We associate a planar graph, GSto the partition generated in the previous step. If the center of the cube is inside the solid and the surface triangulationis fine enough, there isone compatiblesubtriangulation for each Voronoi cell, i.e. there is one head subtriangulation Ti0 S, vertex of the graph GS, with the same connection as the vertex associated to the i-th cube face in GC. In other case, the subtriangulation with the greatest number of elements is chosen as Ti0 S. c) Reduction of the graph. In order to achieve a decomposition of TS, we propose an iterative procedure to reduce the current graph GS. In each step all triangles of Tjk Sare included in the head subtriangulation Ti0 Sif: –Ti0 Sis the head subtriangulation with the fewest triangles. –Tjk Sand Ti0 Sare connected. –kis higher than zero. Then, the vertex Tjk Sis removed from the graph and its connectivities are inherited by Ti0 S. The connectivity of the graph is updated. After this process, Ti0 Scould be connected to other subtriangulations Til Sof the same i-th cell. In this case, the triangles of all Til Sare included in Ti0 S, the graph vertices Til Sare removed from the graph, their connectivities are inherited and the graph connectivitiesare updated. Therefore, the connected subtriangulations are always maximal in all algorithm steps. This procedure continues iteratively until the graph comprises only six head vertices. Although the proposed algorithm obtains a six vertex graph GS, its compatibility with the cube graph GCcan notbe ensured. As the computational costof this algorithm islow, a movementinthe cube center inorder toobtain a compatiblepartition{T i S}5 i=0, does not affect the efficiency of the meccano technique. In fact, this procedure could be guided by a genetic algorithm [43]. In what follows Σi Sis the solid surface patch defined by the triangles of Ti S. 3.2.2 Parametrization of the Solid Surface Triangulation Once the given solid surface ΣSis decomposed into six patches Σ0 S,...,Σ5 S, we map each surface patch Σi Sto the corresponding cube face Σi Cby using the parametrization of surface triangulations proposed by M. Floater; in [14] the author gives a method to compute a planar triangulation on a convex domain that is isomorphic to a simply connected surface triangulation. In our case the convex domain is a cube face Σi C, and the surface triangulation is Ti S. Then, we define Πi−1:Σ i S→Σi C as the parametrization introduced by Floater, and we denote τi F=(Π i)−1(Ti S)as the planar triangulationof Σi Cassociated to Ti S. To obtain τi F, Floater parametrizationfixes their boundary nodes and the position of their inner nodes is given by the solution of a linear system based on convex combinations. Let {Pi 1,...,Pi n}be the inner nodes and {Pi n+1,...,Pi N}be the boundary nodes of Ti S, respectively, where Ndenotes the total number of nodes of Ti S. Once the boundary nodes {Qn+1,...,Q N}of τi Fare fixed, the position of the inner nodes {Qi 1,...Q i n}is given by the solution of the system: Qi k= N  l=1 λklQi l,k=1,...,n. The values of {λkl}l=1,...,N k=1,...,n are the weights of the convex combinations, such that λkl =0,if Pkand Plare not connected λkl >0,if Pkand Plare connected N  l=1 λkl =1,for k=1,...n. In [14] three alternatives are analyzed: uniform parametrization, weighted least squares of edge lengths and shape preserving parametrization. Another choice, called mean value coordinate, is presented in [16]. The goal is to obtain an approximation of a conformal mapping. In order to ensure the compatibility of {Πi}5 i=0, the boundary nodes of {τi F}5 i=0 must coincide on their common cube edges. The six transformations {Πi}5 i=0 define an admissible mapping between ΣCand ΣS, i.e. the cube boundary triangulation τF=5 i=0 τi Fis a global parametrization of the solid surface triangulation TS. Two important properties of mapping Πare: (a) the triangulations τFand TShave the same topology, (b) each triangle of τFis completely contained in one face of the cube. We note that usual polycube-maps [40, 28] verify property (a), but they do not verify property (b), i.e., a triangle belonging to TScan be transformed by a polycube-map into a triangle whose vertices are placed on different faces of the polycube. The proposed mapping Πis used in a following step of the meccano algorithm to map a new triangulation τK(obtained on ΣCby application of the refinement algorithm of Kossaczky [26]) to the solid boundary. Several problems can appear in the application of this transformation and their solutions will be commented on next. Basically, these problems are due to the fact that a valid triangulation τK=τFon ΣCcan be transformed by Πinto a non-valid one on the solid surface. 3.3 Coarse Tetrahedral Mesh of the Meccano We build a coarse and high quality tetrahedral mesh by splitting the cube into six tetrahedra [26]. For this purpose, it is necessary to define a main diagonal on the cube and corresponding diagonal on its faces, see Figure 3(a). The resulting mesh can be recursively and globally bisected [26] to fix a uniform element size in the whole mesh. Three consecutive global bisections for a cube are presented in Figures 3 (b), (c) and (d). The resulting mesh of Figure 3(d) contains 8cubes similar to the one shown in Figure 3(a). Therefore, the recursive refinement of the cube mesh produces similar tetrahedra to the initial ones. (a) (b) (c) (d) Figure 3: Refinement of a cube by using Kossaczky’s algorithm: (a) cube subdivision into six tetrahedra, (b) bisection of all tetrahedra by inserting a new node in the cube main diagonal, (c) new nodes in diagonals of cube faces and (d) global refinement with new nodes in cube edges 3.4 Local Refined Tetrahedral Mesh of the Meccano The next step in the meccano mesh generator includes a recursive adaptive local refinement strategy, by using Kossaczky’salgorithm[26], of those tetrahedra with a face S 6 abc a’ c’ b’ P P’ Figure 4: Mapping Πfrom the boundary nodes of the meccano a, b, c, P to the points a,b ,c ,Pof the solid surface placed on a boundary face of the initial coarse tetrahedral mesh of the cube. The refinement process is done in such a way that the given solid surface triangulation TS is approximated by a new triangulation within a given precision. That is, we seek an adaptive triangulation τKon the cube boundary ΣC, so that the resulting triangulation after node mapping Π(τK)is a good approximation of the solid boundary. The user has to introduce as input data a parameter ε, which is a tolerance to measure the separation allowed between the linear piecewise approximation Π(τK)and the solid surface defined by the triangulation TS. At present, we have considered two criteria: the first related to the Euclidean distance between both surfaces and the second attending to the difference in terms of volume. To illustrate these criteria, let abc be a triangle of τKplaced on the meccano boundary, and abcthe resulting triangle of Π(τK)after mapping the nodes a,band con the givensolid surface ΣS, see Figure 4. We define two different criteria to decide whether it is necessary to refine the triangle (and consequently the tetrahedron containing it) in order to improve the approximation. For any point Qin the triangle abc, we define d1(Q)as the euclidean distance between the mapping of Qon ΣS,Q, and the plane defined by abc. This definition is an estimateof the distance between the surface of the solid and the current piecewise approximation Π(τk). We also introduce a measure in terms of volume and then, for any Qin the triangle abc, we define d2(Q)as the volume of the virtual tetrahedron abcQ. In this case, d2(Q)is an estimate of the lost volume in the linear approximation by the face abc of the solid surface. The threshold of whether to refine the triangle or not is given by a tolerance εifixed by the user. With the previous definition, d1(Q)<ε 1for all Qin the boundary on the cube implies that the distance between the surface of the solid and its new piecewise alternative is to optimize the initial parametrization τF, allowing the node movement between cube faces, as is usual in polycube maps applications [40, 28, 44]. 4 Applications We have implemented the meccano technique using: •The parametrization toolbox of the geometrygroup at SINTEF ICT, Department of Applied Mathematics. •The module of 3D refinement of ALBERTA code. •Our optimization mesh procedure describes in section 3.7. The parametrization of a surface triangulation patch Ti Sto a cube face Σi Cis done with GoTools core and parametrization modules from SINTEF ICT, available on the website http://www.sintef.no/math software. This code implements Floater’s parametrization in C++. Specifically, in the following applications we have used the mean value method for the parametrization of the inner nodes of triangulation, and the boundary nodes are fixed with chord length parametrization [14, 16]. ALBERTA is an adaptive multilevel finite element toolbox [36] developed in C. This software can be used for solving several types of 1-D, 2-D or 3-D problems. ALBERTA uses the Kossaczky refinement algorithm [26] and requires an initial mesh topology [36]. The recursive refinement algorithm could not terminate for general meshes. The meccano technique constructs meshes that verify the imposed restrictions of ALBERTA relative to topology and structure. They can be refined by its recursive algorithm, because they are loop-free, and the degeneration of the resulting triangulationsafter successive refinements isavoided. The minimumquality ofrefined meshes is function of the initial mesh quality [30, 42]. The performance of our novel tetrahedral mesh generator is shown in the following applications. All of them are complex genus-zero solid that are classical test in mesh generation analysis. The first example corresponds to a screwdriver and the second to the armadillo. In addition, we present an application of the meccano method to construct volume Tmeshes for isogeometric analysis. Particularly, we consider the volume parametrization of the Stanford bunny. We have obtained a surface triangulation of these objects from internet. 4.1 Screwdriver The original surface triangulation of the screwdriver has been obtained from CNR IMATI Shapes Repository, http://shapes.aim-at-shape.net,and it is shown in Figure 8 (a). It has 27150 triangles and 13577 nodes. The bounding box of the solid is defined by the points (x, y, x)min =(−15,−26,−15) and (x, y, z)max =(15,44,15). (a) (b) Figure 8: (a) Compatible partitionofthe originalsurface triangulation,(b) Screwdriver meccano, the Floater’s parametrization is showed on the prism faces We consider a regular prism as a meccano (see Figure 8(b)). Its dimensions are 4×48 ×4and its center is placed inside the solid at the point (0,14,0). Although a unit cube could be considered as meccano, we have decided to use a narrow meccano because this choice slightly improves the quality of the final mesh. In any case, both meccanos (cube or prism) have associated the same planar graph. We obtain an initial subdivisionof screwdriver surface in six connected subtriangulations using Voronoi diagram of prism face centers. It is a compatible decomposition of the surface triangulation, according 3.2.1. We then remove the dividing edges and smooth the interfaces between patches by minimizing the zigzag shapes. Figure 8(a) shows the resulting compatible partition {T i S}5 i=0. We map each surface patch Σi Sto the prism faces Σi Cby using the Floater parametrization [14]. The definition of the one-to-one mapping between the cube and screwdriver boundaries is straightforward once the parametrization of the screwdriver surface triangulation is built, see Figure 8(b). A detail of the two more significant parametrizations are shown in Figure 9. (a) (b) Figure 9: Detail of Floater’s parametrization of two screwdriver patches to square meccano faces The prism is divided into 12 cubes, and then in 72 tetrahedra. We now fix several tolerances: ε2=1,0.1,0.01,0.001 and generate the corresponding tetrahedral meccano meshes. In Table 1 we report the main features of them. For example, for a tolerance ε2=0.001, our method applies 51 Kossaczky recursive bisections to generate a local refined mesh that contains 168834 tetrahedra and 39617 nodes, with 36968 triangles and 18486 nodes on its boundary. ε2=1 ε2=0.1 ε2= 0.01 ε2= 0.001 # nodes 1955 4430 12783 39617 # tetrahedra 8118 18814 54276 168834 # nodes on boundary 941 2047 5979 18486 # triangles on boundary 1878 4090 11954 36968 Kossaczky refinement 42 45 48 51 Table 1: Main features of the screwdriver tetrahedral meshes generated by meccano method for tolerances ε2=1,0.1,0.01,0.001 In order to obtain a tetrahedral mesh of the screwdriver we have to apply the procedures: external node mapping, relocation on inner nodes and mesh optimization. The Floater’s parametrizations allow us to map the meccano external nodes to the screwdriver surface, but the resulting tetrahedral mesh is complety tangled. For example, in Figure 10(a) we show the tangled mesh for the tolerance ε2=0.001. No relocation is applied for the coarser tolerance (ε2=1), and we use the volume parametrization (from the screwdriver to the meccano) given by this coarse approximation to relocate the inner nodes in the other cases. The relocation procedure significatively reduces the number of inverted tetrahedra but it does not solve the problem: in the case ε2=0.001 the number of inverted tetrahedra decreases from 20875 to 4961 (see Table 2). We now use the tetrahedral mesh optimization. Although no relocation is applied in the coarse mesh, the optimization procedure untangles it in only 13 iterations. The other meshes are untangled in no more than 6 iterations. All data are reported in Table 2. We could improve the behaviour of this procedure, if we relocate the inner nodes using the volume parametrization defined by the previous value of the tolerance, i.e the inner nodes of mesh ε2=0.01 relocated with the volume parametrization obtained for ε2=0.1. However, we have decided to use the coarser approximation in order to show the robustness of the optimization algorithm introduced in [10]. Finally, we apply 5iterations to smooth the meshes. In all cases, the resulting mesh quality is improved to a minimum value about 0.2and an average about qκ=0.7.We note that the meccano technique generates a high quality tetrahedral mesh, see Figure 10(b). In order to show the efficiency of the mesh optimization technique inside the screwdriver we display in Figure 11 two sections before (a) and after (b) its application for ε2=0.001. In addition, we show the sequence of screwdriver approximation for different values of ε2in Figure 12. We also note that, due to the high quality surface triangulation obtained with our method, the mesh improvement is not significant if we previouslyapply the smoothing surface triangulation algorithm introduced in [11]. The computations have been done on a Dell precision 690, 2 Dual Core Xeon pro- cessor and 8Gb RAM memory. The CPU times for constructing the final meshes of the screwdriver are also reported in Table 2. The most demanding example requires approximately 42 second. More precisely, the CPU time in this case of each step of the meccano algorithm is: 0.5seconds for the subdivision of the initial surface triangulation into six patches, 0.9seconds for the Floater parametrization, 18.6seconds for the Kossaczky recursive bisections, 2.3seconds for the external node mapping and inner node relocation, and 19.7seconds for the mesh optimization. (a) (b) Figure 10: Screwdriver tetrahedral meshes after (a) external node mapping and after (c) the application of the mesh optimization procedure for a tolerance ε2=0.001 4.2 Armadillo We now consider the Armadillo. The original surface triangulation has been obtained from http://graphics.stanford.edu/data/3Dscanrep/ , i.e. the Stanford Computer Graphics Laboratory. It has 30000 triangles and 15002 nodes. The bounding box of the solid is defined by the points (x, y, x)min =(−60,−50,−26) and (x, y, z)max =(68,66,90). ε2=1 ε2=0.1ε2=0.01 ε2=0.001 Inverted elements before/after relocation 1263 2407/11 6358/417 20875/4961 Untangling iterations 13 1 3 6 Smoothing iterations 5 5 5 5 Minimum quality 0.16 0.21 0.21 0.18 Mean quality 0.68 0.69 0.71 0.73 Optimization time (seconds) 2 1 5 19 Total time (seconds) 14 15 21 42 Table 2: Relocation and optimization data for the screwdriver meshes based on the volume parametrization obtained with ε2=1 (a) (b) Figure 11: Section of the screwdriver tetrahedral meshes before (a) and after (b) the application of the mesh optimization procedure for a tolerance ε2=0.001 (a) (b) (c) (d) Figure 12: Screwdriver approximations for tolerances (a) 1, (b) 0.1, (c) 0.01, (d) 0.001 (a) (b) Figure 13: (a) Original surface triangulation of the Armadillo, (b) resulting valid tetrahedral mesh generated by the meccano method We consider a unit cube as meccano. Its center is placed inside the solid at the point (7.5,17.5,55.5). We obtain an initial subdivision of Armadillo surface in eleven maximal connected subtriangulations using Voronoi diagram. We reduce the surface partition to six patches and construct the Floater parametrization from each surface patch to the corresponding cube face. (a) (b) (c) Figure 14: Cross sections of the Armadillo tetrahedral meshes before (a) and after (b) mesh optimization, and (c) same section of the meccano mesh Fixing a tolerance ε2=0.1, the meccano methodgenerates a tetrahedral mesh with 144964 tetrahedra and 33889 nodes; see Figure 13(b). This mesh has 31316 triangles and 15660 nodes on its boundary and it has been reached after 68 Kossaczky refinements from the initial subdivision of the cube into six tetrahedra. The mapping of the cube external nodes to the Armadillo surface produces a 3-D tangled mesh with 10871 inverted elements. The relocation of inner nodes by using volume parametrizations reduces the number of inverted tetrahedra to 413. We use the tetrahedral mesh optimization, presented in [10], such that the mesh is untangled in 3iterations. The mesh qualityis improved to a minimumvalueof 0.08 and an average qκ=0.71 after 6 smoothing iterations. In this case, we also note that the meccano technique generates a high quality tetrahedral mesh, Figure 13(b): only 4tetrahedron has a quality lower than 0.1,150 lower than 0.2and 1046 lower than 0.3. In Figure 14, we display cross sections of the cube and Armadillo meshes before and after the mesh optimization. The location of the cube is shown in Figure 14(a). The CPU time for constructing the final mesh of the Armadillo is approximately 54 seconds on a Dell precision 690, 2 Dual Core Xeon processor and 8Gb RAM memory. More precisely, the CPU time of each step of the meccano algorithm is: 5seconds for the subdivision of the initial surface triangulation into six patches, 1seconds for the Floater parametrization, 31 seconds for the Kossaczky recursive bisections, 2seconds for the external node mapping and inner node relocation, and 15 seconds for the mesh optimization. Finally, we summarize a comparison between our method and standard tetrahedral mesh generation techniques [37, 38, 39]. On the one hand, one of the most important contributions of the meccano method is the resulting volume parametrization of the solid. It can have interesting applications in solid modeling and numerical simulation. For example, the application of isogeometric analysis [2, 8, 3] can be easier. On the other hand, our volume meshes can be utilized in adaptive finite element applications by using the Kossaczky’s algorithm. The local refinement steps are very fast because the sequence of solid meshes is defined from the coarse mesh of the meccano, i.e., the dividing edge for tetrahedron bisection is known straightforward. In addition, we note that the minimum mesh quality is bounded during all the mesh adaptation process, because similar elements appear after three consecutive bisections. For a given solid surface triangulation, we have checked that a constrained Delaunay tetrahedralization [39] can be faster than our method, but the resulting meshes have lower quality. If the admissible minimum quality is increased, many points can be added and the number of tetrahedra increases to reach the objective. If only a conforming Delaunay tetrahedralization is desired, the number of tetrahedra can be much higher than in our method. A great advantage of our method is that the adaptive node distribution on the boundary and the inner of the solid is more structured, because the positions are fixed with a sequence of nested meshes. In the case of advancing front [37], the resulting mesh depends on the quality of the given surface triangulation. We have seen that the optimization of this triangulation, changing the mesh topology, can be too costly. 4.3 Isogeometric Model of the Bunny The volume parametrization presented in this paper has applications in other fields different from tetrahedral mesh generation. For example, it can be used to construct a volume T-mesh for isogeometric analysis [8]. The key lies in using the mapping, provided by the volume parametrization, to transform a T-mesh defined on the parametric domain (a unit cube in our case) into the physical domain. The T-mesh of the parametric domain is the parametric space in which the set of T-splines are defined [3]. The technique to construct a T-mesh starts by dividing the parametric cube in lower cubes by using an octree in such a way that each leaf of the octree is divided in eight children (eight cubes). The division continues until each terminal cube of the octree does not contain a node of the Kossaczky’s mesh in its inner. The octree defines a T-mesh in the parametric space that it is used to determine the local knot vector and the anchors of the T-splines following the description of [3]. Thus, the image of a point (u, v, w)in the parametric domain is given by S(u, v, w)= α∈A PαRα(u, v, w) where Rα(u, v, w)= Bα(u,v,w) β∈ABβ(u,v,w)is the T-spline blending function and Bαthe corresponding B-spline associated to the sαanchor. The index set Aruns over all the anchors of the T-mesh. The control points Pαare calculated by imposing the conditions S(sβ)=α∈APαRα(sβ)for all the anchors of the T-mesh. Here, we have also used the anchors as interpolation points. The image of each interpolation point sβin the physical space, S(sβ), is determined by the volume parametrization. We have applied our technique to the Stanford bunny. The original surface triangulation has been obtained from http://graphics.stanford.edu/data/3Dscanrep/, i.e. the Stanford Computer Graphics Laboratory. Figure 15(a) shows the cube tetrahedral mesh obtained by the meccano method. Figure 15(b) shows the parametric T-mesh. Figure 16(a) shows the tetrahedral mesh of the bunny constructed by the meccano method. Figure 16(b) shows the T-mesh of the Standford bunny generated by the application of the mapping S(u, v, w)to the parametric T-mesh. (a) (b) Figure 15: Meccano tetrahedral mesh (a) and corresponding parametric T-mesh (b)