# Abstract

Nowadays large part of the time needed to perform a numerical simulation is spent in preprocessing, especially in the geometry cleaning operations and mesh generation. Furthermore, these operations are not easy to automatize because they depend strongly on each geometrical model and they often need human interaction. Many of these operations are needed to obtain a watertight geometry. Even with a clean geometry, classical unstructured meshing methods (like Delaunay or Advancing Front based ones) present critical weak points like the need of a given quality in the boundary mesh or a relatively smooth size transition. These aspects decrease their robustness and imply an extra effort in order to reach the final mesh. Octree based meshers try to relax some of these requirements.

In the present work an octree based mesher for unstructured tetrahedra is presented. The proposed mesher ensures the mesh generation avoiding most of the geometry cleaning operations. It is based in the following steps: fit an octree onto the model, refine it following given criteria, apply a tetrahedra pattern to the octree cells and adapt the tetrahedra close to the contours in order to represent accurately the boundary shape. An important and innovative aspect of the proposed algorithm is it ensures the final mesh preserves the topology and the geometric features of the original model.

The method uses a Ray Casting based algorithm for the identification of the inner and outer parts of the volumes involved in the model. This technique allows the mesh generation of volumes even with non-watertight boundaries, and also opens the use of the mesher for immersed methods only applying slight modifications to the algorithm.

The main advantages of the presented mesher are: robustness, no need for watertight boundaries, independent on the contour mesh quality, preservation of geometrical features (corners and ridges), original geometric topology guaranteed, accurate representation of the contours, valid for immersed methods, and fast performance. A lot of time in the preprocessing part of the numerical simulation is saved thanks to the robustness of the mesher, which allows skipping most of the geometry cleaning operations.

A shared memory parallel implementation of the algorithm has been done. The effectiveness of the algorithm and its implementation has been verified by some validation examples.

# Resum

En l'actualitat gran part del temps emprat per córrer una simulació numerica esta dedicat al preprocés, especialment a les operacions de neteja de geometria i generació de malla. A més, aquestes operacions no són facils d'automatitzar degut a la seva forta dependencia del model geometric i sovint necessiten d'interacció humana. Moltes d'aquestes operacions són necessaries per aconseguir una definició topologicament hermetica de la geometria. Inclús amb una geometria neta, els metodes classics de mallat (com els basats en Delaunay o avancament frontal) presenten punts febles crítics com la necessitat d'una certa qualitat de les malles de contorn o una transició de mides relativament suau. Aquests aspectes disminueixen la seva robustesa i impliquen un esforc extra a l'hora d'obtenir la malla final. Els metodes de mallat basats en estructures octree relaxen alguns d'aquests requeriments.

En aquest treball es presenta un mallador basat en octree per tetraedres no estructurats. Un dels aspectes claus d'aquest mallador és que garanteix la generació de malla evitant moltes de les operacions de neteja de geometria. Es basa en els següents passos: encaixar un octree al model, refinar-lo seguint certs criteris, aplicar un patró de tetraedres a les cel${\textstyle \cdot }$les de l'octree i adaptar-los a les zones properes als contorns a fi i efecte de representar acuradament la forma del domini. Un aspecte important i innovador de l'algorisme proposat és que manté la topologia del model a la malla final i preserva les seves característiques geometriques.

El metode presentat utilitza un algorisme basat en la tecnica Ray Casting per la identificació de les parts interiors i exteriors dels volums del model. Aquesta tecnica permet la generació de malla de volums inclús amb contorns que no tanquen hermeticament, i també obre l'ús del mallador a metodes immersed aplicant només petites modificacions a l'algorisme.

Els principals avantatges del mallador presentat són: robustesa, no necessitat de definicions hermetiques dels contorns, independent de la qualitat de la malla de contorn, preservació de característiques geometriques (cantonades i arestes abruptes), topologia original de la geometria garantida, representació precisa dels contorns, valid per metodes immersed i rapid rendiment. L'ús del mallador estalvia molt de temps en la part del preprocés de la simulació numerica gracies a la seva robustesa que permet obviar la majoria d'operacions de neteja de geometria.

S'ha dut a terme una implementació paral${\textstyle \cdot }$lela amb memoria compartida de l'algorisme. L'efectivitat del mateix i la seva implementació ha estat verificada mitjancant exemples de validació.

# 1 Introduction

Numerical simulations try to reproduce virtually a physical behavior by solving given equations in a specific domain. They are nowadays essential to understand some complex physical problems in scientific and engineering field. Although experimental setups can be build to study the specific behavior of a given phenomena, sometimes it is hard for these experiments (or even impossible depending on the scale of the tackled problem) to represent it accurately. The increasing advances in terms of computer science technology allow to treat larger and larger problems virtually (in the computer), so each time more and more numerical methods have been developed in the scientific field in order to capture the physics of complex problems.

Together with these developments, the adoption of numerical simulations in industrial processes has became a reality, as it can save a lot of time and effort when evaluating possible solutions for a given problem. The use of numerical simulation tools by the industry requires a software with very high level of robustness, efficiency and performance.

The process to run a numerical simulation involves three main parts: pre-processing, calculation and post-processing. The pre-processing part includes all the operations needed to define and discretize the geometrical domain, assign the required data to it so that the solver can solve the corresponding equations representing a given physical problem (in the calculation part). The post-processing part tries to analyze and visualize the results from the solver in a smart way so that they can be correctly interpreted. This monography is focused on the pre-processing part of the numerical simulations, and specifically on the discretization of the geometrical domain.

The pre-processing operations can be summarized as follows:

• Geometrical definition of the domain. This part of the process is not always performed by the person responsible to run the simulation, and often a third part (a designer) creates the CAD definition of the geometry. Most of the times, this geometry creation does not take into account the requirements of numerical simulation, so a geometry cleaning process is needed in order to be able to perform the discretization of the domain.
• Discretization of the domain. Depending on the nature of the numerical method used for the simulation, different kinds of discretization of the geometrical domain are required. Meshless methods use a collection of points (nodes) as a discretization of the domain. Discrete Element Methods (DEM) use a collection of topologically unconnected objects for this purpose. In many occasions (like in the Finite Element Method (FEM) [1,2], Finite Volumes (FV) or Finite Differences (FD)) the result of the discretization is a mesh: a collection of polygons (in 2D) or polyhedra (in 3D) occupying the space where the domain is. This work focuses on this type of discretization, so in this document the discretization part of the pre-processing operations will be referred as mesh generation.
• Assignation of the data needed to run the simulation. This data can be of different nature depending on the kind of numerical simulation to be run, but most of the times includes the material properties of the different parts of the domain, the initial and boundary conditions for the equations to be solved, and some general parameters for the simulation.

The presented work proposes a new algorithm for mesh generation: a mesh generator (or simply a mesher). For designing a mesh generator, the basic requirements to be covered must be very clear. Missing the right requirements for the numerical simulation may lead to several limitations in the use of the mesher. There are three basic requirements to be covered by any mesh generator:

• Input data requirements: these requirements fix the characteristics of the input data able to be processed by the mesher. In this sense the input data is understood as the geometrical definition of the domain into which the simulation will be run. Typical requirements of this nature are: allow only a given mathematical definition of the geometry (NURBS, mesh, etc.), a specific topology of the input geometrical entities (watertight geometries, contact entities, etc.) or sizes of input entities inside a given range, among others. Considering the large number of CAD systems available and their different ways for geometrical definition of the domain, these kinds of requirements are in practice really important. Many times the person in charge of the numerical simulation may have no control on the geometrical definition. The situation where the mesher is very restrictive with the quality or the topology of the geometrical definition of the domain, often leads to a considerable increase of the effort in the geometry cleaning operations before the meshing itself. These requirements are not theoretically taken into account to evaluate if a mesher is suitable to generate meshes for given kinds of simulations or not, but in practice they can limit the use of the mesher as the effort dedicated to geometry cleaning operations can grow exponentially when the input geometries become more complex.
• Final mesh requirements: these requirements focus on the specific characteristics of the mesh required by the simulation to be run successfully. Typical requirements of this nature are: uniform elements sizes distribution, control of the mesh sizes in different zones of the domain or precise representation of the contours among others. Each kind of numerical simulation (structural dynamics, computational fluid mechanics (CFD), electromagnetism, etc.) has specific requirements for the mesh to be used. If a numerical simulation has one requirement of this nature for the final mesh and it is not covered by a mesher, this mesh generator is not suitable for the simulation.
• Mesher requirements: these requirements deal with the mesher behavior itself. The mesher can cover the two kinds of requirements presented before in several ways (with different implementations). The use of the mesher in practice often fixes some requirements in terms of speed, usability or memory availability which are also important.

Often, a mesh generator is focused in one type of mesh. Different kinds of meshes can be identified depending on their nature:

• 1D, 2D or 3D; depending on the dimensions of the space used to define the geometry.
• Line, surface or volume meshes; depending on the hierarchy of the geometrical entity to be meshed.
• Structured, semi-structured, unstructured or cartesian; depending on the nature of the topology of the mesh.
• Isotropic or anisotropic; depending on the aspect ratio of the elements to be generated following specific directions.
• Embedded or body-fitted; depending on how they fit or not the boundaries of the domain to be meshed. Body-fitted meshes (used, for instance, in FEM or FV simulations) match perfectly the boundaries of the domain, as embedded ones (used, for instance, in FD, immersed or embedded methods) have faces crossing that boundaries.
• Conformal or non-conformal; depending on the continuity between neighbor elements across edges or faces. In conformal meshes edges and faces match perfectly between neighbor elements, as oposed to non-conformal ones.
• Depending on the element type to be used and the quadratic type of the element (linear or different degrees of quadratic types). Although there are methods using arbitrary polyhedral shapes for the mesh elements [3], it is common to use the same simple polyhedral element for all the domain (triangle, quadrilateral, tetrahedral, etc.).

In this monography 3D unstructured isotropic conformal tetrahedral meshes are considered. Both embedded and body-fitted cases will be covered.

## 1.1 Motivation

In industrial simulations, the pre-processing operations represent the most time consuming part of the whole process. Among the pre-processing operations, the geometry cleaning and mesh generation parts are the ones which consume more time, due to their specific characteristics:

• The geometry cleaning process implies all the operations needed to allow the mesher the correct discretization of the domain. It often needs much human interaction, as it is really difficult to automatize it following general criteria for all kinds of simulations. Most of the meshers need watertight geometries to be able to generate the mesh. In this context, the concept of watertight (for volumes) is used to define a closed collection of surfaces sharing their contour lines. Make watertight an input geometry implies to avoid surface overlapping and gaps in the geometrical definition of the volumes. For many complex 3D geometries (like the one shown in Figure 1), to reach this goal requires a huge effort. Furthermore, classical meshers are not totally independent on the mathematical definition of the geometry. Often the characteristic size of the geometrical entities defining the boundary must be similar to the final mesh size desired for the simulation. This aspect enforces a modification on the definition of the input boundaries depending on the final mesh desired, although it may represent well the shape of the domain.
• Once the geometry is clean enough to generate the mesh, the mesh generation itself can consume a considerable amount of time. Furthermore, the mesh generation process often involves an iterative loop where the person in charge of the simulation defines the desired mesh sizes in the different parts of the domain, and several meshes have to be generated, as the mesher parameters are tuned in order to get an optimum mesh. It also has to be taken into account that nowadays the volume unstructured meshers are not fully robust in the sense that, when dealing with complex geometries, it is not easy to tune its parameters and prepare the input geometry to generate the mesh successfully: often the user must try several configurations (geometry and parameters) until a mesh is generated. This relativizes the speed in the mesh generation, as the mesher itself can be very fast, but the user may need several tries to get a suitable mesh for the simulation.
 Figure 1: Example of a 3D input boundary (represented with a mesh) where the quality of the triangles is very bad.

Much effort has to be spent in order to generate a volume mesh of a complex geometry considering the existing meshing algorithms. Some of them are really fast and robust, but they require several geometry cleaning operations. If the time needed for them is added to the total meshing time, they became not so fast in practice. Contribute to reduce this extra effort is the motivation of this monography.

## 1.2 Objectives

The main objective of this monography is to develop an algorithm for isotropic unstructured volume mesh generation robust enough to be able to generate a mesh from non-cleaned input geometries, using as less input data as possible. This will lead to a drastic reduction of the time consumed in the pre-processing part, and will overcome the actual bottleneck in the whole simulation process.

The mesh generator must be flexible enough in terms of mesh adaptation to the solver requirements considering specific input data. The idea in this work is that the algorithm should be able to generate always a mesh from a given geometrical domain, almost without any specific meshing property assigned to it.

The meshing algorithm presented in this monography has been designed with the following objectives in mind:

• It should be prepared to use the more common geometrical definitions as input data reducing as much as possible its preparation (CAD cleaning operations).
• It must be able to run on a simple PC, but also it should be able to take profit from more complex computer architectures, which enable the use of parallel processing.
• It should be useful for a wide range of application fields. The main characteristics for the mesh to be generated are taken from the requirements for the FEM, but they are general enough to be applied to other numerical methods.

These characteristics lead to a set of requirements to be covered by the meshing algorithm developed in this monography. The main ones are: robust and fast mesh generation, and ability to mesh from non-watertight geometries. The explanation of all the requirements covered by the new mesher is detailed below in Sections 1.2.1, 1.2.2, and 1.2.3.

### 1.2.1 Mesher requirements

In this section, the requirements of the mesher developed in the monography regarding its behavior are detailed:

• Robustness. This is one of the key objectives and it can also be seen as a requirement for the input data, as this robustness refers to the capability of generating the mesh independently from the input data quality. Some of the existing unstructured volume meshers accomplish all (or almost all) the requirements presented in this work, but they require a very specific tuning of their parameters for successfully generating the mesh. Independently of the sizes assigned by the user, the quality of the input geometrical data or the general parameters chosen, the mesher should be able to generate always a mesh suitable for the simulation. The idea of the new algorithm developed in this monography is that it should generate always a mesh accomplishing the requirements defined in this section without a special training in the use of the mesher.
• Fast mesh generation. One of the key motivations of this work is to reduce the time consumed in the pre-processing operations for a numerical simulation. For this purpose the mesher has to be fast. Moroves, it has to be naturally fast, in the sense that it should imply fast and simple operations. Furthermore, its implementation should take profit from parallel computing strategies.
• Prepared to generate really large meshes. As the computer technology is evolving very fast, each time the solvers are capable to solve larger and larger algebraic systems. Thinking about the requirements from different simulations fields for the next years, the mesher should be able to generate billions of tetrahedral elements without problems. This characteristic is really important when deciding some implementation aspects of the algorithm.

### 1.2.2 Input data requirements

This section focuses on the requirements to be covered by the mesher concerning the input data:

• Accept non-cleaned input geometry. This is another key point of the mesher. The idea is that it should be able to generate the mesh of the model although its geometrical contour comes from a non-cleaned input geometry. In this context, a non-cleaned input geometry is defined as the one containing:
• Non-watertight volumes. This means volumes with contours presenting gaps or overlapping entities.
• Non-coherently oriented contour entities. This means that the contour entities of a volume are not needed to be oriented (all of them) towards the same part of the volume (inner or outer).
• Low quality definition. The geometrical contour can be defined by geometrical entities (like NURBS) or by a surface mesh. In this context, a geometrical contour with low quality definition means a highly distorted parametrization of the surfaces (in case of geometrical representation), or a mesh containing low quality elements (in case of surface mesh representation). Examples of low quality input geometries and no water-tight domains are shown in Figures 1 and 2. It has to be noted that the algorithm proposed in this document is for 3D volume meshing. However, for sake of simplicity, 2D examples are shown in some figures.
 Figure 2: Example of a 2D non watertight input boundary with gaps and overlapping entities.
• Accept input data in geometry and mesh format. The natural input for a mesher is the geometrical definition of the contours of the model. This geometrical definition can be represented in several ways, but the most common ones are CAD or mesh entities. In this document, when talking about CAD entities, the NURBS surfaces and curves [4] (trimmed or not) will be considered, as they are the most general mathematical representation able to represent all the geometrical shapes.

Although the mesher should generate a mesh from a non-cleaned input geometry, some minimum criteria concerning topology or shape definition may be needed.

While mesh entities are simpler to be defined, CAD entities are a more precise way of defining a geometry. Meshes often loose continuity in the geometrical definition depending on the smoothness or curvature of the shape to be represented. It also has to be considered that, almost always, the mesh presents a chordal error compared to the smooth original geometry (the chordal error is the distance between a point on the mesh and the original smooth shape the mesh is trying to represent). In the Figure 3 a graphical representation of this chordal error is shown.

 Figure 3: Graphical representation of the chordal error in a 2D case. The black thick line represents the smooth geometry, and the gray thin line represents its mesh. The distance ${\displaystyle d}$ is the chordal error of the element ${\displaystyle e}$

The simpler definition and treatment of mesh entities has made them the most common used as input for the existing volume meshers. There are cases where a surface (or a line) cannot be meshed because of its bad parametrization. In these cases, a volume mesher requiring a mesh as input cannot be used. This is the reason why the proposed volume meshing algorithm should be prepared to get as input data mesh entities, as well as NURBS surfaces and curves. Actually, it should be prepared not only to get these kinds of entities as an input, but also to work with them in the geometrical operations needed during the whole meshing process. The methodology presented in this monography is prepared to work either with CAD and mesh entities. Hence, the general form of geometrical entities will be used to refer the input entities defining the contours of the domain. Only if some specific operation is needed for just one of the representations it will be specified if the geometrical entity is a CAD or a mesh one. Following this duality, in this document a single term will be used to refer the different natures of geometrical entities:

• Surface entities for surface mesh elements (especially triangles) and geometrical surfaces (especially NURBS surfaces and patches of connected NURBS surfaces).
• Line entities for line mesh elements and geometrical curves (especially NURBS lines and groups of connected NURBS lines).
• Point entities for nodes of a mesh and geometrical points.

### 1.2.3 Final mesh requirements

The requirements to be covered by the mesher regarding the final mesh generated are:

• Maintain volume topology. The topology of the initial domain to be meshed must be preserved in the final mesh. For most meshing algorithms this requirement is automatically satisfied, but it is not automatically guaranteed in all meshing techniques. Considering the final mesh, it has to be possible to identify the volume each tetrahedron belongs to. Then, if getting all the tetrahedra belonging to a given volume, the concept of the mesh of that volume can be defined. The maintenance of the volume topology of the initial domain means that there must be a mesh for each of the initial volumes, and each volume mesh must represent the same relationship between the others. For instance, if two volumes are connected by a line, the meshes of these two volumes must be neighbors sharing line elements. Furthermore, the mesh of a volume itself must represent the topology of the volume (for instance, if the volume has a hole, the mesh must have a hole). In Figure 4 two examples are shown with meshes preserving or not the topology of the original model.
 (a) (b) (c)
Figure 4: Topology preservation between two volumes and their meshes. (a) Two volumes in contact sharing a line. (b) Meshes of the volumes shown in (a) sharing the line elements corresponding to the mesh of the line shared by both volumes: the topology is preserved. (c) Meshes of the volumes shown in (a) not sharing line elements: the topology of the original model is not preserved.

A special case which evidences the importance of maintaining the initial topology is the situation where the domain has very thin parts representing relevant details of it. A 2D case of this kind is shown in Figure 5, where a thin channel-like part can be identified in a surface (Figure 5(a)). Independently on the mesh desired size required by the simulation, the mesher must generate elements small enough not to close this channel-like zone, as the mesh shown in Figure 5(b). The mesh depicted in Figure 5(c) is not acceptable, as it does not preserve the topology of the domain.

 (a) (b) (c)
Figure 5: Example of surface mesh preserving or not the topology of the initial domain. (a) 2D geometrical domain formed by two surfaces (colored as blue and gray). (b) Mesh of the surfaces preserving the topology of the initial domain. (c) Mesh of the surfaces not preserving the topology of the initial domain.
• Maintain a representative lines and surfaces topology. Often it is not required to preserve the topology of all the surface and line entities of the domain to be meshed. As the proposed algorithm is a volume mesher, the surface or line elements in the final mesh are of interest just in order to assign properties to them for the simulation. These properties can be of different nature: material properties, boundary conditions, etc. At the end, the mesher should maintain the topology of some patches of lines and surfaces on user demand, or following some automatic criteria.

Note that this requirement is not as hard as the typical constrained condition in the boundaries of the domain. Especially when the input data for a mesher is a surface mesh, it is common to require the mesher to be constrained at the boundary. This means that the contour mesh of the final tetrahedral mesh (the triangles representing the skin of the generated tetrahedra) must be topologically identical to the input surface mesh defining the contours of the domain.

 (a) (b) (c)
Figure 6: (a) Contour mesh of a volume with a set of triangles highlighted in red representing the region ${\displaystyle A}$. (b) A partially constrained tetrahedra mesh generated from the contour mesh shown in (a) where the highlighted set of tetrahedra faces corresponds to the region ${\displaystyle A}$. (c) A not constrained tetrahedra mesh of the contour mesh shown in (a) (it cannot be identified a set of tetrahedra faces representing the region ${\displaystyle A}$).

In Figure 6(a) an example of the contour of a volume is shown, where a set of surface entities (in this case triangles)are colored in red. Let us call ${\textstyle A}$ the region represented by those triangles. If the simulation requires a set of triangles representing the region ${\textstyle A}$, the mesh generator should provide with a tetrahedral mesh whose skin should have a set of triangles representing that region, but it is not needed for those triangles to be identical as the ones in the input geometry. A tetrahedral mesh accomplishing this criterion is shown in Figure 6(b). A totally unconstrained mesh is depicted in Figure 6(c), where it cannot be identified a set of triangles corresponding to the region ${\textstyle A}$. It is obvious that the mesh in Figure 6(b) is not constrained with the input boundary shown in Figure 6(a), as the triangles representing the region ${\textstyle A}$ are different. However, it can be seen that the requirement of maintaining a representative topology of the surface and line entities of the input data allows a bijective relationship between groups of surface and line entities in the input data and surface and line mesh elements in the final mesh. This requirement is enough to allow an automatic assignment of data from the input boundary geometry to the final mesh entities.

• Generate automatically the contour mesh of the volumes. Although the mesher is thought to generate volume elements (tetrahedra), it should give as an output also the triangular elements corresponding to the contours of the volumes, and the line elements corresponding to the boundaries of certain patches of triangles. This surface and line elements correspond to faces and edges of the final tetrahedral mesh.
• Preserve geometrical features. This is automatically achieved by the boundary constrained meshers, but it is not guaranteed at all by the types of meshers. This requirement is crucial for the final mesh to represent the shape of the domain precisely. For several kinds of numerical simulations, the presence of sharp edges and corners in the domain affects drastically the results, as they often govern the physical behavior of the process to be simulated. The importance of preserving sharp edges and corners to represent some domains can be appreciated in Figure 7. Part of the contours of a volume are shown in Figure 7(a). The tetrahedral mesh of that volume generated without preserving the sharp edges of the input geometry is shown in Figure 7(b).
 (a) (b) Figure 7: (a) View of a part of the contours of a mechanical piece. (b) Tetrahedral mesh of the mechanical part generated without preserving the sharp edges of the input geometry.

In terms of the mesher, this requirement means that it should be partially constrained to some line entities from the input geometrical data, and must include fixed nodes in the final mesh (corresponding to specific point entities in the input data). In this context, partially constrained means that if some specific line elements in the input geometrical data must be preserved, a collection of edges from the final mesh should follow the path of those line elements. This is not as restrictive as totally constrained condition.

As it will be explained later on (Section  3.1), the strategy to cover this requirement will be also useful to reach the requirement of maintaining a representative lines and surfaces topology.

• Allow to skip given details of the domain. This requirement is somehow complementary to the previous one. What it means is that the mesher must be not constrained in some regions. The definition of the geometrical domains includes often very thin or small entities in given parts. The reason for the size of these entities may be the size of the part of the domain represented itself, but the presence of these small entities often responds to the application of a given tangency criteria for the geometrical definition, or they are just the result of some geometrical operation previously done (intersections, Boolean operations, etc.). The basic idea of this requirement is that the sizes of the geometrical entities used to define the contours of the domain are not necessary related to the mesh size needed for the simulation. In the Figure 8, an example is shown where the triangles generated from a patch of surfaces skipping the inner lines between them are much larger than the size of some of the surfaces of the patch. Clearly the simulation should not need such thin triangles in the regions where the surfaces are so thin.
 (a) (b) Figure 8: (a) Zoom of a patch of surfaces representing a mechanical part. (b) Triangle mesh of the patch of surfaces shown in (a).

Even when the sizes of the input entities are related to the representative size of a specific shape to be represented, the user may desire a bigger mesh size to skip the representation of given details because they are not of interest for the simulation. Instead of forcing the user to modify the input geometry which defines the domain, this requirement makes the mesher skip the detail within the mesh generation process.

The accomplishment of this requirement will let the user to define a desired mesh size distribution in the final mesh just because of the simulation requirements, and not because of the initial way of defining the geometry. In terms of the mesher, this requirement means that it should be able to skip some of the line and point entities from the input geometrical definition of the domain. The entities to be skipped should be defined automatically following some given criterion, or by the user.

• Precise representation of the contour geometry. This requirement focuses on the need of the nodes of the contours of the final mesh to lay exactly on the geometrical entities defining the contours of the domain. Some mesh generators only approximate the shape of the domain to be meshed, or represent its contours in a staircase manner (specially the meshers based on octree or bin structures). Several numerical methods (specially the FEM) need a precise representation of the contour, as its shape and smoothness may affect drastically the result of the simulation. This is the main reason for this requirement.
• Final elements quality ensured. The acceptable quality of the mesh elements is a relative parameter, as the different methods used for the numerical simulation are applicable for different ranges of element qualities [5]. Usually, most numerical methods perform better as the shape of the element is closer to the corresponding regular polyhedron, however this cannot be extended to all numerical simulations. This requirement will try to reach a minimum dihedral angle in the final tetrahedra.
• Allow a 3D spatial size assignment. The use of a non-uniform mesh is crucial for some simulations. In some specific regions of the domain the size of the elements may have an upper bound, but if that size is used in the whole domain the final mesh should have too many elements. Typically, in the regions where there are higher gradients of the variable under study, the simulation needs a finer mesh than in other regions. This requirement is also important in the accurate representation of the shape of the domain, as the curvature of the domain geometry is strongly related with the mesh size needed to accomplish a given chordal error criteria. Just for geometrical purposes, a refined mesh in some regions would be needed.
• Control of sizes transitions. As the final mesh can be non-uniform size distributed, the size of the elements should vary from a region of the domain with a desired size to another. Some numerical simulations do not only require a specific mesh size in some specific zones, but also a given growing law between neighbor elements. The mesher should let the user to control how the element sizes grow (or decrease) from a region of the domain, depending on the distance from the elements to that region.
• Applicable for immersed methods. A family of methods for numerical simulations are the so-called immersed methods and embedded methods [6,7]. The main characteristic of these methods is that they do not require a body-fitted mesh to represent the domain. The shape is implicitly introduced by distance from the nodes to the boundaries and a casting process (determine whether a node is inside or outside the domain). In some sense, generating a mesh for a simulation using an immersed method is less restrictive as the geometrical features should not be preserved, but in contrast, the final mesh must have extra information attached. This information is the distance from each node of the mesh to the boundary of the domain.

### 1.2.4 Surface meshing

Although the main objectives of the monography are related to volume meshing, there is a secondary objective for the new mesher, which is to be able to mesh 3D surfaces and lines which are not boundary of any volume accomplishing the requirements defined above. This capability would give several advantages to the method:

• Generate conformal meshes of volumes containing inner surfaces or lines. These are entities totally inside the volume, but not belonging topologically to the boundaries of the volume. An example of this case could be the representation of a sail (3D surface) into a control volume in order to perform a CFD simulation. The sail is totally inner to the control volume, and the mesh of it must be conformal with the tetrahedra representing the control volume. A 2D example of this case is shown in Figure 9. In this example a surface (bounded by solid black lines) has an inner line (dotted black curve), and the mesh of the surface (in gray) is conformal with this line.
 Figure 9: 2D example of a surface mesh (gray lines) conformal with an inner line of the surface (dotted line).
• Mesh in a conformal way a model with volumes and surfaces or lines connected to them. A case of two volumes connected by a surface is depicted in Figure 10. The mesh of the volumes and the surface must be conformal in the contact lines between them. This case can be also meshed by generating separately the volume meshes with the octree mesher, and afterwards mesh the surface using other meshing algorithms. This would not take profit on the advantages of meshing the whole model at once.
 (a) (b) Figure 10: (a) Example of a model containing two volumes (in blue) and a surface (in grey) connecting them. (b) A conformal mesh of the model considering the volumes and the surface mesh.
• Mesh patches of surface entities together preserving features, but skipping the inner line entities between them if needed. The entities to be skipped could be detected automatically following a given smoothness criteria, or under user demand. Figure 8 shows an example of triangle mesh representing a collection of surface entities, but skipping the inner line entities.
• If the patches of surfaces are not topologically connected, the proposed method could act as a CAD cleaning tool, as it could extract a continuous manifold surface mesh from a collection of surface entities with gaps or overlaps. The extraction of surface skin meshes from a volume discretization has been used in [8] and [9], and it is interesting for CAD cleaning, as well as for boolean operations involving non-watertight geometries. The preservation of geometrical features is essential for these purposes.

## 1.3 Structure of the monography

The monography is structured in the following chapters:

• Chapter 1: Introduction. This chapter focuses on introducing the context of the presented work. Its motivation is explained (Section 1.1), and its objectives are detailed (Section 1.2).
• Chapter 2: State of the art on mesh generation. In this chapter the different mesh types are presented (Section 2.1) as well as the state of the art of different meshing techniques (Sections 2.2, 2.3, 2.4 and 2.5). The advantages and drawbacks of each meshing method is studied, and the proposed solution for the meshing algorithm is justified (Section 2.6).
• Chapter 3: Basic concepts of the new mesher. Several concepts and some auxiliary algorithms are defined in order to allow a better understanding of the meshing algorithm presented. Its main parts are:
• Definition of the input data for the presented meshing algorithm and a proposal of how the mesher can interact with CAD data, integrated within a pre-processing system (Section 3.1).
• Octree structure characteristics (Sections 3.2 and 3.3). The definition and main characteristics of the octree structure (which plays a key role in the presented work) are presented. Special interest is given to its specific properties for the meshing process. The notation is defined which will be referred to in the following chapters.
• Geometrical intersections (Section 3.4). As the geometrical intersections are crucial for some of the algorithms presented, this section analyzes them deeply and studies the pathological configurations that may occur.
• Chapter 4: Coloring algorithm. The coloring process determines where a point in space is topologically: it determines if it is outside the domain, inside a given volume, or onto an interface between volumes. This is one of the main operations involved in the presented mesher and can be understood as an auxiliary algorithm for the meshing itself, so a whole chapter is devoted to it. In this chapter different approaches to solve the coloring problem are studied (Section 4.1), and a new algorithm based in the ray casting technique is exposed (Section 4.2). The implementation of this algorithm is presented in Section 4.3.
• Chapter 5: Octree based mesher. This chapter focuses on the new meshing algorithm itself. After an introduction, where its main idea is highlighted (Section 5.1), the algorithm for embedded (Section 5.2) and body-fitted (Section 5.3) meshes is detailed.
• Chapter 6: Implementation aspects. In this chapter all the implementation details for the meshing algorithm developed are explained. After some general aspects of the implementation (Section 6.1), the implementation of the octree structure is analyzed (Section 6.2). Sections 6.3, 6.4, and 6.5 focus on the implementation of specific parts of the meshing algorithm itself. Some considerations on the parallel implementation of the mesher are pointed out in Section 6.6. A list of the values of the parameters relevant for the mesher used in the presented implementation of the algorithm is detailed in Section 6.7.
• Chapter 7: Examples. In this chapter the results of some validation examples highlighting specific characteristics of the presented mesher are shown (Section 7.1). Sections 7.2 and 7.3 analyze the results of the meshing algorithm applied to real complex geometries under different configurations.
• Chapter 8: Conclusions and future research lines. This chapter lists the conclusions of the presented work (Section 8.1) and proposes some futures lines of research (Section 8.2).
• Appendix A: Profiling tables and complete data of examples. In this appendix, tables with all the profiling data concerning times and memory for the whole configurations used in the examples run to validate the meshing algorithm are provided.

It is highlighted that, although the meshing algorithms presented are directly applied to volume meshing, most of the geometrical operations are applicable (with slight modifications) to surface meshing in 2D cases. Taking into account that some concepts are much more clear to understand from a 2D scheme (especially when dealing with geometry), in this monography 2D examples are used sometimes to illustrate some of the concepts explained.

# 2 State of the art on mesh generation

In this monography isotropic volume meshers are considered. Isotropic meshers can be defined as the ones trying to generate elements as much regular as possible, understanding a regular element as the one whose edges have the same length.

## 2.1 Mesh types

There are two main families of meshers depending on the kind of mesh they generate: structured and unstructured [10]. Actually, a third kind of meshers can be classified as semi-structured ones. A structured mesh is defined as a mesh which all inner nodes have the same degree (the degree of a node is the number of elements owning it), while the nodes of an unstructured mesh have different degrees. Semi-structured meshes can only be applied to topologically prismatic geometries, and they basically repeat the structure of an unstructured mesh (in the tops of the prismatic shape) in different layers following the structured direction. An example of this kinds of mesh is shown in Figure 11.

Unstructured meshers  [11,5,12] can be divided in three main families: advancing front, Delaunay and space decomposition methods.

In the following sections, the main characteristics of these methods as well as their main advantages and drawbacks are detailed, focusing on the requirements defined in section 1.2. The aim of this chapter is to highlight which of those requirements are covered by each meshing method, so the algorithms are not deeply detailed and only their main characteristics are pointed out.

 (a) (b) (c)
Figure 11: Examples of different types of mesh(a) Structured triangle mesh. (b) Unstructured triangle mesh. (b) Semi-structured prism mesh.

## 2.2 Structured meshers

Structured and semi-structured meshers often get as input data the position of the nodes in the contours of the domain, and generate the inner nodes positions from a given interpolation  [10,12,4]. The quality of the final mesh obtained is directly related to the kind of interpolation used and the degree of distortion of the contours of the domain, so a minimum level of element quality cannot be guaranteed for arbitrary domains. However, for good shaped volumes and uniform sizes distributions these methods provide with very good quality meshes.

The main advantages of structured meshers are:

• Fast and robust. As they are based on a given predefined interpolation, they are really fast and robust.
• Parallelizable. Apart from being fast, these kind of meshers are very parallelizable.

On the other hand, these meshers have some important drawbacks:

• Need of specific topology for the input data. The main problem of this kind of meshers remains in a requirement for the input data: they need a specific topology for the geometrical definition of the contours of the domain. As an example, to generate a structured mesh of hexahedra, the input geometry must be topologically an hexahedra. This means it has to be a volume with 6 contour surfaces, and each one of their must have 4 contour lines. This requirement makes impossible to use this kind of meshers for arbitrary geometries with complex topology. A family of methods has been proposed that are able to skip this problem  [13]. They are basically based in decomposing the original domain to be meshed into different parts which do accomplish with the topology demanded by the mesher. Then, each of these parts can be meshed in a structured manner, so they are often referred as structured by blocks or multi-block meshers. Several implementations of these methods have been carried out [11,10], but in practice they are meaningful for very specific domains only. Depending on the topology required by the mesher, some geometries are impossible to be decomposed in such this way, and even when this decomposition is possible, from complex original geometries is not obvious to generate this decomposition automatically. This often implies an extra pre-processing operation before using the mesher, which is the splitting (manually) of the domain in different parts.
 Figure 12: Example of a 2D structured quadrilateral mesh with a non-uniform size distribution. A high level of element distortion can be appreciated.
• High distortion in non-uniform meshes. As the nodes of structured meshes have the same degree, a predefined topology of the mesh is forced. It generally tends to generate uniform meshes. Structured meshers can relatively distort the elements in some region in order to provide with a non-uniform mesh, but the fixed topology of the mesh does not allow a sharp sizes transition without an excessive distortion in the elements. A 2D example illustrating this problem is depicted in Figure 12.

The advancing front method  [14,15,16,17] is a common technique for generating unstructured meshes. It gets a closed and oriented mesh of the boundary of the domain as input and mesh its inner part. The surface elements of this mesh are the ones in the active front, and the algorithm can be summarized with the following points:

• From each element in the front (face) a new volume element is generated. An optimal position of the node needed to build the element is obtained and the nodes and faces of the active front close to it are considered. The new element is build using a new node in the optimal position or an existing close node, depending on the resulting configuration: the new element should not intersect any one of the close already existing faces.
• The active front is updated adding the new faces created, and removing the ones already shared by two elements.
• These two steps are repeated until there are no faces in the active front. Then all the domain has been filled with volume elements.

Although the advancing front is mainly used for generation of tetrahedral (or triangle in 2D) elements, some adaptations of the method have been done to generate other types of elements [18], and even for generating particles for DEM simulations [19]. The present work is focused in tetrahedra mesh generation. There are several possible implementations of the advancing front method [11] depending on the way of creating the new elements from a face of the front, the way of considering the mesh desired size in the inner part of the domain, or the order in which the faces of the front are processed, among others. In special, much work have been done in order to improve the efficiency in evaluating the desired mesh size in a specific region of the domain, and control a smooth size transition in the final mesh. This is one of the strong points of advancing front techniques, as the creation of each element can fit very well the desired mesh requirements. Different approaches use a background grid to set the mesh desired size  [16,15,20,21], or sources from which the desired mesh size vary following a given function [22], which provides with a very smooth transition between element sizes. A combination of different methods can be used in order to improve accuracy and reach an efficient implementation of the method [23]. However, independent of the implementation, one can identify some general advantages and disadvantages of advancing front based techniques. The main advantages are:

• Good quality meshes. In the advancing front method, each new element is created connecting a face of the active front with a new node or an existing one. This choice is done based on an optimal position for this new node in order to get the best element (in terms of quality) taking into account the desired size of the mesh in that region. This methodology for creating the elements gives high quality elements in the final mesh, even before any smoothing process.
• Automatic preservation of topology. The advancing front method is naturally constrained, as each element of the surface mesh of the contours of the domain to be meshed will become a face of the final tetrahedra.
• Good control on the size transitions. As explained before, each new element is created from an existing face and an optimal position of the node to create the element from the face. Considering the desired size for the final mesh in the region where the face is, and the size of the face itself, the size transitions from one part of the mesh to another is totally controlled.
• Parallelizable. The creation of each new element using the advancing front method is clearly local. The geometrical and topological checks needed to create the new element only take into account the nodes and faces inside a specific radius of influence of the face the new element is created from. This aspect makes totally independent the creation of an element from the creation of another element which is far enough. Although it is not obvious to implement an efficient parallel version of the advancing front method, at least the element generation method is local enough to make the implementation affordable. [24] proposed a first octree subdivision of the domain, and then apply the advancing front technique in each of the cells of the octree, seaming afterwards the interface parts between cells.

The main drawbacks of the advancing front method are:

• Strong dependency on the contour mesh quality. Considering that the surface elements of the initial contour mesh will be faces of the final volume elements, the quality of these contour elements is totally related with the quality of the tetrahedra in contact with them: if a triangle of the contour of the domain has a very small angle in some of its nodes, the resulting tetrahedron from it will have a very small dihedral angle. This implies a great effort in the mesh generation of the contours of the volumes to be meshed. Often this effort takes most of the part of time in the pre-processing operations of the whole simulation when using the advancing front technique.
• Requires watertight input geometry. The reason for this is that the initial surface mesh representing the boundary of the domain will be the skin of the tetrahedra in the final mesh.
• Not very robust. Theoretically, the advancing front is a robust method. However, the locality of the element creation process and the checks needed to evaluate if the new element created intersects or not the active front require some tolerances. These tolerances must be tuned precisely in order not only to avoid the crossing of the active front, but also to avoid different parts of the front to get very close. This situation can lead to geometrical configurations of the active front which make almost impossible the creation of a new element. The tuning of these tolerances and other parameters can be done as general as possible, but as the active front evolves during the meshing process, these geometrical configurations are not under control, and cannot be predicted. This aspect, together with the hard dependency on the contour mesh quality makes the advancing front method not so robust in practice.
• Not very fast. Every candidate new element to be created needs to be checked in order to verify if the active front is intersected by some of its faces and edges. This involves the use of several geometrical intersections operations in the whole meshing process. This aspect, together with other characteristics of the advancing front method, makes it not very fast in comparison with other meshing algorithms.

## 2.4 Delaunay method

A Delaunay mesh is defined as a mesh which elements accomplish the Delaunay condition: the circumcircle (in 2D case) or the circumscribed sphere (in 3D case) of any element has no node from the mesh inside  [10]. Given a cloud of points, a Delaunay triangulation can always be created from their Dirichlet tesselation connecting them with a set of triangles (in 2D) or tetrahedra (in 3D) without adding any extra node [10].

The Delaunay meshing methods [12,25,5,11] depart from the contour of the domain and generate the Delaunay triangulation of its nodes. This mesh is the convex hull of the domain to be meshed. Although it is already a mesh, its elements may have a low quality, or may not fit with the desired mesh size in that region of the domain: these are the bad elements. The following strategy is applied recursively to all the bad elements:

• A node is created in the centroid of the element. Some strategies allow the creation of the new node at the edges or faces of the elements [26] in order to improve the quality of the elements respecting the contours of the domain.
• All the elements of the mesh which circumscribed sphere (or circumcircle in 2D) includes the new node are deleted. Note that a void region is created containing the new node inside.
• A Delaunay triangulation is created with the new node and the contour nodes of this void region.

This procedure leads to a Delaunay mesh (accomplishing the Delaunay condition), but this does not guarantee a given level of quality by itself. Often, some elements can present very low quality (especially in 3D cases). These elements may have null volume and are called slivers. Even if its volume is equal to zero, they can accomplish the Delaunay condition. For this reason, it is common to relax the Delaunay condition in some regions in order to avoid quality problems [27].

The main advantages of this method are:

• Robust. Mathematically, we can always obtain a Delaunay mesh from a cloud of points. However, sometimes the modifications of the method to cover the requirements of the simulation make it not so robust.
• Fast. Although it depends strongly on the implementation, it can be said that Delaunay methods are naturally faster than the advancing front ones.

The main drawbacks of the Delaunay method are:

• Naturally not constrained. Delaunay methods use the representation of the contour of the domain as an initial configuration for the insertion of nodes procedure, but the final mesh boundaries are not guaranteed to be constrained with that contour. Some strategies can be applied [28,29,30] in order to guarantee that the nodes in the initial contours will lay on the contour of the final mesh, or even a constrained condition in the faces of the initial contour, but these modifications reduce the performance and robustness of the method. Furthermore, if there are huge different mesh desired sizes in the domain, these strategies are not very robust.
• Requires watertight input geometry. As in the advancing front method, the contours of the domain are required to be watertight. Actually, a strategy could be followed to treat non-watertight geometries. It is based on generating the first mesh (the convex hull in the traditional Delaunay methods) using only the nodes of the contours, with an automatic recognition of the boundaries: the so called alpha-shape method [31]. This would indeed generate a volume mesh from a non-watertight geometry, but there is no enough control in the boundaries recognition to ensure that mesh correctly represents the topology of the domain.
• Not naturally parallelizable. The check of the Delaunay condition of an element requires taking into account all the mesh entities participating in the circumscribed sphere of the element. As its radius depends on the position of the nodes of the element (cannot be bounded a priori), the treatment of one element can involve the whole mesh in some configurations. This implies a hard dependency ranging from an element to the whole mesh, so a parallel implementation of the method is not obvious. However, some parallel implementations (mainly for shared memory paradigm) of the Delaunay method have been carried out, [32] or [33].

## 2.5 Space decomposition methods

Space decomposition-based methods follow a different philosophy than the methods explained before. To generate the mesh, they basically subdivide the space into cells providing with a spacial decomposition covering the space where the domain is (overlapping the domain). These cells can be thought in a general way, but it is common to use one of the following main structures which govern their configuration:

• Bin: a bin structure provides with an homogeneous grid as the space decomposition formed by regular cells (squares in 2D or cubes in 3D). A graphical view of a 2D bin is shown in Figure 13(a).
• Octree: an octree (quadtree for the 2D case) is basically a hierarchical spatial structure that partitions the 3D space into regular cells [34]. These cells can be refined in given zones of the domain. A graphical view of a quadtree is shown in Figure 13(b). A more detailed definition of the octree structure is given in Section 3.2.
 (a) (b) Figure 13: 2D examples of typical structures used in space decomposition based meshing algorithms. (a) Bin. (b) Quadtree.

The bin structure is suitable for homogeneous discretizations. It is common to use an octree structure as the regular grid (which gives the name of the family of methods), because it is more flexible for mesh generation purposes, as domains to be meshed and desired mesh sizes are commonly non-homogeneous. Octree-based meshers were pioneered by Yerry and Shephard in  [35] and, since then, several approaches have been proposed [36,37,38,39,40,41,42].

The octree structure was thought for the first time for space searching purposes [34], and the specific topology of the spatial decomposition it represents (detailed in Section 3.2) gives several advantages for mesh generation. Somehow, an octree can be considered a mesh itself (it can be thought as a non-conformal hexahedra mesh), so it is really natural to build a mesh from it.

Although several algorithms have been proposed parting from the octree-based family of methods, almost all of them follow three main steps:

• Generate the regular grid for the space decomposition.
• Generate the elements of the final mesh directly using given patterns from the cells of the regular grid.
• Fit somehow the boundaries of the domain.

Concerning the first step, as it has been pointed out before, the octree is the most common structure used for the space decomposition. From the theoretical point of view, all octrees are similar, but depending on the way the octree will be used, different implementations have been proposed by several authors [43,34] in order to improve the efficiency of the octree, the performance for searching processes, the optimization considering the memory needed to store it, etc... More details on the implementation of the octree are presented in Section 6.2.

The generation of the elements of the final mesh from the octree is a simple process. It is based in creating the mesh elements directly from the octree cells (the definition of octree cell, as well as other octree related basic concepts are explained in detail in Section 3.2). Some of the existing methods apply different splitting patterns from the cells to get tetrahedra  [35,44,37]. Other methods can get directly the cells of the octree as hexahedra elements of the final mesh (in cases where the final mesh is not needed to be conformal), or create transition elements when two neighbors present hanging nodes [42].  [40] proposes a different approach: create special cells where the octree is refined (where the neighbor elements are not conformal) and build the dual of the octree. The cells of it are directly the elements of the final mesh. With this approach a final mesh of conformal hexahedra is obtained automatically.

The key difference of each method remains in the third step, which is the most complex one. As the octree is a regular space decomposition, their cells do not fit exactly the contours of the domain to be meshed. Getting only the elements coming from the cells which are in the inner part of the domain (or even the ones intersecting its boundaries), the contours of the final mesh are staircase-like, so they are not able to represent smooth shapes with a given curvature.

Considering the inner cells (the ones totally inside a volume) and the interface ones (the ones colliding with the contours of a volume), different strategies have been proposed to fit the contours:

• Get only the inner cells and fill somehow the empty space between them and the real contour of the domain. This family of methods profits from the main advantages of octree-based methods for the inner part of the domains, but presents the same limitations near the contours as the methods used to fill those empty spaces.
• Project onto the contour of the domain the boundary nodes of the elements created from the inner cells using different techniques [36,39]. This strategy warps the nodes coming from the octree, so the shape of the octree changes, but its topology remains the same. [42] proposes a pillowing technique near the contours to avoid bad shaped elements (hexahedra) after mapping the nodes.
• Move some nodes and split the elements intersecting the contours of the domain in order to represent it precisely [37].

Although these methods achieve the smooth representation of the contours, they have to follow specific strategies to preserve the geometrical features (corners or sharp edges). [41] proposes a re-tetrahedralization of the octants of the octree containing sharp edges using advancing front or Delaunay technique, taking into account the intersection points between the sharp edges and the octree cells.  [40] follows a strategy based on detecting which triangle from the input boundary the final nodes lay onto, and assuming that if two nodes lie onto two triangles connected by a sharp edge, there should be a sharp edge between those nodes of the final mesh. This strategy is not so robust, as it assumes that the sizes of the final elements is quite similar to the sizes of the triangles of the contour, so the triangles where two neighbor nodes of the final mesh lie onto are supposed to be neighbors connected by an edge. This is not a general situation.

As it has been explained, several approaches have been proposed departing from the octree-based family of methods. Although each approach has its own characteristics, some common advantages can be detected:

• Robust. The operations involved in most common octree-based methods are designed to be robust independently on the tuning of their control parameters. This is applicable to the operations dealing with the octree itself, but some processes dealing with the mesh generated (especially for body-fitted meshes) may not be so robust.
• Fast. As mentioned above, the octree is a structure used to make faster the searching algorithms in space. Indeed, almost all the meshing algorithms use an octree as a tool for searching purposes. In this case, this structure is also the base of the mesh generation itself, so the algorithm takes profit from it to improve its performance. Furthermore, the inner elements of the mesh are created directly based on a division pattern of the octree cells, which is really fast in comparison with other methods which have to take care on geometrical aspects when generating each element.
• Naturally parallelizable. The octree itself is a partition of the space, so it is really easy to identify the parts of the domain affected when some part of the mesh is modified or generated. Furthermore, typical operations involved in an octree-based mesher are performed at octree cell level, and most of them are independent from one cell to another, so the parallelization of some parts of the algorithm is almost automatic. Several parallel implementations of the octree structure have been carried out taking advantage on these features [TU:2005:SPO:1105760.1105769,CamataCoutinho].
• Good quality meshes. As the mesh elements in these methods come from a predefined pattern of decomposition of an octree cell, all the inner elements are in a range of qualities known a priori. However, the elements in the contour of the domain are typically more distorted (depending on the method used).
• Allow not cleaned geometry as input. Most octree-based methods only take care on the position in space of the boundaries of the domain rather than a topological relationship between the contour entities. This is a key point to allow some of the octree-based meshers to work with non-watertight input geometries. This property also makes octree-based meshers independent from the quality of the contour meshes used as input.

The main drawbacks of octree-based methods are:

• Naturally not constrained. The space decomposition of the octree does not respect the topology of the model. Octree-based methods are often not constrained, and the topology preservation is usually linked to the sizes of the octree cells used. This ends with a need for the user of the mesher to assign properly the desired sizes of the mesh, taking into account some kind of characteristic size of the model in specific regions of the domain.
• Hard to preserve geometric features. This drawback is strongly linked with the previous one: as the mesher is not constrained, it is hard to preserve sharp edges or corner points from the input data in the final mesh. Furthermore, cases where the sharp edges involve a small dihedral angle are very unfavorable for these methods, as the octree cells are very regular by its nature: it is hard to represent a small angle parting from a regular shape.
• Predefined size transitions. As most elements come from a direct pattern from the octree cells, the size of an element in comparison with the size of its neighbor is strongly related with the sizes ratio between neighbor cells in the octree. This link makes quite stiff the possibility of applying a smooth size transition in the final mesh. For some simulation methods, the meshes obtained by octree-based meshers present (locally) too sharp variations in the sizes of neighbor elements.
• Alignment of edges in a preferred direction. Because of the regular partition of space the elements come from, these meshes tend to generate meshes aligned with it. This alignment may influence the solution of the numerical simulation depending on the method used.
• Creation of too refined meshes. Octree-based meshers tend to generate a large number of elements, specially when the model has sharp features. This is because these methods often try to guarantee a minimum element quality by refining the octree near them.

## 2.6 Proposed solution

The strategy chosen in this work to cover all the requirements described in Section 1.2 is to develop an octree-based mesher. The election of an octree-based mesher in this work has been made taking into account the main advantages and disadvantages of the different methods:

• Structured and semi-structured meshers have been disregarded because the requirement for the input geometry to present a given topology is so restrictive that it makes them non applicable for many complex cases.
• Advancing front like methods have not been taken into account because the hard dependency they have on the quality of the input mesh for being successful in the mesh generation. Furthermore, they are directly not applicable if the input geometry is not watertight, unless some geometry repairing process is applied prior to the meshing itself.
• Delaunay methods have not been taken into account because they are not naturally constrained, and the proposals for them to be constrained blur some of its advantages. Furthermore, Delaunay methods are not naturally parallelizable, and they often generate localized bad quality elements.
• Octree-based methods are fast, robust and lead to very good quality meshes. Their main disadvantage are the inability to preserve the input topology, and their not constrained nature. However, the advantages of this family of methods, together with the ideas presented in this monography to overcome these drawbacks have tipped the balance to these methods.

As explained in Section  1.2.4, although the main objective of this work is to develop a new volume mesher, the methodology proposed is applicable to generate meshes of 3D surfaces and lines not belonging to any volume. The case of lines is automatically solved with the special treatment of line elements in the volume mesher (Section 5.3.1), and some adaptations are maid to the volume mesher in order to mesh surfaces as it is explained on Section 5.3.10.

# 3 Basic concepts of the new mesher

This chapter focuses on defining the different concepts involved in the proposed meshing algorithm, as well as some auxiliary algorithms needed to understand it. The following concepts will be described:

• Input data (Section 3.1). The input data needed for the mesher is defined, as well as the way of interaction of the mesher with a CAD system.
• Octree structure (Section 3.2). As the octree is the base structure for the mesher, an introduction to it is carried out in this section, highlighting its main characteristics.
• Octree properties for mesh generation (Section 3.3). Some specific properties of the octree are used in the meshing algorithm. A detailed explanation of them, as well as the key notation is introduced.
• Geometrical intersections (Section 3.4). The way the new algorithm deals with the geometrical intersections is presented. Also, the notation of pathological intersection types is introduced. This issues will be referred to in other parts of the work.

## 3.1 Definition of input data

The essential input for the mesher is the geometrical definition of the boundaries of each volume of the domain. As indicated in Section 1.2.2, this definition can be carried out using CAD or mesh entities. In this document the general concepts of surface, line and point entities will be used for both representations.

At this point, it has to be commented that the mesher considers the outer part of the domain as another volume to be meshed. It takes the name of outer volume, or volume number zero. Of course, the outer volume extends until infinite and it has no sense to consider it as a closed volume, so it is treated in a special manner. Later on it will be explained in more detail how the mesher deals with it.

Together with the surface entities an extra information is needed: the identification of the volumes each surface entity is interfacing. Note that considering the outer part of the domain as a virtual volume, all the surface entities defining the domain are interfacing two volumes. It may have sense for a surface entity to interface more than two volumes, but this would imply overlapping definitions of the 3D space (parts of space belonging to more than one volume). These kinds of topology are not considered in the present work.

Apart from the geometrical definition of the boundaries of each volume of the domain, extra information can be given to the mesher in order to specify some characteristics of the final mesh. It is important to note that this extra information is optional, as the mesher should generate the mesh of the domain with or without it. This information is given by the mesh size entities, the forced point entities, the forced line entities and the general parameters. Hereafter the characteristics of these data are detailed:

• Mesh size entities. The mesh size entities are geometrical entities used to provide the mesher with the desired size of the mesh in different regions of the 3D space. These geometrical entities can be point, line, surface or volume entities (both in mesh or CAD form), and they have a desired mesh size associated to them. These entities can be part of the input boundaries or not. If a desired mesh size is needed to be assigned to a volume, it can be provided just with the identification of the volume, without the need of creating the corresponding geometrical entities (volume entities in this case).
• Forced point entities. These are positions in space where the final mesh is forced to have nodes. They can be part of the input boundaries or not.
• Forced line entities. These are line entities to be preserved in the final mesh. This means that the final tetrahedra mesh will have a connected path of edges identifiable as a linear mesh of each forced line entity. As an example, the line entities in Figure 14(a) can be considered as forced line entities, as a path of connected edges representing it can be identified in the mesh shown in Figure 14(b). Forced line entities can be part of the input boundaries or not.
 (a) (b) Figure 14: (a) Contours of a volume highlighting some of its forced line entities. (b) View of the tetrahedra mesh of the volume highlighting the sharp edges corresponding to the forced line entities in (a).
• General parameters. These are some global parameters which govern the behavior of the mesher and are not attached to geometrical entities. These parameters are optional, as the mesher should work without them, but they can help to fit the final mesh into the requirements of the simulation. Hereafter there is a list of these parameters:
• General mesh size. This is the mesh size desired in the regions of the domain where there are no mesh size entities. If this parameter is not entered, the general mesh size is taken as the size of the bounding box of the domain.
• Size transition factor. It controls the behavior of the size transition function. The definition and use of this function is explained in Section 5.2.2, and it basically determines if the size of the elements grows faster or slower from the areas where the mesh size is smaller to the ones where it is larger.
• Maximum angle for sharp edges. This is the dihedral angle below which an edge of the input geometry is considered to be preserved in the final mesh.
• Maximum chordal error. Two values may be needed: a relative and an absolute limit for the chordal error allowed for the mesh to be generated.
• Maximum distance between overlapping boundaries. This parameter is named ${\textstyle tol_{nw}}$. It is only applicable in case of non watertight geometries. It represents an estimation of the maximum distance between overlapping surfaces interfacing the same volumes. The definition and use of this parameter is explained in Section 3.4.

Apart from all these data, extra information can be attached to the entities defining the input boundaries: the external data. The external data can be of any nature and it has no relevance for the meshing process, but the mesher will transfer it to the corresponding entities of the final mesh. As an example, if a surface entity defining part of the boundary of a volume has some external information, the nodes of the final mesh placed on that surface entity, or the faces of the tetrahedra with all its nodes onto it will have the same external data attached to them.

### 3.1.1 Integration with CAD data

To prepare all the data needed for a numerical simulation it is common to use a software tool: the pre-processor. As explained in Chapter 1, part of these data is the mesh representing the geometry of the domain. As the geometry of the domain is often provided in a CAD format, pre-processors typically are forced to work with CAD data and generate the meshes for the simulation. This is the reason why pre-processors are often CAD systems, and they include several meshers inside.

From the point of view of a general pre-processor or a CAD system willing to use the presented mesher, there are some interesting aspects to be considered. Typically the geometrical entities inside these systems have several information attached related with the simulation data (boundary conditions, material properties, etc...) or to the CAD system structures (layers grouping the entities, topological information, etc...). This information must be transferred to the generated mesh, which implies the following requirement for the mesher: it should return specific meshes representing given geometrical entities (not only volumes, but also curves and surfaces).

The mesher should provide not only with the tetrahedra generated, but also with triangular meshes (made of triangles which are faces of the tetrahedra) and linear meshes (made of linear elements which are edges of the tetrahedra). Actually point entities can also be identified with a final node in the mesh, just by setting the corresponding point entity as a forced point entity.

The way to make the mesher returns (apart from the tetrahedra meshes of the volumes of the domain) the mesh of some line or surface entities from the input data is explained hereafter.

How to get the mesh of a line entity:

• Set the line entity as a forced line entity in the input data.
• Assign a specific external data to the line entity.
• The mesher will return a set of linear elements with the same external data attached.

How to get the mesh of a surface entity:

• Considering the line entities which are contour of the surface entity, set them as forced line entities in the input data. Note that the contour line entities of a surface entity form a closed path of line entities (actually, more than one set of closed paths can be present when the surface entity has holes).
• Assign a specific external data to the surface entity.
• The mesher will return a set of triangular elements enclosed by the forced line entities with the same external data attached. This situation can be seen in Figure 6: the patch of triangles in Figure 6(c) is the surface mesh of the surface entity represented by the highlighted triangles in Figure 6(a) (region ${\textstyle A}$). To achieve this situation, the line entities surrounding the region ${\textstyle A}$ should be set as forced line entities.

Following this mechanism, the pre-processor can obtain not only the mesh of the volumes of the domain, but also the mesh of any line or surface entity. Once a mesh of an entity can be identified, all the information of the input entity can be transferred to the corresponding part of the final mesh.

## 3.2 Octree structure

As the octree is a key structure in the proposed mesher, a brief introduction is given in this section to highlight its characteristics.

 (a) (b) (c)
Figure 15: Example of a quadtree structure. (a) The root cell of a quadtree. (b) Root cell subdivided in 4 cells. (c) Example of a quadtree refined 4 levels.

An octree (quadtree in the 2D case) is basically a hierarchical spatial structure that partitions the space [34]. The basic structure of an octree is the cell, which is a cubic portion of space (square in the 2D case). Actually, the cell can be a parallelepiped (or parallelogram in 2D). In this work the octree used is an homogeneous one, which implies that the cells are regular parallelepipeds (cubes). From a first cell which is the bounding box of the space to be partitioned (the so called root of the octree), a successive subdivision can be performed, where each cell is subdivided in eight cells (four in 2D case). These eight cells are the sons of the cell they come from, which is their father. Cells with no sons are called leaves.

In the Figure 15 a graphical view of the root of a quadtree and different levels of refinement are shown.

Considering the root of the octree, it can be equilateral or not. This, together with the way the cells are subdivided, leads to different configurations of the octree (Figure 16). In particular, for notation purposes, an octree accomplishing the following two properties receives the name of isotropic octree:

• The octree root is equilateral.
• The division of a cell results in eight sons, and the division criterion is equidistant (the eight sons of a cell are identical). In other words, the division of a cell is done parallel to its faces passing by its center.

To fix the notation, some interesting concepts related with the octree are detailed hereafter:

• In a cell there can be identified vertices, edges and faces (analogously as the nodes, edges and faces of an hexahedral mesh element).
• The cell size is the length of the edge of the cell.
• Two cells are neighbors by face if they are in contact by a face.
• Two cells are neighbors by edge if they are in contact only by an edge.
• Two cells are neighbors by vertex if they are in contact only by a vertex.
• An octree is refined if some of its leaves is subdivided.
• The degree of neighborhood between two cells is defined as the minimum number of cells needed to go from a cell to the other traveling by neighbor cells. As an example, the degree of neighborhood between two neighbor cells is one.
• A maximum depth of cells subdivision can be defined. Lets call it ${\textstyle N_{L}}$. A level can be assigned to each cell, so as: the root has level equal to ${\textstyle N_{L}-1}$, and the sons of a cell has one level less than their father. Note that the smallest cell of the octree has a level equal to 0.

## 3.3 Specific octree properties for mesh generation

As explained in previous sections, the octree structure was thought for the first time for space searching purposes [34]. In this section, the adaptations to the octree structure done in order to use it for mesh generation and its main properties to understand the algorithm are explained.

The decision of using the octree for isotropic mesh generation leads to use an isotropic octree. This decision has an important relevance at the time of implementing the algorithm (as it will be seen in Section 6.2).

The proposed meshing algorithm should be valid using other kinds of octree (like the analogous quadtree examples shown in Figures 16(a) and (b)), but it would lead to non isotropic meshes.

 (a) (b) (c) (d)
Figure 16: Different kinds of quadtree. (a) Quadtree with a non-equilateral root, with a equidistant cell division criterion. (b) Quadtree with an equilateral root, with a non-equidistant cell division criterion. (c) Isotropic non-balanced quadtree. (d) Isotropic balanced quadtree.

Another important characteristic of the octree chosen for the method is the so called constrained two to one condition. This is a widely used condition in octree based meshers, and limits the number of neighbors of a cell. The two to one name comes from the two dimensional case (quadtree), and limits the maximum number of neighbors of a cell to two. In the octree case (3D), this condition implies that a cell cannot have more than four neighbor cells by face, or two by edge. In the present document an octree accomplishing the constrained two to one condition is referred as a balanced octree. Two configurations of an isotropic quadtree (non-balanced and balanced) are shown in Figures 16(c) and (d).

The main reason to use a balanced octree for the proposed meshing algorithm is to simplify the patterns to build the tetrahedra from the octree cells, ensure a better quality in the final tetrahedra and avoid a very strong sizes transitions in the final mesh. The tetrahedra generation process from the octree has the following characteristic: the more difference between sizes of neighbor cells, the worse aspect ratio will have the tetrahedra generated from them.

### 3.3.1 Octree cell types

As it has been pointed out in previous sections, octree cells are the result of the space partitioning by the octree structure. Besides this, they have some information attached as the input boundary entities colliding with them.

 Figure 17: 2D example where the three kinds of cells can be identified. The black curved line defines the contours of a domain formed by two surfaces (which are in contact), and the light gray lines represent the octree. Outer cells are the white ones, interface cells are marked with dots, and inner cells are colored in gray.

In the frame of this work, the octree cells are classified in three categories:

• Interface cells. These are the cells colliding some input boundary entity. This means that some input boundary entity is inside them, or intersects them. Also the cells containing a forced node (Section 5.3.2) are considered as interface cells.
• Outer cells. These are the cells which are completely out of the domain to be meshed and they do not contain any forced node.
• Inner cells. These are the cells which are not interface or outer ones. These cells can be identified with some volume of the domain, as they are totally inside one volume.

Figure 17 shows a 2D example where the three kinds of cells can be identified.

### 3.3.2 Octree positions and nodes

As it is explained in Section 3.3.3, tetrahedral elements will be generated following a pattern from the octree. The nodes of these tetrahedra are called octree nodes and they are assigned to some predefined positions in space: the octree positions. Each cell has 27 octree positions corresponding to the vertices of the cell (8), the center of the cell (1), the center of its edges (12) and the center of its faces (6). A graphical view of the octree positions of a cell is shown in Figure 18.

 Figure 18: Octree positions of an octree cell. The cell is represented by the black lines.

It has to be noted that an octree position can be shared by more than one cell, as it can be seen in Figure 19.

 Figure 19: Linear octree positions of a part of an octree. White dots are center of cells, and black dots correspond to vertices of cells.

The linear positions of an octree cell are defined as the ones corresponding to the vertices and the center of the cell. The other positions are called quadratic positions. As an extension, the term of linear or quadratic octree node can be used to refer an octree node associated to a linear or quadratic cell position.

When referring to the whole octree, a linear position (or octree node) is the one which is linear in some cell. It has to be noted that an octree node can be linear regarding one cell, but quadratic from the point of view of another cell containing it. The linear octree positions of a part of an octree are shown in Figure 19 as an example.

It is important to note that not all the octree positions are forced to have an octree node associated, but all the octree nodes are linked to an octree position.

In Section 5.3.2 the concept of forced node is introduced. It basically corresponds to a node linked to an octree position, but occupying a different position in space.

### 3.3.3 Tetrahedra patterns

Parting from a given octree configuration, there are several possible tetrahedra patterns to be applied in order to split it. One important consideration is that the pattern chosen must fill the space with tetrahedra in a conformal way: it must not leave hanging nodes. A hanging node is a node lying on an element not being a vertex of it (it is on an edge or a face).

The option chosen in this work is based on the body centered cubic (BCC) lattice, which lead to a space-filling tetrahedra [45]. The use of BCC patterns linked to octree structures is quite natural considering the spacial distribution of the octree, and was proposed by [DBLP:conf/imr/Fuchs98,NME:NME616]. This option fills the space in a conformal way with a set of identical high quality tetrahedra: dihedral angles are ${\textstyle 60}$ or ${\textstyle 90}$ degrees, and edge lengths are 1 and ${\textstyle {\dfrac {\sqrt {3}}{2}}}$ times the cell size. This provides the tetrahedra with an edge ratio of ${\textstyle 1.155}$ (the ratio of the longest and shortest edges of the element).

A BCC lattice based pattern is local in the sense that each tetrahedra generated only depends on one cell and one of its neighbor. A pattern only depending on each cell (independent from the neighbor ones) would be more efficient for parallelization. However, this kind of patterns often provides with a lower quality tetrahedra (although their quality is acceptable) and a larger number of them [46].

The BCC lattice is defined in a regular grid (an octree with all its leaves equal-sided). As the octree used in the proposed method is not regular (it has leaf cells in different levels), other tetrahedra patterns must be defined.

To reach a tetrahedra mesh with no hanging nodes, all the linear octree nodes are used, and (in case they exist) the forced nodes linked to quadratic positions (forced nodes are defined in Section 5.3.2).

Basically, the tetrahedra pattern proposed focuses the tetrahedra generation cell by cell, and focusing on one cell, using its six faces independently. The only information required to generate the tetrahedra from a face ${\textstyle F}$ of a cell ${\textstyle C}$ is: the octree nodes of the face ${\textstyle F}$, the center octree node of the cell ${\textstyle C}$ and the neighbor cell of ${\textstyle C}$ from the face ${\textstyle F}$, denoted by ${\textstyle C_{neig}}$.

As explained previously, the octree used is balanced (it accomplishes the constrained two to one condition). This ensures that the level of ${\textstyle C_{neig}}$ is one less, equal, or one more than the level of ${\textstyle C}$. These three configurations lead to the different tetrahedra patterns:

 (a) (b) (c) (d) Figure 20: Tetrahedra pattern in case where cell ${\displaystyle C}$ has the same level as ${\displaystyle C_{neig}}$. Tetrahedra generated from edge ${\displaystyle e}$ (${\displaystyle {\overline {e_{1}e_{2}}}}$) of common face between ${\displaystyle C}$ and ${\displaystyle C_{neig}}$ in different situations. (a) No quadratic octree nodes involved: one tetrahedron is generated (${\displaystyle c_{1}e_{2}c_{2}e_{1}}$). (b) Octree node in the center of face ${\displaystyle F}$: two tetrahedra are generated(${\displaystyle c_{1}e_{2}f_{c}e_{1}}$ and ${\displaystyle f_{c}e_{2}c_{2}e_{1}}$). (c) Octree node in the center of edge ${\displaystyle e}$: two tetrahedra are generated(${\displaystyle c_{1}e_{2}c_{2}e_{c}}$ and ${\displaystyle c_{1}e_{c}c_{2}e_{1}}$). (d) Octree nodes in the center of face ${\displaystyle F}$ and edge ${\displaystyle e}$: four tetrahedra are generated(${\displaystyle c_{1}e_{2}f_{c}e_{c}}$, ${\displaystyle c_{1}e_{c}f_{c}e_{1}}$, ${\displaystyle f_{c}e_{2}c_{2}e_{c}}$ and ${\displaystyle f_{c}e_{c}c_{2}e_{1}}$).
• ${\textstyle C}$ and ${\textstyle C_{neig}}$ have the same level (this means they have the same size). In this configuration the center node of ${\textstyle C_{neig}}$ is used to build the tetrahedra. Taking into account that the face ${\textstyle F}$ (the face in contact between ${\textstyle C}$ and ${\textstyle C_{neig}}$) has four sides (edges), from each edge ${\textstyle e}$ from the face ${\textstyle F}$ one, two or four tetrahedra are generated, depending on the possible presence of octree nodes in quadratic positions. The tetrahedra generated under each situation are shown in Figure 20.

The minimum dihedral angle of the tetrahedra generated in this configuration is ${\textstyle 45}$ degrees, and it occurs when some node in a quadratic position is used.

This case involves the creation of tetrahedra in two cells. Hence, only one of them should create them to avoid repetitions. The criteria used is that cell ${\textstyle C}$ creates the tetrahedra only if it is lower than ${\textstyle C_{neig}}$, otherwise the tetrahedra will be created by ${\textstyle C_{neig}}$. Considering two leaves (${\textstyle A}$ and ${\textstyle B}$), ${\textstyle A}$ is lower than ${\textstyle B}$ if the ${\textstyle x}$ coordinate of the center of ${\textstyle A}$ is lower than the ${\textstyle x}$ coordinate of the center of ${\textstyle B}$. If the ${\textstyle x}$ coordinates are equal, the ${\textstyle y}$ coordinate is checked, and if it is also equal, the ${\textstyle z}$ coordinate is used to compare the cells. It has to be noted that, because of the characteristics of the octree structure, there cannot be two leaves with their centers in the same position.

• ${\textstyle C}$ has one level more than ${\textstyle C_{neig}}$ (this means ${\textstyle C}$ is bigger than ${\textstyle C_{neig}}$). In this configuration eight tetrahedra are created from the face ${\textstyle F}$. These tetrahedra are shown in Figure 21(a). It has to be considered that if ${\textstyle C}$ is bigger than ${\textstyle C_{neig}}$, the constrained two to one condition ensures that there will be four cells neighbors of ${\textstyle C}$ in the other side of ${\textstyle F}$.

The minimum dihedral angle of the tetrahedra in this configuration is ${\textstyle 45}$ degrees.

 (a) (b) Figure 21: Tetrahedra pattern in case where cell ${\displaystyle C}$ has different level than ${\displaystyle C_{neig}}$. (a) ${\displaystyle C}$ is bigger than ${\displaystyle C_{neig}}$: 8 tetrahedra are generated (${\displaystyle f_{1}f_{5}f_{c}c_{1}}$, ${\displaystyle f_{5}f_{2}f_{c}c_{1}}$, ${\displaystyle f_{2}f_{6}f_{c}c_{1}}$, ${\displaystyle f_{6}f_{3}f_{c}c_{1}}$, ${\displaystyle f_{3}f_{7}f_{c}c_{1}}$, ${\displaystyle f_{7}f_{4}f_{c}c_{1}}$, ${\displaystyle f_{4}f_{8}f_{c}c_{1}}$ and ${\displaystyle f_{8}f_{1}f_{c}c_{1}}$). (b) ${\displaystyle C}$ is smaller than ${\displaystyle C_{neig}}$: two tetrahedra are generated(${\displaystyle f_{1}f_{2}f_{3}c_{1}}$ and ${\displaystyle f_{3}f_{4}f_{1}c_{1}}$).
• ${\textstyle C}$ has one level less than ${\textstyle C_{neig}}$ (this means ${\textstyle C}$ is smaller than ${\textstyle C_{neig}}$). In this case two tetrahedra are created from the face ${\textstyle F}$, shown in Figure 21(b). Note that this configuration could build the two tetrahedra crossing the face ${\textstyle F}$ by diagonals ${\textstyle f_{1}f_{3}}$ or ${\textstyle f_{2}f_{4}}$. To ensure the final mesh to be conformal, the diagonal chosen is the one passing by the octree node which is the center of face in ${\textstyle C_{neig}}$.

The minimum dihedral angle of the tetrahedra in this configuration is ${\textstyle 45}$ degrees.

It can be seen that the quality of the tetrahedra generated using these patterns depends on the configuration chosen, but it is always very good: the minimum dihedral angle of any tetrahedra coming from this predefined patterns is ${\textstyle 45}$ degrees.

The case where a cell has no neighbor by one of its faces occurs in the cells in contact with the boundaries of the octree root. This case is not considered because the building of the octree root (Section 5.2.3) ensures these kind of cells are totally out of the domain to be meshed, so they are not used to generate any tetrahedron.

## 3.4 Geometrical intersections

Geometrical intersections play a key role in the new meshing algorithm and affect several parts of it. This section focuses in defining how this intersections are considered. Their treatment has special interest when the input boundaries are non-watertight.

The intersection operation involved in the algorithm is typically the one between a segment and the surface entities defining the contours of a volume. This segment is understood as a portion of a straight line between two points. Depending on the part of the algorithm, it could be referred with different names such as edge (when talking about tetrahedra or triangle edges) or ray (when talking about ray casting operations). In this document some figures are based in 2D examples for clarity purposes. The equivalent intersection operation in these cases is between segment and line entities defining the input boundaries.

 (a) (b) (c) (d) (e) (f)
Figure 22: Types of intersections between a segment and a surface entity. Crosses are the intersection points. (a) No intersection. (b) ${\displaystyle S}$ type intersection. (c) ${\displaystyle C}$ type intersection. (d) ${\displaystyle T}$ type intersection. (e) ${\displaystyle PT}$ and (f) ${\displaystyle PI}$ type intersection; the red thick part of the segment is co-planar with the surface entity.

Five situations are distinguished when evaluating the intersection between a segment and a surface entity (a graphical interpretation of these cases using simple examples is depicted in Figure 22):

• No intersection: the segment and the surface entity do not intersect (Figure 22(a)).
• Single intersection (${\textstyle S}$ type): they intersect in one point within the surface entity (Figure 22(b)).
• Contour intersection (${\textstyle C}$ type): they intersect in one point laying on the contours of the surface entity (Figure 22(c)).
• Tangent intersection (${\textstyle T}$ type): they intersect in one point tangentially to the surface entity (Figure 22(d)) and the segment remains in the same side of the surface entity at both sides of the intersection point.
• Co-planar intersection (${\textstyle P}$ type): the segment and some part of the surface entity are parallel and the distance between them is zero. This case can be split in two:
• Tangent (${\textstyle PT}$ type): before and after the intersection, the segment remains at the same side of the surface entity (Figure 22(e)).
• Intersected (${\textstyle PI}$ type): if the ${\textstyle PT}$ case is not accomplished (Figure 22(f)).

Analytically, the ${\textstyle P}$ type of intersection has an infinite number of intersection points. As it will be seen in further sections, these intersection types only take relevance in the node coloring process (Section 4), that determines inside which volume a node is. For this reason the treatment of ${\textstyle P}$ type intersections is detailed in Section 4.2.

Types ${\textstyle S}$, ${\textstyle C}$ and ${\textstyle T}$ have only one intersection point. However, as the contours of a volume may be formed by more than one surface entity, each intersection may involve a different number of intersection points (one for each intersected surface entity). This situation is illustrated in the 2D examples of Figure 23 (as a 2D example, line entities play the role of surface entities in 3D). In Figure 23(b) a ${\textstyle C}$ intersection type has two intersection points (one for each intersected line entity), as well as in Figure 23(c) with the ${\textstyle T}$ type intersection.

 (a) (b) (c)
Figure 23: Line entities enclosing surface ${\displaystyle A}$ intersected by a segment (bounded by two dots) presenting different intersection types. Intersected line entities are drawn in dotted line. (a) ${\displaystyle S}$ intersection type. (b) ${\displaystyle C}$ intersection type. (c) ${\displaystyle T}$ intersection type.

These cases are characterized by the fact that all the involved intersection points are really close one from each other. Theoretically, all the intersection points should be in the exact same position, but because of the tolerances involved in the numerical computation of intersections this cannot be guaranteed. The tolerance ${\textstyle tol_{c}}$ is defined in order to determine if a collection of close intersection points corresponds to the same intersection: if all of them are within a distance lower than ${\textstyle tol_{c}}$, they are collapsed into a multiple intersection point (MIP). The value of ${\textstyle tol_{c}}$ is a portion of the model bounding box size:

 ${\displaystyle tol_{c}=\alpha _{bb}\cdot s_{box}}$
(3.1)

being ${\textstyle s_{box}}$ the length of the minimum side of the model bounding box, and ${\textstyle \alpha _{bb}}$ a real value between zero and one. The value used in the present work for ${\textstyle \alpha _{bb}}$ is detailed in Section 6.7.

The position of the MIP corresponding to a collection of intersection points is the mean position between all of them. In the present work, when evaluating the intersection point between a segment and the contours of a volume, in cases where there are intersection points close enough, their MIP is considered instead of them. This matches with the theoretical number of intersection points corresponding to each intersection type: ${\textstyle C}$ and ${\textstyle T}$ intersection types have only one intersection point (a MIP).

### 3.4.1 Non-watertight geometries

For non-watertight geometries, special cases may be treated when evaluating the intersection between a segment and the contours of a volume. There are basically two specific intersection types (a 2D example of this cases is depicted in Figure 24):

• Overlap intersection (W type). This situation happens when the part of the boundaries of the domain where the segment intersects is defined by more than one overlapped surface entity (Figure 24(a)).
• Gap intersection (G type). This situation happens when the part of the boundaries of the domain where the segment intersects has a gap (Figure 24(b)).
 (a) (b) Figure 24: Intersection between a segment and the line entities enclosing surface A in non-watertight situations: (a) overlap (W intersection type) and (b) gap (G intersection type).

For the W intersections type, the corresponding MIP of the intersection points involved is considered. For this purpose, in cases where non-watertight geometries define the domain, ${\textstyle tol_{c}}$ takes the value of ${\textstyle tol_{nw}}$. This parameter must be provided in the input data and is a characteristic length of the gaps and the overlappings of the model. ${\textstyle tol_{nw}}$ must be an upper limit of the distance between overlapping entities. The MIP corresponding to the 2D example of Figure 24(a) is shown in Figure 25(a).

G intersections do not provide with any intersection point. However, as it will be seen in further sections, there are cases where the extremes of a segment are known to be inside different volumes. Situations where both extremes belong to different volumes indicate that there should be an intersection point although it is not detected. In this cases the gap intersection point (GIP) is created. Considering the surface entities surrounding the segment, the closest point from each surface entity to the segment can be computed. The GIP takes the position of the closest one. The GIP corresponding to the 2D example of Figure 24(b) is shown in Figure 25(b).

 (a) (b) Figure 25: Treatment of intersections of non-watertight geometries depicted in Figure 24. (a) MIP corresponding to the W intersection depicted in Figure 24(a). (b) GIP corresponding to the G intersection shown in Figure 24(b).

In the present work, the intersection point between a segment and the contours of a volume could be a single intersection point, a MIP or a GIP.

# 4 Coloring algorithm

As explained in Section 2.5, space decomposition meshing methods use a regular grid (in our case, an octree) over the domain to be meshed. The regular grid is larger than the domain, so at some point of the algorithm, this family of methods are forced to use a strategy to know if a given position of the grid is inside or outside the domain.

The coloring operation consists in determining where the entities are in the topological sense. This is, to determine whether each entity is inside a volume, out of the domain or laying on an interface entity between volumes. Applied to points, this operation is known in the literature as the point-in-polygon (PIP) problem [47,48,49].

In the presented meshing method, the coloring operation is applied to nodes (the case explained in this section) and to tetrahedra (as explained in Section 5.3.7).

This is one of the key points of the meshing algorithm, as it is not obvious to determine if a point in the space is inside a volume or not when the contours of the volume are non-watertight. Actually, if the contour of the volume has gaps, the concept of interior or exterior of the volume is not even defined from the topological point of view.

It has to be noted that several existing octree-based methods are focus on meshing a domain formed by a single volume. In this work arbitrary domains with several volumes are meshed at once, so coloring a node is not reduced to identify if it is inside or outside the domain: it has to be determined inside which volume (or interface entity) is.

## 4.1 Coloring strategies

There are several ways of coloring points considering a watertight definition of the volumes. These cases have some clear advantages as the contours of each volume can be oriented coherently (towards the interior or the exterior of the volume). This orientation provides with valuable information at the time of determining if a point near the contour entity is inside or outside the volume.

However, this work is focused on non-watertight geometries. This implies that a coherent orientation of the contours of a volume cannot be guaranteed. To deal with such geometries it is common to work with a voxelization of the model  [50,51,52]. These strategies aim to converting the non-watertight geometries into water-tight ones in order to be able to apply the coloring strategies of points. A voxel representation of a model is a regular grid (typically axis aligned and isotropic) where each voxel contains a topological information. In the case of study, it contains the color of the voxel. If a voxel is intersected by a surface entity it has the color of that surface (lets call it a contour voxel), otherwise it has the color of the volume it is into (inner voxel), or the color of the outer part of the domain if it is not inside any volume (outer voxel). All the points inside an inner voxel can be considered as interior to the corresponding volume. An example of voxelized 2D model is shown in Figure 26.

 (a) (b) (c)
Figure 26: Example of three voxelizations of a 2D model. The model has a gap in its contour and it is represented by the black curved line. Contour voxels are drawn in red. (a) The size of the voxels is large enough to close the gap of the domain: the topology of the voxelized model is watertight. (b) The size of the voxels is too small, so the gap of the domain is not closed. (c) The size of the voxels is too large: the gap is closed, but the final topology does not represent correctly the domain.

The advantage of working with voxelized models is that, depending on the size of the voxel, the topology of non-watertight geometries can be improved. If the gaps of the input boundary or the distance between overlapped contour entities is lower than the voxel size, the voxeled model may close the gaps or join the overlapped entities. This situation is shown in Figure 26(a), where the voxelized model presents a closed set of inner voxels (with no gap in its contour). As a drawback, the voxelized model can neglect some important topological information of the domain: the topology of the voxelized model shown in Figure 26(b) represents exactly the topology of the model (with its gap). Despite the size of voxels used in Figure 26(c) is large enough to close the gap (which is desirable), it represents an undesired topology, as it includes different sets of unconnected inner voxels.

The problem is that the voxel size needed to represent correctly the topology of the model cannot be estimated a priori and it may not be valid for the whole domain. Figure 26 shows that three different voxel sizes lead to three different topologies of the same domain.

If the voxelized model is watertight, one strategy for the coloring of the voxels is by propagation [50]. This strategy consists in the following steps:

• Determine all the contour voxels (the ones colliding a contour entity).
• Determine one (and only one) voxel interior to each volume solving the corresponding PIP problem.
• From a known inner voxel, assign its color to its neighbour voxels if they are not contour ones.
• Repeat the last step until all the voxels have color.

This propagation method is totally robust in watertight representations and has the advantage of solving very few times (one for each volume) the PIP problem, which can be computationally expensive. However, it has the drawback that is not naturally parallelizable. Only the voxels which are neighbour of a determined voxel (a voxel with a known color) can be colored at a time.

Voxelizing the model can help the coloring strategy as it can make the model watertight, enabling the use of the propagation method (which requires the solving of PIP problem only in a few points). However, it has been seen that even the voxelized model may be non-watertight. In this case the color propagation cannot be done, so the PIP problem has to be solved for each one of the points to be colored (in our case, the octree nodes).

### 4.1.1 Ray casting

The solution chosen in this work for coloring the octree nodes is based on the ray casting technique. The ray casting algorithm was first developed by Arthur Appel for rendering purposes in 1968 [53]. It proposes a solution for determining the visibility of a 3D object from a given point of view, and uses this information to paint a representation of the 3D object in a 2D image (made of pixels). The idea behind ray casting is to shoot rays (straight lines) from the point of view (one per pixel) and find the closest part of the object intersecting the ray. The algorithm needs to compute first all the intersections between a ray and the contours of the object (surface entities), and then get the closest to the point of view.

 Figure 27: 2D example of the ray casting technique to solve the PIP problem. Point ${\displaystyle P_{1}}$ is considered outside of the polygon because ray ${\displaystyle r_{1}}$ has an even number of intersection points (2). Points ${\displaystyle P_{2}}$ and ${\displaystyle P_{3}}$ are considered inside of the polygon because rays ${\displaystyle r_{2}}$ and ${\displaystyle r_{3}}$ have an odd number of intersection points (1 and 3 respectively).

Despite the first applications of ray casting were focused on rendering of 3D objects, its use has been generalized for several purposes following the philosophy of analyzing the intersection points between the ray and given 3D objects.

A common use of the technique is to solve the PIP problem. Among the several existing methods to solve this problem [47,48,49], a classical approach is the ray intersection method [47]. It consists in tracing a random ray from the point of analysis and compute the number of intersections between it and the contours surface entities. If the number is even, the point is outside the domain, and if it is odd, it is inside.

A 2D example illustrating the application of ray casting technique to solve the PIP problem is depicted in Figure 27. In this example, ${\textstyle P_{1}}$ is considered outside of the polygon because the ray ${\textstyle r_{1}}$ has two intersection points (an even number). Points ${\textstyle P_{2}}$ and ${\textstyle P_{3}}$ are considered inside of the polygon because their rays (${\textstyle r_{2}}$ and ${\textstyle r_{3}}$) have an odd number of intersection points (1 and 3 respectively).

## 4.2 Ray casting-based proposed technique

As pointed out previously, in the present work the coloring process of a node must identify not only if it is inside or outside the domain, but also the specific volume of the domain it is into. To fit this requirement, the ray casting method in this work focuses more on the topological information of each intersection rather than in the number of them. Following this philosophy, the proposed method consists in tracing rays in space and take care about the intersection of these rays with the contours of the volumes: if a ray intersects the interface between volumes A and B, it can be ensured that, near the intersection point, at one side of the intersection point the ray will be inside A, whereas at the other side it will be inside B. Following this principle, if we move along a ray from a point which color is known, we can color all the points of the ray just looking at the intersections of it with the contours of the volumes. Figure 28 shows a 2D example where a ray beginning in the outer part of the domain is colored in different parts corresponding to the contours it is intersecting.

 Figure 28: A 2D example of coloring parts of a ray (represented by the black arrow). The black dot is the beginning of the ray which color is known: ${\displaystyle 0}$. The crosses are the intersection points between the ray and the contours of the surfaces A and B (which are the ones forming the domain). In the upper line the color of the different parts of the ray is shown.

It is important to note that the presented technique only need the topological information of the surface entities: which volumes they are intersecting. Other techniques require the use of normal vectors in some points, or a given connectivity among the surface entities defining the contours of the domain. This is crucial considering the proposed algorithm must work with non-watertight definitions of the volumes. In these cases it is not obvious to define the oriented normal of a surface entity pointing towards the inner or the outer part of a volume.

### 4.2.1 Pathological configurations

The ray casting technique presents some pathological configurations [49] for the PIP problem. One of these is the case when the intersection between the ray and the volume boundaries is done tangentially (T type intersection defined in Section 3.4). In this situation the ray intersects the boundaries of the geometry, but both sides of the intersection are in the same volume. To solve this problem, this kind of intersections (tangentially to the boundaries of the volumes) are not taken into account for coloring the regions of the ray. Figure 29(a) shows an example of this pathological case: both sides of the ray from the dark intersection point are inside surface B, although there is an intersection point.

Another pathological configuration occurs when the intersection point is a point of the boundary interfacing more than two volumes (in 2D case, more than two surfaces). This is the case shown in Figure 29(b). For notation purposes, this kind of intersections are referred as a M intersection type.

If a ${\textstyle M}$ intersection is detected in a ray, a color has to be chosen for the following part of the ray (from the ${\textstyle M}$ intersection point on). The strategy followed to decide this color is explained in the Implementation chapter (Section 4.3.1). It is based on a try and error approach considering all the possible colors.

 (a) (b) Figure 29: 2D examples of pathological configurations for the ray casting technique. (a) Both sides of the ray from the dark cross intersection are colored equally (surface B) although there is an intersection point because the ray is tangent to the contours in it. (b) M point is interface between surface A, B and ${\displaystyle 0}$ (outer part of the domain), so the color of the right part of the ray cannot be set.

The case of co-planar intersections (${\textstyle PT}$ and ${\textstyle PI}$ intersection types defined in Section 3.4) presents also a pathological configuration for the ray casting technique. In these cases, the part of the ray co-planar with the boundaries must be colored with the color of the interface itself. In this case, the color does not correspond to a volume, but to an interface between volumes. 2D examples for ${\textstyle PT}$ and ${\textstyle PI}$ intersection types between the ray and the boundaries of a model is depicted in Figure 30. The end point of the co-planar part of the ray is considered as an ${\textstyle M}$ point from the coloring point of view.

 (a) (b) Figure 30: 2D examples of co-planar intersections of the ray. (a) ${\displaystyle PT}$ intersection type. The color at the right part of the ${\displaystyle M}$ point could be ${\displaystyle A}$ or zero. (b) ${\displaystyle PI}$ intersection type. The color at the right part of the ${\displaystyle M}$ point could be ${\displaystyle A}$, ${\displaystyle B}$ or zero.

For non-watertight contours of the domain volumes, two more pathological configurations affect the ray casting technique. They are related to the possible gaps and overlapping entities in the contours (G and W intersection types defined in Section 3.4):

• If a ray passes through a gap of the contours of a volume (G type intersection), it does not detect that is entering into the volume, so it would assign the following regions an invalid color. This is the case shown in Figure 31(a), where the ray casting technique only identifies two zones in the ray separated by the intersection point: the left region as outer, and the right one as inside surface A.
• When a part of the contour of a volume is overlapped (W type intersection), if a ray intersects the contour by that region, it intersects more than one time the contour. The case where the distance ${\textstyle d}$ between the overlapping entities is smaller than ${\textstyle tol_{c}}$ presents no problem, as only one intersection point is considered: the MIP. This can be understood as a local voxelization of the model, as the topological information of different entities is collapsed into one data (voxel). The case where ${\textstyle d}$ is greater than ${\textstyle tol_{c}}$ disturbs the coloring operation as the color of the ray changes two times instead of one. This is the case shown in Figure 31(b), where the part of the ray inner to surface A is colored as outer because of the two close intersections highlighted in the figure.
 (a) (b) Figure 31: 2D examples of pathological configurations for the ray casting technique with non-watertight domains. (a) G type intersection of the ray. The part of the ray inside the surface A is not colored as A because the ray has not intersected its boundary entities. (b) W type intersection of the ray. The part of the ray inside surface A is colored as zero because the ray intersects two times the boundary in the same region (crosses).

Both pathological cases (gaps and overlappings greater than ${\textstyle tol_{c}}$) can lead to set the ray as invalid. An invalid ray is a ray with some coloring contradictions. There are two kinds of contradictions, depicted in Figure 32(a) and (b):

• The last part of the ray (from the last intersection point to the outer part of the bounding box of the model) has a color different than zero (see Figure 32(a)). Considering the principle that the parts of the ray outside the bounding box of the domain must be outer, the beginning and end of the ray must be outer (color zero).
• A part of the ray colored as ${\textstyle i}$ ends in an intersection point with a surface entity not interfacing the ith volume (Figure 32(b)).
 (a) (b) (c)
Figure 32: Types of invalid rays (crosses are intersection points). (a) The last part of the ray is colored as ${\displaystyle A}$, but it should be zero.(b) At the dark cross intersection point the ray should arrive with color ${\displaystyle A}$ or ${\displaystyle B}$, as these are its interfacing surfaces, but its color is zero. (c) Rays to be canceled because they provide with different colors to a given position.

Situations of invalid rays can also happen even if the contours of the domain are watertight, taking into account that the intersection operations in 3D are done numerically, so they depend on tolerances. Some specific configurations of the ray and the contours may lead to an invalid ray.

### 4.2.2 Adaptations of the method

A special kind of rays are the Cartesian rays, which are straight and parallel to the ${\textstyle x}$, ${\textstyle y}$ and ${\textstyle z}$ directions. The ray casting technique can be applied to any kind of ray, but in the proposed algorithm Cartesian rays are used because they simplify the intersection operations, and they take profit on the spacial distribution of the nodes provided by the octree structure (this aspect will be seen in more detail in Section 4.3).

Taking into account that the domain to be meshed is totally inside its bounding box, the octree nodes which are outside it can be directly colored as zero (they are out of the domain) without the need of any coloring process. All the other octree nodes are colored using rays which begin out of the model bounding box, so as the color of the initial point of the ray is known: zero. The color of a node is directly the color of the part of the ray the node is into.

As there are situations in which a ray is considered as invalid to color a point because of a pathological configuration, the following strategy is applied in order to color all the octree nodes:

• For each octree node, three Cartesian rays (one for each direction ${\textstyle x}$, ${\textstyle y}$ and ${\textstyle z}$, from the lower to the higher coordinate) passing by the position of the node are shooted from the outer part of the bounding box of the domain.
• Each one of these three rays can be invalid or valid. If it is valid it provides with a candidate color for the octree node. This opens the following cases:
• There is only one valid ray. The node is colored with the color given by the ray.
• There is more than one valid ray and they give the same color. The node is colored with that color.
• There is more than one valid ray and they give different colors. A 2D example of this case is depicted in Figure 32(c). In this case the three rays are canceled. To cancel a ray means not to use its color information. It has to be taken into account that a ray can be used to color more than one node, as some nodes may be aligned with the Cartesian directions. In this case, canceling a ray also affects the other nodes it passes throw.
• The three rays are invalid. In this case the node is not colored.
• At this point of the process there may be nodes with no color. These ones are colored using a local ray casting. This operation consists in launch a ray (not a Cartesian one) from a neighbor node with a known color to the node itself. The color is set evaluating the intersections of this local ray with the input boundaries. A 2D example of this local ray casting operation is shown in Figure 33. It has to be noted that the situation where the node cannot be colored with the Cartesian rays is rather improbable, as it should imply the coincidence of pathological situations in the three directions of the space aligned with the node.
 Figure 33: A 2D example of local ray casting. A non watertight representation of surface A is shown. The two Cartesian rays passing by the black dot (drawn with dotted arrows) are not valid. The color of the white dot is known (inside surface A). As the ray from the white dot to the black one has no intersection with the boundaries of the domain, the color of the black dot is also set to A.

This strategy have been proved in several examples, and it solves successfully the coloring process for the octree nodes with pathological configurations for the ray casting technique. The use of three rays for each node provides with a redundancy minimizing the possible effect of an invalid undetected ray.

## 4.3 Implementation of nodes coloring algorithm

This section is focused in the implementation of the algorithm explained in Section 4.2.2 for nodes coloring. It requires three Cartesian rays for each one of the nodes in order to color them (determine where they are topologically). Applying this algorithm to a general cloud of ${\textstyle n}$ nodes would imply to compute the intersections of ${\textstyle 3\cdot n}$ rays. However, as the nodes of the tetrahedra mesh treated in this work are octree nodes, most of them are placed on regular positions in space, aligned with the Cartesian directions. To take advantage on this configuration, one ray can be used in the coloring of several nodes: all the ones aligned with a Cartesian direction. Following this strategy an important time saving is achieved in the nodes coloring operation.

It is important to note that only the nodes with unknown color must be colored following the presented algorithm. The octree nodes with known color before the coloring operation are:

• The forced nodes (section 5.3.2); as they are linked to an entity coming from the input data, they directly have the color of the corresponding input geometrical entity.
• The octree nodes outside the bounding box of the domain. These ones are outer nodes, so their color is zero.

A 2D example of the rays used to color the nodes of a given configuration of a quadtree is shown in Figure 34. In the example, there are ${\textstyle 28}$ nodes to be colored and ${\textstyle 22}$ rays are used (${\textstyle 12}$ for ${\textstyle x}$ direction and ${\textstyle 10}$ for ${\textstyle y}$ direction).

 (a) (b) (c)
Figure 34: Rays used to color the quadtree nodes of a 2D example. The contour of the domain is the solid black line and the dotted line represents its bounding box. Rays are represented with red arrows. (a) Quadtree configuration with all the octree nodes represented. The black nodes are forced nodes. (b) The rays used in ${\displaystyle x}$ direction to color the nodes. Only the nodes to be colored are plotted in order to clarify the figure. (c) The rays used in ${\displaystyle y}$ direction.

The implementation of the coloring operation is based in the following steps:

• Creation of all the rays. For each one of the three directions of the space (${\textstyle x}$, ${\textstyle y}$ and ${\textstyle z}$), the coordinate of each one of the nodes to be colored is considered. Then, one ray is created for each one of the coordinates. This process generates a collection of rays. Each one of the rays has the following information:
• The alignment direction.
• The two coordinates corresponding to the plane perpendicular to its direction.
• The nodes it passes through, sorted from lower to upper coordinate in the direction of the ray.
• Color the rays and set the contribution to each node. In this step each one of the rays is colored. It has to be noted that this process is fully parallel, as each ray is totally independent from the others. In this process all the intersections between the ray and the input boundaries are computed. They are sorted from lower to upper coordinate, dividing the ray in different parts corresponding to the regions between intersection points. Looking at the interfaces intersecting at each intersection point the color of the parts of the ray can be set. Pathological situations in the intersection operations are detailed in Section 4.3.1. As explained in Section 4.2.2, each node has three possible colors, corresponding to the three Cartesian rays passing through it. At this point, the color contribution of the ray is set to all the nodes it passes through.
• Determine the invalid rays. As explained in Section 4.3.1, the invalid rays are identified. Then, the rays with contradictions in one node are canceled.
• Set the color of the nodes. The process described in Section 4.2.2 is used to determine the color of the node considering its three color contributions: if there is some valid ray the color is set automatically, otherwise, the local ray casting is performed. In case the local ray casting provides with different candidate colors (from the local rays from different neighbors), a voting process is performed: the color set is the one given by more neighbors.

### 4.3.1 Pathological situations

In this section the strategy followed to solve the pathological situations for the ray casting technique defined in Section 4.2 is detailed. These situations are related with the ${\textstyle T}$, ${\textstyle C}$, ${\textstyle P}$, ${\textstyle W}$ and ${\textstyle G}$ intersection types and the ${\textstyle M}$ intersection points.

In all these situations the intersection point is a ${\textstyle MIP}$ (Section 3.4). Actually, the ${\textstyle G}$ intersection type (situation where the ray passes through a gap in the boundaries) involves a ${\textstyle GIP}$, but this kind of intersections is not solved in the present approach because trying to find the corresponding ${\textstyle GIP}$ to all the parts of the ray not intersecting the input boundaries would be too computationally expensive. Leaving this situation as unsolved could lead to rays declared as valid, but giving wrong information in terms of coloring. It has to be noted that if this situation occurs, it does not mean necessary that an invalid color will be assigned to the nodes. As the three Cartesian rays passing by a node contribute to the decision of its color, only in cases where the three rays are considered valid, provides with wrong information and this information is compatible (within the three rays), the color of the node would be wrong. This is a rather improbable situation, which have never happened in the examples run in this work.

A MIP in a ray indicates the presence of a pathological configuration. Each intersection point has the information of the two volumes interfaced by the surface entity where it is (the interfacing volumes of the intersection point). The interfacing volumes of a MIP are the union of the interfacing volumes of each intersection point involved, so a MIP can have two or more interfacing volumes.

In Section 4.2 the situation of ${\textstyle P}$ intersection types (co-planar ones) is explained. They correspond to a ${\textstyle M}$ intersection point at the end part of the co-planar part of the ray. Considering that the presented algorithm uses linear triangle elements as surface entities for defining the boundaries of the domain, the situation where the ray is co-planar can be detected easily. If geometrical surface representations (line NURBS) are used for defining the contours of the volumes, solve this situations becomes much more difficult. For coloring purposes, the co-planar part of the ray takes no relevance, as it is ensured there will be no node to be colored there. If a node is so close to a surface entity it is a forced node in interface (Section 5.3.2), so its color is already set without the need of ray casting technique. The problematic situation comes at the time to color the following part of the ray (next to the co-planar one), where a color have to be assigned among the corresponding interfacing volumes.

When a MIP (a pathological configuration) is detected, its type has to be determined (${\textstyle T}$, ${\textstyle C}$, ${\textstyle W}$ or ${\textstyle M}$) in order to apply the right strategy to color the following parts of the ray from that point on. For this purpose, a sort of auxiliary local rays are built: the surrounding segments. These are segments parallel to the ray at a distance ${\textstyle tol_{c}}$ of it. In the present implementation, a number of ${\textstyle 4}$ surrounding segments has been chosen following the Cartesian directions. Longitudinally, the surrounding segments are centered in the MIP position and have a length of ${\textstyle D}$, being ${\textstyle D}$ the distance to the closest intersection point of MIP (in the ray). An example of the surrounding segments of a ray is shown in Figure 35. In this example the ray intersects two triangles by its contour edge, so it is a ${\textstyle C}$ intersection type.

 Figure 35: Surrounding segments of a ray around a MIP. The two triangles represent a part of the input boundaries. The black dot is the MIP and the white ones are the intersection points between the surrounding segments and the input boundaries. The black cross is the nearest intersection point of the ray to the MIP.

Then, the intersection points between the surrounding segments and the input boundaries are computed. Considering them and the MIP interfacing volumes, the following cases should be accounted for distinguishing the type of pathological configuration:

• Case I: there are more than two interfacing volumes among the ones of the MIP and the ones of the intersection points of the surrounding segments. This case corresponds to the ${\textstyle M}$ intersection type.
• Case II: case ${\textstyle I}$ is not accomplished (so only two interfacing volumes are involved) and some surrounding segment has no intersection point. This case corresponds to the ${\textstyle T}$ intersection type.
• Case III: when cases ${\textstyle I}$ and ${\textstyle II}$ are not accomplished. This means that all the surrounding segments have intersection points, and there are only two interfacing volumes involved. The following sub-cases are considered:
• Case IIIa: each one of the surrounding segments have only one intersection point (MIP or single one). This case corresponds to the ${\textstyle C}$ intersection type.
• Case IIIb: all the surrounding segments have only one intersection point and some of them is a MIP (this means all the intersection points of the surrounding segment are enclosed in a distance ${\textstyle tol_{c}}$). This case corresponds to the ${\textstyle W}$ intersection type.
• Case IIIc: there are surrounding segments with more than one intersection point. This case corresponds to the ${\textstyle T}$ intersection type.

Actually, cases IIIa and IIIb could be merged, as case IIIb is included in IIIa. The differentiation has been made only to highlight that a ${\textstyle W}$ intersection type always implies a MIP, and the ${\textstyle C}$ one can present a MIP or not in the surrounding segments. This is not relevant because the treatment of ${\textstyle C}$ and ${\textstyle W}$ intersection types is analogous for coloring purposes.

A graphical interpretation of these cases in a 2D example is illustrated in Table 4.3.1. As a 2D cases, only two surrounding segments take sense for a MIP.

 Case I: there are three interfacing surfaces (${\textstyle A}$, ${\textstyle B}$ and zero) involved in the intersection points of the surrounding segments. ${\textstyle M}$ intersection type. Case II: all the intersected entities have the same interfacing surfaces (${\textstyle A}$ and zero), and there is one surrounding segment with no intersection. ${\textstyle T}$ intersection type. Case IIIa: each of the surrounding segments intersects the boundaries in one point and all the interfacing surfaces are the same (${\textstyle A}$ and zero). ${\textstyle C}$ intersection type. Case IIIb: all the surrounding segments can create a MIP (its intersection points are closer than ${\textstyle tol_{c}}$) and all the interfacing surfaces are the same (${\textstyle A}$ and zero). ${\textstyle W}$ intersection type. Case IIIc: Some surrounding segment has more than one intersection point. ${\textstyle T}$ intersection type.

Different cases for surrounding segments for identification of intersection types corresponding to a MIP (white dot) in a 2D example.

Once the type of pathological intersection has been determined, the following strategy is carried out to color the following part of the ray (the one next to the MIP):

• ${\textstyle T}$ type intersection. In this case the color of the ray is not changed in the following ray part.
• ${\textstyle C}$ and ${\textstyle W}$ type intersections. In this case the color of the ray is changed accordingly with the interfacing volumes: as there are two interfacing volumes (${\textstyle A}$ and ${\textstyle B}$), if the previous part of the ray was colored as ${\textstyle A}$, the next one is ${\textstyle B}$ and viceversa.
• ${\textstyle M}$ type intersection. This case is the more complicated one. Considering the ${\textstyle n}$ interfacing volumes of the corresponding MIP and the color of the previous part of the ray, the following part of it can have any color among the rest ${\textstyle n-1}$ colors. The strategy chosen is based on a try and error approach: ${\textstyle n-1}$ rays are shooted from the MIP, one for each possible color, and if one of them is valid, it is taken as the good one. If no one of them is considered valid, then the whole ray is set as invalid.

There is one situation not solved by the presented strategy. It is a mix between the ${\textstyle T}$ and the ${\textstyle M}$ intersection types. It occurs when the ray passes tangentially to two parts of the boundaries of a volume at the same point. A 2D example of this situation is shown in Figure 36. If the corresponding surrounding segments have one MIP each one, it corresponds to the case IIIb, so the intersection would be considered as ${\textstyle W}$ and the color of the ray would change when it should remain with the same color ${\textstyle A}$ (as if it was a ${\textstyle T}$ intersection).

 Figure 36: Pathological situation where the analysis of the surrounding segments can lead to erroneous classification of intersection type if ${\displaystyle d.

### 4.3.2 Computation of distances

As explained in previous sections, this work proposes a meshing algorithm for body fitted meshes and for embedded ones. The later ones are used in embedded and immersed methods, and require for the nodes of the mesh their distance to the input boundaries in addition to their color (Section 5.2). This section focuses in the computation of these distances for the embedded meshing algorithm.

As explained in Section 5.2.1, the computation of distances from the nodes of the tetrahedra mesh to the input boundaries for embedded methods takes advantage on the ray casting technique used in the coloring algorithm. This is mainly because the Manhattan distance [54] has been chosen for approximating the distance for far nodes (nodes further than a given distance of the boundaries of the domain). The Manhattan distance between two points is the sum of the distances of the three Cartesian components of the vector defined by the points. As the rays used in the coloring algorithm are Cartesian, its use fits perfectly for this purpose.

To improve the efficiency in a parallel implementation of the method, each ray contributes with a distance to the nodes it passes through. This implies each node has three contributions of distances, as each node has three rays passing trough it (one for each direction: ${\textstyle x}$, ${\textstyle y}$ and ${\textstyle z}$).

The process of computing the distances from the nodes of the tetrahedra mesh to the input boundaries follows the steps detailed hereafter:

1. Compute the exact distance to the input boundaries of the octree nodes belonging to an interface leaf of the octree. The exact distance is the one from the octree node to its closest coordinate onto the input boundaries. The strategy chosen to compute this distance profits from the octree itself. The surface entity from the input boundaries the closest coordinate is onto must collide an octree leaf containing the octree node. So only the surface entities colliding these octree leaves are considered at the time of finding the closest coordinate to a node.
2. Set to the maximum possible distance the three distance contributions of the rest of octree nodes.
3. Proceed with the ray casting technique corresponding to the coloring process. Apart from the color assignment to the nodes, an extra operation is performed in order to set the distances. Each ray launched in the coloring process has the information of its intersection points with the input boundaries, and the octree nodes it passes through (the ray nodes). Being ${\textstyle d_{0}}$ the existing distance contribution of a ray node and ${\textstyle d_{i}}$ the distance from the node to its closest intersection point in the ray, ${\textstyle d_{0}}$ is equalized to ${\textstyle d_{i}}$ if ${\textstyle d_{i}.
4. Propagate the distance of each node within each ray. This step only affects the rays with more than one node. Let us consider the nodes of a ray sorted out by the coordinate corresponding to the direction of the ray. Being ${\textstyle n}$ one of the nodes and ${\textstyle n_{-1}}$ and ${\textstyle n_{+1}}$ its previous and next nodes, ${\textstyle d_{n}}$, ${\textstyle d_{n-1}}$ and ${\textstyle d_{n+1}}$ are defined as their corresponding distance contributions of the ray. (For the first node only the ${\textstyle n_{+1}}$ is considered, and for the last one only the ${\textstyle n_{-1}}$). Two propagated distances are defined for each node:
5.  ${\displaystyle pd^{-}(n)=d_{n-1}+d(n_{-1},n)}$
(4.1)
 ${\displaystyle pd^{+}(n)=d_{n+1}+d(n,n_{+1})}$
(4.2)

where ${\textstyle d(a,b)}$ is the distance between nodes ${\textstyle a}$ and ${\textstyle b}$. Note that, as the nodes are aligned in the ray, the distance between them is directly the difference between their coordinates corresponding to the ray direction.

Considering each node ${\textstyle n}$ of the ray, its distance contribution is set to the minimum between ${\textstyle d_{n}}$, ${\textstyle pd^{-}(n)}$ and ${\textstyle pd^{+}(n)}$.

6. Update the distance contributions of each node. This operation consists in equalizing the three distance contributions of each node to the minimum one of them.
7. Repeat Steps 4 and 5 (propagate and update) two more times. This two steps must be done a number of times equal to the dimension of the problem (in our case, three times). This is to allow a complete propagation of the distances from the nodes. As this propagation is done using Cartesian rays, the information provided by one direction is not passed to the other ones until the distance contributions are updated. A 2D example is depicted in Figure 37 where the iso-lines for distance are plot before and after the second propagate and update operation (as a 2D case, propagate and update operation must be done only two times). It can be seen (in Figure 37(b)), that the distances in the lower left and right part of the surface are not well computed. As the hole of the surface is not aligned with them, the second propagation and update operation is needed to compute correctly the distances to the boundaries of all the nodes (Figure 37(c)).
 (a) (b) (c)
Figure 37: 2D example for computation of distances to the boundaries. (a) Surface representing the domain. The surface is a square with a squared hole inside. (b) Iso-lines of distance to the boundaries of the surface when only one propagation and update operation have been done. (c) Iso-lines of distance performing the propagation and update operation two times.

# 5 Octree based mesher

## 5.1 Introduction

This chapter present the octree-based mesher developed in this monography.

One of the requirements defined in Section 1.2.3 refers to the immersed [6] and embedded [7] methods. This family of methods uses the so called embedded meshes. An embedded mesh is a mesh of a part of the 3D space containing in its interior the domain to be simulated, but the mesh is not limited to the inner part of the domain and, furthermore, does not fit the contours of it. The effect of the contours of the domain in the results of the numerical simulation to be run is captured by the assignment of special boundary conditions in the mesh elements and nodes. Because of its nature, it takes no sense to talk about preserving features or topology preservation for embedded mesh generation.

On the other hand, there are the body-fitted meshes, which must capture precisely all the contours of the domain. They need a faithful representation of the boundaries preserving the geometrical features, the topology of the volumes of the domain, and forcing the nodes and elements to follow the shape of their contours.

The algorithm proposed in this monography considers the embedded meshing as a particular case of mesh generation where there is no need to apply any strategy for preserving geometrical features or volume topology. In this sense, the generation of embedded meshes will be easier than the generation of body-fitted ones.

The detailed explanation of the meshing algorithm requires the definitions of specific concepts and auxiliary algorithms explained in Chapter 3. In Sections 5.2 and 5.3 the meshing algorithm itself is presented for embedded and body-fitted meshes.

As a general view, the main steps of the mesher are pointed out hereafter. Let us consider that the input data for the mesher is the definition of the contours of the volumes of the domain to be meshed (from now on, input boundaries) and a given list of parameters (detailed in Section 3.1). The new embedded meshing algorithm can be summarized in five steps:

1. Create an octree enclosing the domain.
2. Refine the octree following some given criteria.
3. Generate a tetrahedra mesh from predefined patterns from the octree.
4. Find the volume containing each node of the mesh.
5. Compute the distance of the nodes of the resulting mesh to the input boundaries.

The first three steps are directly related to the octree structure. The two other ones involve mesh operations that can be applied to any unstructured tetrahedra mesh independently on its generation method, although the octree structure can be used as an auxiliary tool in order to improve the efficiency of the algorithms. They are detailed in Chapter 4 and Section 5.2.1.

Concerning the body-fitted mesher, the four first steps are the same ones as for the embedded mesher (adding some extra refinement criteria to refine the octree in the second step), and the fifth step is replaced by the following three ones:

• Apply given mesh edition operations to fit into the tetrahedra mesh the sharp edges and corners to be preserved and the surface entities of the domain. These mesh edition operations are detailed in Sections 5.3.5 and 5.3.6, and they basically involve nodes movement and element splitting.
• Identify the volume each tetrahedra of the mesh is into. This part of the algorithm is detailed in Section 5.3.7.
• Improve the mesh quality if needed using make-up and smoothing operations. This step is optional, as it is only needed if the mesh quality is poor. It is detailed in Section 5.3.8.

All these specific steps related to the body-fitted mesher are also applicable to any unstructured mesh, independently on the generation method chosen.

An important aspect of the new mesher is that it generates the mesh of the whole domain at once, including all the volumes which are part of it. Other existing meshers are designed to generate the mesh of just one volume, so they generate the mesh of the whole domain on a volume by volume manner.

## 5.2 Embedded mesher

One of the objectives of the presented meshing algorithm is to generate meshes for embedded methods. This section is focused in the specific aspects of the algorithm for this kind of methods.

The main characteristic of embedded methods is that the mesh used is not body-fitted, and its nodes have the information of:

• where they are topologically (inside which volume or onto which interface between volumes);
• the distance they have to the contours of the domain.

Knowing if a node is inside or outside a volume is solved by the proposed coloring algorithm (Chapter 4), and the computation of distances also takes profit of that algorithm, as explained in the following section. Furthermore, it is common to apply these methods to domains where only one volume is involved, so only two colors take part: inside and outside.

Embedded meshes are mainly used in Computational Fluid Dynamics (CFD) simulations, where the behaviour of a fluid around solid bodies is studied. If the solid bodies are in movement, rather than remeshing the whole domain at each time step, this family of methods maintains the volume mesh static and updates the nodes information (color and distance).

### 5.2.1 Computation of distances

The distance function plays a key role in embedded methodology in order to apply the boundary conditions. It is common to combine the coloring and distance function into one signed distance function where negative distance indicates inside and positive distance indicates outside nodes. A direct result of this definition is the fact that the iso-surface representing the zero distance defines the approximated embedded boundary condition, as shown in Figure 38.

 (a) (b) Figure 38: (a) Surface entities defining a cube. (b) The cube representation via a zero iso-surface of the distance function in the tetrahedral embedded mesh generated.

To ensure the accuracy of the method the distance near the boundary and especially in cut elements must be calculated exactly. For the nodes far from the boundary the exactness of the distance becomes less important so, in order to reduce the computational cost, it is convenient to calculate the exact distance only for the points inside a given distance range from the boundary, and leave the rest with a maximum value (an upper distance limit). However, in cases with moving boundaries one may convect the distance function given the interface motion. In this procedure, having sharp gradients in distance functions can lead to numerical error and having a constant maximum distance near the boundary may affect the convergence and results. In order to deal with this problem an approximation of the distance function from a given distance of the boundary would be interesting.

In this work the exact Euclidean distance is calculated for the octree nodes belonging to the interface leaves. For the rest of them the Manhattan distance [54] is computed. It is an estimation of the exact distance given by Equation 5.1:

 ${\displaystyle d_{M}(\mathbf {p} ,\mathbf {q} )=\sum _{i=1}^{n}\parallel \mathbf {p} (i)-\mathbf {q} (i)\parallel }$
(5.1)

where ${\textstyle d_{M}(\mathbf {p} ,\mathbf {q} )}$ is the Manhattan distance between positions ${\textstyle \mathbf {p} }$ and ${\textstyle \mathbf {q} }$, and ${\textstyle i}$ represents the ${\textstyle i}$-coordinate of the vector. ${\textstyle n}$ is the dimension of the space considered, in this case: three.

The Manhattan distance provides with the desired continuity and smoothness of the distance function in order to deal with moving objects with minimum distance calculation overhead. The reason to use this distance measure is to take profit on the coloring algorithm, as explained in Section 4.3. Cartesian rays are used for the coloring of the nodes, so using the distances onto the rays provides directly with the Manhattan distance.

### 5.2.2 Octree refinement criteria for embedded meshes

As explained in Section 5.1, one of the main steps of the meshing process is refine the octree following given criteria. The basic idea is that the octree must be refined in such a way that the application of the further steps of the meshing process ensures the accomplishment of the requirements to be covered.

In this section all the refinement criteria (RC) needed for embedded meshes are defined. It is important to highlight that all the octree refinement criteria defined in this section (the ones needed for embedded meshes) only depend on the desired mesh sizes and the topology of the octree itself (the balance criterion). This implies that the application of these refinement criteria does not need to consider the input boundaries. For notation purposes, the collection of all the refinement criteria defined in this section is called size refinement criteria.

#### 5.2.2.1 Uniform size refinement criterion

Let us consider the bounding box of the model (${\textstyle Bbox_{m}}$). It is a parallelepiped with its faces parallel to the Cartesian axes (note that ${\textstyle Bbox_{m}}$ is not needed to be coincident with the octree root). ${\textstyle s_{box}}$ is defined as the length of the smaller side of ${\textstyle Bbox_{m}}$.

For given desired mesh sizes (via mesh size points or general mesh size parameter), lets define ${\textstyle s_{dmax}}$ as the maximum of the desired mesh sizes entered in the input data. Then, ${\textstyle S_{max}}$ is defined by Equation 5.2:

 ${\displaystyle S_{max}=min(s_{dmax},s_{box})}$
(5.2)

If no size has been introduced in the input data, then ${\textstyle S_{max}}$ is directly ${\textstyle s_{box}}$.

It has to be noted that because of the tetrahedra pattern definition (Section 3.3.3), tetrahedra sizes are directly related to the octree cells sizes. Actually, the size of the edges of the tetrahedra generated from a cell is always equal or smaller than the cell size. Refining the octree implies dividing by two some of the cells, reducing accordingly the corresponding tetrahedra size. This aspect makes impossible to reach exactly a given size for a tetrahedron. The algorithm tends to fit it by subdividing cells until their size is close to the desired one approximating it above or below. To incorporate in the algorithm the difference between the desired mesh size and the cell size and avoid an excessive level of refinement, the parameter ${\textstyle \alpha _{ms/cs}}$ is defined. This is a real value ranging between 1 and 2. The value used for ${\textstyle \alpha _{ms/cs}}$ is specified in Section 6.7.

Considering the isotropic octree structure, the longest edge of the final mesh of the model must be always smaller than ${\textstyle S_{max}}$. This leads to the first refinement criterion:

RC 1: If an octree cell collides with ${\textstyle Bbox_{m}}$ and its size is greater than ${\textstyle (\alpha _{ms/cs}\cdot S_{max})}$, the cell must be subdivided.

#### 5.2.2.2 Mesh size entities refinement criterion

In order to account with the possible desired mesh sizes required by the simulation defined with the mesh size entities, the following refinement criterion is defined:

RC 2: If an octree cell collides with a mesh size entity with a desired size ${\textstyle s}$ and its size is greater than ${\textstyle (\alpha _{ms/cs}\cdot s)}$, the cell must be subdivided.

At this point the concept of generalized mesh size points needs to be introduced. To make easier the implementation of the method, and to automatize the methodology, the use of points instead of generic entities (lines, surface and volume ones) helps. The idea is to replace (only for mesh size purposes) each mesh size entity by a collection of mesh size points with the same desired mesh size: its generalized mesh size points. Actually only mesh size lines, surfaces and volumes are involved in this process, as the mesh size points are already points.

If we consider a mesh size entity with a desired mesh size ${\textstyle s}$, its generalized mesh size points are located onto the mesh size entity in such a way that there is no point onto the entity further than ${\textstyle s}$ from a generalized mesh size point. An example of the generalized mesh size points of a surface mesh size entity (a triangle) is shown in Figure 39.

 Figure 39: Representation of a surface mesh size entity (a triangle) with a set of its generalized mesh size points (black dots). It can be seen that all the coordinates inside the triangle have at least one mesh size point closer than ${\displaystyle s}$.

It has to be noted that, following this definition, there can be infinite sets of generalized mesh size points. The process used to obtain them is detailed in Section 6.3.

The implementation of the mesh size entities criterion is simplified if the generalized mesh size points are considered instead of the mesh size entities.

#### 5.2.2.3 Volume mesh sizes refinement criterion

For a given mesh size desired for each volume (considering ${\textstyle s_{vol,i}}$ as the desired size for the ith volume) the following criterion is defined:

RC 3: If an octree cell is an inner cell of the ith volume and its size is greater than ${\textstyle (\alpha _{ms/cs}\cdot s_{vol,i})}$, the cell must be subdivided.

#### 5.2.2.4 Balance refinement criterion

The following refinement criterion is widely used in octree based meshers. It limits the number of neighbors of one cell, and it is often referred as constrained two to one condition (Section 3.3). It is an essential criterion for the pattern used for the creation of tetrahedra from the octree leaves:

RC 4: If an octree cell has more than four neighbor cells by face or two by edge, the cell must be subdivided.

#### 5.2.2.5 Size transition refinement criterion

The balance criterion gives an upper limit for the size transition: as the size of two neighbor cells differs (at maximum) by a factor of 2, the difference between the size of the tetrahedra generated from those cells is bounded as well. However, it may be required for the final mesh to present a more smooth size transition between regions with small and large elements.

Considering the desired mesh size for given regions of the domain (from the input data), and a given size transition function, an envelope function can be defined to set an upper limit for the element size allowed in each position of the space ${\textstyle S_{u}({\vec {x}})}$. The definition if this function is detailed in Section 6.4, and leads to the following refinement criterion:

RC 5: Being ${\textstyle {\vec {c}}_{i}}$ the center of an octree cell, if the size of the cell is greater than ${\textstyle (\alpha _{ms/cs}\cdot S_{u}({\vec {c}}_{i}))}$, the cell must be subdivided.

### 5.2.3 Octree root

As explained in Section 3.2, the octree root is the bounding box of the octree. Clearly it should contain totally the domain to be meshed in its interior, but considering the meshing process, it has to accomplish some other requirements.

For the creation of the octree root, the extended bounding box of the model is needed (from now on ${\textstyle Bbox_{m}^{+}}$). ${\textstyle Bbox_{m}^{+}}$ is the minimum bounding box containing the bonding box of the model and all the generalized mesh size points. It has to be noted that the mesh size entities can be outside the domain. The octree root is built centered in the center of ${\textstyle Bbox_{m}^{+}}$, and with a size equal to the maximum size of ${\textstyle Bbox_{m}^{+}}$ plus ${\textstyle S_{max}}$ (Equation 5.2). This offset of the octree root with respect to ${\textstyle Bbox_{m}^{+}}$ is needed to ensure the tetrahedra creation by the tetrahedra pattern, considering that in some cases this creation involves one cell and its neighbor. A graphical 2D example (quadtree instead of octree) of the ${\textstyle Bbox_{m}^{+}}$ and the octree root is shown in Figure 40.

 Figure 40: 2D example where the domain is the solid surface. Black dots are the generalized mesh size points. The bounding box of the model is represented by the dotted line. The ${\displaystyle Bbox_{m}^{+}}$ is the gray line, and the quadtree root is the black square.

### 5.2.4 Meshing algorithm

Hereafter, the steps of the meshing algorithm for embedded meshes are presented.

1. Process input data and create the octree root. This includes:
• Creation of generalized mesh size points (in case they exist).
• Compute ${\textstyle S_{max}}$.
• Creation of the octree root (Section 5.2.3).
2. Refine the octree accomplishing the size refinement criteria (defined in Section 5.2.2):
• RC 1: Maximum size.
• RC 2: Mesh size entities.
• RC 3: Volume size.
• RC 4: Balance.
• RC 5: Sizes transition.
3. Classify the input boundary entities into the octree. This step implies the identification of the octree leaves colliding with the input boundary. Until now the octree has been refined, but each one of the cells had no information about its type: interface, inner or outer. Prior to this step the input boundaries only have been taken into account for building the bounding box of the model. At this point is where the input boundaries (surface entities from the input data) are analyzed in order to determine the interface cells. Furthermore (as it will be explained in Section 6.2) in the implemented octree, each interface cell stores the input boundary entities colliding with it.
4. Create and color the linear octree nodes (Section 3.3.2) and compute the distances to the input boundaries (Section 5.2.1).
5. Apply the tetrahedra pattern.(Section 3.3.3).

A 2D example of the algorithm is depicted in Table 5.2.4, where each step is illustrated by a figure. The model is formed by two surfaces in contact (drawn in orange and blue), and it contains only two mesh size points (represented by a cross in the figures) in the input data.

 1- Process input data and create the octree root. The crosses are the mesh size points. The arrow beside each mesh size point represents its associated size. The octree root is the black square, and the ${\textstyle Bbox_{m}^{+}}$ of the model is drawn with dotted lines. ${\textstyle S_{max}}$ in this example coincides with ${\textstyle S_{2}}$, as it is the largest size coming from the generalized mesh size points, and it is smaller than ${\textstyle s_{box}}$. 2- Refine the octree accomplishing the size refinement criteria. In this example the size transition factor is equal to one (the transition corresponds to the constrained two to one condition) to make the figure clearer. 3- Classify the input boundary entities into the octree. It is important to note that, until now, the input boundaries only have been considered to build the bounding box of the model, but they have not been implied in the octree refinement process. 4- Create and color the linear octree nodes and compute the distances to the input boundaries. The nodes of the figure are painted with the corresponding color: orange or blue if they are inside the corresponding surface, black if they are onto the input boundaries, and white if the are outside the domain. 5- Apply the tetrahedra pattern. It can be appreciated in this figure that the final tetrahedra (triangles in this 2D case) are not body-fitted, as they do not preserve the original shape of the domain.

Steps followed by the embedded meshing algorithm applied to a 2D example with two surfaces and two mesh size points.

## 5.3 Body-fitted mesher

The concept of body-fitted mesh is applied to a mesh which contours represent precisely the shape of the contours of the domain. Considering the presented algorithm is a volume mesher generating a body-fitted mesh implies, from the geometrical point of view, that the final tetrahedra mesh must take into account lower level entities: point, line and surface ones.

The adaptation of the final mesh to the surface entities defining the different volumes of the domain plays a key role to ensure the preservation of the topology of the model and the accuracy of the mesh near the contours. This process is detailed in Section 5.3.6.

In the literature, the concept of geometrical features is used to refer the point and line entities relevant for the shape definition of the model, which must be preserved in the final tetrahedra mesh. Point entities to be preserved are the ones where the surface normal has multiple discontinuities and are referred as corners. Line entities to be preserved are the ones shared by two surfaces forming a sharp angle in it and are referred as ridges.

An example of two meshes of the same model (one preserving and the other one not preserving the geometrical features) is shown in Figure 41.

 (a) Geometrical model. (b) Non body-fitted mesh. (c) Body-fitted mesh. Figure 41: Example of body-fitted and non body-fitted mesh of a geometrical model of a cone. The non body-fitted mesh does not preserve the apex of the cone (corner), nor the line entities (ridges) defining its base.

The preservation of the geometrical features is one of the weak points of the octree-based meshers. As these family of meshers are based in a regular grid (the octree) rather than in the shape of the domain, when the shape of the domain has specific distorted regions in comparison with the octree cells, some strategy has to be followed. One can think that refining the octree in those regions should solve the problem, but this is not true, as some configurations can lead to an infinite refinement process. This is because the geometrical problems governed by angles are reproduced exactly in all the refinement levels.

In this work, not only the geometrical features are taken into account to be preserved in the final mesh, but also forced point and forced line entities coming from the input data, needed to preserve the topology of the model or to provide with specific attached data to the final mesh entities. The generalized concept of forced edges and forced nodes (detailed in Sections 5.3.1 and 5.3.2) is used to include all the line and point entities to be preserved, independently on its purpose: preservation of topology or preservation of geometrical features.

### 5.3.1 Forced edges

In case the input data for the mesher has forced line entities (from the input data), or a minimum angle for sharp edges is defined, the so called forced edges must be created. These are edges the final tetrahedra mesh will preserve. Forced edges plays a key role for preserving sharp edges and representative surface and line topology.

 (a) (b) (c)
Figure 42: (a) Contours of a volume with some of its sharp edges highlighted. (b) Constrained tetrahedra mesh of the volume highlighting the sharp edges corresponding to the ones in figure (a).(c) Partially constrained tetrahedra mesh of the volume highlighting the sharp edges corresponding to the ones in figure (a).

At this point it has to be noted that the mesher is partially constrained. In this context, the term partially constrained means that if there are some forced line entities to be preserved, a collection of edges from the final mesh should follow the path of those line entities. This is not as restrictive as totally constrained condition, which would force to have as much edges in the final mesh as the number of forced line entities, and bounded (each one of the edge) by the same nodes bounding the forced line entities.

This difference can be appreciated in the Figure 42, where the definition of a domain is shown with some of its sharp edges highlighted (Figure 42(a)), and two different final meshes of the domain being partially or totally constrained are shown. In this example forced edges come only from the sharp edges of the domain. As it can be seen in Figure 42(b), generating a totally constrained mesh yields a mesh with the same sharp edges present in the input data. In Figure 42(c) it can be appreciated that a different number of edges are generated to represent the initial forced edges: generating a partially constrained mesh, a collection of sharp edges follows the path of the sharp edges present in the input data, so as the shape of the domain to be meshed is well captured.

The creation of the forced edges is based on three steps: identification of the base line entities, creation of the polyline entities and linear mesh generation from the polyline entities:

1. Identification of the base line entities. The base line entities are line entities coming from three different sources:
• Forced line entities: each forced line entity from the input data is directly a base line entity.
• Line entities shared by more than two surface entities from the input boundaries. A T-junction is an example of this case, where a line entity is shared by three surface entities.
• Line entities shared by two surface entities from input boundaries which form an angle (along the line entity) smaller than the maximum angle for sharp edges (in case it is defined in the input data). The angle between two surface entities is not well defined, as it can only be evaluated point by point. In the present work, the angle between two surface entities along their common line entity is defined as the mean of the angles formed by the normal vectors of each surface entity in some sampled points of the line entity. The number of points used depends on the curvature of the surface and line entities.

Note that the base line entities is a collection of line entities in the space, and they can be related or not to the input boundaries of the domain.

In case of non watertight input boundaries, a previous collapse of nodes and edges may be done to the input boundaries in order to be able to capture the sharp edges, as two surface entities forming a sharp edge may not be in contact topologically. This aspect is treated later on.

2. Creation of the polyline entities. In this step, the collection of base line entities is clustered in different groups called polyline entities in such a way that two base line entities are in the same polyline if:
• they share a point entity,
• and the shared point entity does not belong to any other base line entity,
• and the angle formed by the tangent vectors of each line entity in the shared point is smaller than the maximum angle for sharp edges.

The result of this clustering operation leads to a collection of polyline entities having each one of them one or more base line entities. Two point entities can be identified as extremes of each polyline entity.

The polyline entities corresponding to the base line entities shown in Figure 43(a) are depicted in Figure 43(b). In this example base line entities ${\textstyle l_{2}}$ and ${\textstyle l_{3}}$ are not part of the same polyline entity because they form a small angle in the common point (smaller than an hypothetical maximum angle for sharp edges).

 (a) (b) (c)
Figure 43: Example illustrating the entities involved in the creation of the forced edges. (a) A collection of base line entities ${\displaystyle l_{1}}$, ${\displaystyle l_{2}}$, ${\displaystyle l_{3}}$, ${\displaystyle l_{4}}$, ${\displaystyle l_{5}}$, ${\displaystyle l_{6}}$, ${\displaystyle l_{7}}$, ${\displaystyle l_{8}}$ and ${\displaystyle l_{9}}$. (b) The polyline entities ${\displaystyle pl_{A}}$, ${\displaystyle pl_{B}}$, ${\displaystyle pl_{C}}$, ${\displaystyle pl_{D}}$ and ${\displaystyle pl_{E}}$ created from the base line entities shown in Figure (a) considering the angle between ${\displaystyle l_{2}}$ and ${\displaystyle l_{3}}$ smaller than the maximum angle for sharp edges (from in input data). (c) A possible distribution of forced edges (shown in dotted line) created from the polylines present in Figure (b).
3. Generate a mesh from the polyline entities. In this step a linear mesh is generated from each polyline entity taking into account three criteria for the linear element creation:
• If a polyline entity is closed (both extreme nodes are the same one), it must have at least three elements. This condition avoid the creation of zero-volume elements in the final mesh.
• If the base line entities forming the polyline entity has some desired size assigned (from the input data), a desired size ${\textstyle s_{m}}$ can be computed for each linear element. This ${\textstyle s_{m}}$ is the mean of the desired sizes of the base line entities in the positions of the nodes of the element. Each element of the mesh must have a length equal or lower than its ${\textstyle s_{m}}$.
• The chordal error of the elements of each polyline meshes is limited by the most restrictive value between the chordal error required by the simulation (in case it is defined in the input data) and a given relative chordal error ${\textstyle Ec_{l}}$. The relative chordal error is defined as the chordal error of the element divided by its length. The parameter ${\textstyle Ec_{l}}$ tries to ensure that the forced edges describe correctly the shape of the forced line entities and the input boundaries. The value of this parameter is discussed in Section 6.7.

The discretization of the polyline entities can be done using any mesh generation method. This linear mesh will be used only auxiliary (it will not be part of the final mesh) and its quality is not relevant more than providing a sort of sizes distribution in the octree. In this work a simple recursive splitting method is used to generate the mesh of each polyline entity. The method creates one linear element using the two extreme nodes of the polyline entity, and subdivide it recursively until all the elements accomplish with the three criteria defined above. It has to be considered that each new edge node created when an element is subdivided is mapped onto the polyline entity in order to capture well the shape of it.

The forced edges are directly the elements of the meshes generated from the polyline entities.

Note that the two first steps refer to line entities in general, without specifying if they are CAD or mesh entities. The forced edges created in the third step are always mesh entities.

For notation purposes, the extreme nodes of a forced edge are called edge nodes. Each forced edge has two (and only two) edge nodes, and each edge node can belong to more than one forced edge.

It will be seen later on that the octree cells near a forced edge should have a similar size to it not to produce too distorted tetrahedra. To reach this goal, the forced edge condition has been defined:

Forced Edge ConditionForced Edge Condition

Condition 1: Being cell A the octree leaf containing one edge node of a forced edge, and cell B the octree leaf containing the other edge node of the forced edge, the degree of neighborhood between A and B must be lower or equal than two.

Considering a given configuration of the octree, all the forced edges must accomplish the Forced Edge Condition. If a forced edge violates it, it must be subdivided in other forced edges until all of them accomplish it. The split of the edge is done by its middle point. It has to be taken into account that the subdivision of a forced edge implies the creation of new edge nodes. These ones are mapped onto the polyline entity where the forced edge comes from to yield a better approximation of the shape.

#### 5.3.1.1 Pathological configurations of forced edges

There are some pathological situations for forced edges that can occur due to given configurations of input boundaries:

• Very thin surface entities. The presence of very thin surface entities in the boundary may be a problem due to the distortion in the computation of their normal vector. In case of triangles, the normal vector is computed using the vector product of two of their adjacent sides. If one of the them is very small (almost zero), the normal vector calculated becomes erroneous: it can point to any direction. For this reason, when comparing normal vectors of adjacent triangles in order to set sharp line entities, some fake ones may be detected.
• Very close forced line entities. The presence of very near forced line entities without intersection may lead to very distorted elements in the final mesh. Furthermore, some of the octree refinement criteria (Section 5.3.3) may lead to infinite level of refinement in this situations.

A possible strategy to solve the first situation is to collapse the line entities with a length lower than a given tolerance, or directly exclude these small ones as forced line entities. Collapsing the small line entities may not solve the problem, because they may belong to a longer polyline entity (also fake) which won't be collapsed. The option of not considering the small line entities also may fail, because there may be some cases where a relevant polyline entity is formed by very small line entities.

For the second situation, a possible solution may be to collapse the line entities which are close enough one from each other. However, these geometrical operations are not trivial in some 3D configurations.

In order to detect the sharp line entities, only the ones belonging to two, and only two, surface entities are considered. Cases where more than two surface entities share a line may be important for the topological definition of the domain. In these cases, the corresponding line entity should be set as forced line entity in the input data, not in the sharp edges detection process. The edges surrounding a gap (in case of non-watertight boundaries) are not considered to be preserved.

These automatic strategies may not solve all the pathological situations that can occur in the input boundaries. Hence, it may be needed to preprocess them in order to specify the desired forced line entities following given criteria to accomplish the simulation requirements. This is the case of the example shown in Figure 44. The input boundaries of this example are depicted in Figure 44(a). As it can be seen, there are very thin triangles. The normal vector of these triangles is not well computed numerically. When comparing the normal vector of two adjacent triangles, some edges which are not sharp are detected as so. It is the case of the small edges depicted in Figure 44(b).

 (a) (b) (c)
Figure 44: (a) Input boundaries (triangles) of a part of an example model of a mechanical piece. (b) Zoom view of Figure (a), where the shape is smooth enough not to present sharp edges. (c) Same view of (b), showing only the sharp line entities (black lines) automatically detected considering the normal vectors of the triangles of the boundaries (dotted lines surround them in order to highlight their position). As there are very thin triangles, some normal is not well computed, so some of the detected sharp line entities are fake.

In this example, some of the fake sharp edges (forced line entities) are connected creating the corresponding polyline entity. These polyline entities are not so small, so collapsing the nodes of the small forced line entities will not eliminate these fake ridges: the polyline entity would remain, but with less forced line entities. A possible solution here is to consider only the polyline entities longer than a given tolerance (as the polyline entities useful to define the shape of the domain are much more longer). However, this strategy cannot be applied in general. Some models may not have a so clear separation between the lengths of relevant and not relevant sharp entities.

In cases where the 3D model is not very complex, another option could be to select manually a priori the line entities to be forced ones.

In conclusion: Depending on the quality of the input boundaries of the domain, more information may be needed in order to identify clearly the forced line entities to be preserved by the mesher. In these cases, the geometrical definition of the domain is not enough to set the appropriate forced edges.

### 5.3.2 Forced nodes

The forced nodes are octree nodes (linked to an octree position) with a prescribed position in space: the forced position. The forced position of a forced node is not coincident with the octree position it is linked to. There are three kinds of forced nodes:

• Forced isolated nodes. These ones are created directly from the forced point entities given by the input data. Every forced point entity becomes a forced isolated node whose forced position is the one of the forced point entity and the related octree position is the closest one to it.
• Forced nodes in edge. These forced nodes correspond to the edge nodes (Section 5.3.1). Every edge node become a forced node in edge which forced position is the one of the edge node and the related octree position is the closest one to it.
• Forced interface nodes. These correspond to octree nodes which octree position is very near to the input boundaries. Taking into account an octree node, lets define ${\textstyle \mathbf {cp} }$ as its closest coordinate onto the input boundaries. Considering ${\textstyle D}$ the distance from an existing octree node to its ${\textstyle \mathbf {cp} }$, and ${\textstyle s_{c}}$ the size of the octree leaf the octree node is into, an octree node is a forced interface node if

 ${\displaystyle D<\alpha _{ip}\cdot s_{c}}$
(5.3)

where ${\textstyle \alpha _{ip}}$ is a real positive value. The value of ${\textstyle \alpha _{ip}}$ must be lower than ${\textstyle 0.25}$, otherwise it could lead to inverted tetrahedra (with negative jacobian) when applying the tetrahedra pattern (Section 3.3.3). The tuning of this value is explained in Section 6.7.

The forced position of a forced interface node is the ${\textstyle \mathbf {cp} }$ of the octree node it is generated from. It is important to note that the definition of the forced interface nodes depends on the size of the octree cells, so an octree node can be considered as a forced interface node or not depending on the size of the cells surrounding it. It can be seen that refining an octree can lead to the creation or deletion of some forced interface node.

It will be seen in Section 5.3.3 that the forced nodes are involved in an octree refinement criterion (RC 6). To anticipate the general idea, each forced node is linked to an octree position, and each octree position cannot have more than one forced node associated. So if two forced nodes are very close one from each other, the octree cell containing them should have a size similar to the distance between them. This could lead to an excessive level of refinement of the octree, specially in cases where there are non-watertight definitions of the boundaries.

In extreme cases where the forced nodes can be almost exactly in the same position, the refinement criteria could lead to an infinite level of refinement. To solve this problem, the forced nodes which are closer than a given tolerance (a portion of the mesh desired size) are collapsed. It has to be considered that if two forced nodes are collapsed and they belong to a forced edge, it must be collapsed as well.

### 5.3.3 Octree refinement criteria for body-fitted meshes

All the size refinement criteria (the ones applied for embedded meshes, detailed in Section 5.2.2) are also applicable to body-fitted meshes, as they are based on the desired mesh size.

In this section the specific refinement criteria (RC) needed for body-fitted meshes are defined. Some of them may not be applied in the meshing process depending on the input parameters.

These refinement criteria are a key point for the body-fitted mesher, as they ensure the geometrical features of the domain and its topology will be preserved by the tetrahedra generated from the octree.

#### 5.3.3.1 Forced nodes refinement criterion

The following criterion refers to the forced nodes. As it has been explained in Section 5.3.2, each forced node is associated to an octree cell position.

RC 6: If an octree cell has a forced node inside, and the octree position associated to that forced node is occupied by another forced node, the cell must be subdivided.

In order to ensure that this refinement criterion does not lead to an infinite level of refinement, the distance between forced nodes must be greater than a minimum value, corresponding to the minimum cell size allowed in the octree. The fulfillment of this condition involves a treatment of the forced nodes before the refinement criterion: if two forced nodes are closer than the minimum cell size allowed in the octree they are collapsed into one. If they are part of a forced edge, the forced edge is also collapsed.

#### 5.3.3.2 Tetrahedra distortion refinement criterion

For the refinement criteria defined hereafter, the tetrahedra created from the octree following the tetrahedra pattern are considered. For this reason, these criteria only take sense if the RC 4(balance) and RC 6(forced nodes) are accomplished. From now on the concept of the tetrahedra created from an octree cell is used to refer the tetrahedra result from applying to the cell the patterns defined in Section 3.3.3.

The quality of the tetrahedra generated from cells without forced nodes is very high (and it can be evaluated a priori), but the presence of forced nodes implies that their positions in space are not the octree positions, so the tetrahedra generated may be distorted.

 Figure 45: Tetrahedron with the local node numeration following the right hand rule. Vectors involved in the definition of an inverted tetrahedron are depicted.

The concept of inverted tetrahedron must be defined at this point in order to introduce the following refinement criterion. Considering a tetrahedron with its nodes ${\textstyle n_{1}}$, ${\textstyle n_{2}}$, ${\textstyle n_{3}}$ and ${\textstyle n_{4}}$ sorted in a given way (as the tetrahedron depicted in Figure 45), it is inverted if the scalar product of ${\textstyle (\mathbf {e_{1}} }$x${\textstyle \mathbf {e_{2}} )}$ by ${\textstyle \mathbf {e_{3}} }$ is negative. ${\textstyle \mathbf {e_{1}} }$, ${\textstyle \mathbf {e_{2}} }$, and ${\textstyle \mathbf {e_{3}} }$ are vectors aligned with the directions ${\textstyle n_{1}n_{2}}$, ${\textstyle n_{1}n_{3}}$ and ${\textstyle n_{1}n_{4}}$ respectively, oriented from ${\textstyle n_{1}}$ towards the other nodes. In the case where the scalar product is null, the tetrahedron has zero volume and receives the name of sliver.

RC 7: If a tetrahedron created from an octree cell is inverted or a sliver, the cell must be subdivided.

#### 5.3.3.3 Topological refinement criterion

The following refinement criteria are the basis to ensure the preservation of the topology for the mesher. Some auxiliary definitions are needed here. If we consider an edge of a tetrahedra generated from a cell, a limit distance of the edge ${\displaystyle d_{edge}}$ is defined as

 ${\displaystyle d_{edge}=\alpha _{edge}\cdot l_{edge}}$
(5.4)

where ${\textstyle l_{edge}}$ is the length of the edge and ${\textstyle \alpha _{edge}}$ is a real value between zero and one (the value taken for this parameter is detailed in Section 6.7).

The so called intersection points of an edge are the intersection points between the edge and the input boundaries. As an intersection operation between an edge and surface entities, the following situations must be considered (see also Section 3.4):

• If the intersection between an edge and the the surface entity defining the input boundaries is co-planar (${\textstyle P}$ intersection type), any intersection point is considered.
• For the evaluation of the following refinement criterion, a ${\textstyle MIP}$ is created only if all the intersected surface entities involved interface the same volumes.
• As the extremes of an edge are already colored, if both belong to different volumes and no intersection point has been detected, the corresponding ${\textstyle GIP}$ must be considered as the intersection point.

If we define the portion of a given volume enclosed in an octree cell delimited by the input boundaries and the cell faces, the interface cells can have several portions of volumes. Let us name face limit surfaces to the boundaries of these portions of volumes laying onto the the cell faces. An example of the portions of volumes of a cell is shown in Figure 46.

 Figure 46: Example of portions of volumes enclosed in an octree cell. Volume ${\displaystyle A}$ lies onto one face of the cell creating one face limit surface, volume ${\displaystyle B}$ does not lay onto any face of the cell, volume ${\displaystyle C}$ lies onto three faces of the cell creating three connected face limit surfaces, and volume ${\displaystyle D}$ lies onto two faces of the cell creating two unconnected face limit surfaces.

RC 8: This criterion is split in three levels:

• 8(a) If an edge of a tetrahedra created from an octree cell fulfill some of these conditions:
• both extreme nodes of the edge are forced nodes and there are more than one intersection point,
• one extreme node of the edge is a forced node and the distance between it and an intersection point is lower than ${\textstyle d_{edge}}$,
• there are two intersection points and the distance between them is lower than ${\textstyle d_{edge}}$, or
• there are more than two intersection points,

then the cell must be subdivided.

• 8(b) If a portion of volume enclosed in an octree cell does not intersect any edge from the tetrahedra created from the cell, and it lies onto more than one face of a cell creating unconnected face limit surfaces, then the cell must be subdivided.
• 8(c) If a portion of volume enclosed in an octree cell does not lay onto any face of the cell and does not intersect any edge from the tetrahedra created from the cell, then the cell must be subdivided.

The implementation of these refinement criteria is detailed in Section 6.5.

Examples of different configurations of an edge fulfilling the refinement criterion 8(a) are shown in Figure 47.

 (a) (b) (c) (d) Figure 47: Different configurations of the edge ${\displaystyle {\overline {AB}}}$ that force the refinement of the octree in order to accomplish the refinement criterion 8(a). Dotted line represents the input boundaries and crosses are the intersection points between the edge and the input boundaries. (a) Both ${\displaystyle A}$ and ${\displaystyle B}$ nodes are forced nodes and there are more than one intersection point. (b) ${\displaystyle A}$ is a forced node (${\displaystyle B}$ could be forced node or not), and there is an intersection point closer than ${\displaystyle d_{edge}}$ to it. (c) there are two intersection points and the distance between them is lower than ${\displaystyle d_{edge}}$. (d) There are more than two intersection points.

An example of a pathological configuration where the refinement criterion 8(b) is needed is the portion of volume ${\textstyle D}$ of Figure 46. On the other hand, the portion of volume ${\textstyle B}$ of the same figure evidences the need for the fulfillment of the criterion 8(c) in order to preserve the topology of the input data.

### 5.3.4 Meshing algorithm

Hereafter, the steps of the meshing algorithm for body-fitted meshes are presented:

1. Process input data and create the octree root.
• Creation of forced edges and forced nodes (in case they exist).
• Creation of generalized mesh size points (in case they exist).
• Compute ${\textstyle S_{max}}$.
• Creation of the octree root. The creation of the octree root is done in the same way as for embedded mesher (Section 5.2.3).
2. Refine the octree according to the size refinement criteria (the ones used for the embedded mesher, defined in Section 5.2.2).
3. Classify the input boundary entities into the octree. Here is where the input boundaries (surface entities from the input data) are analyzed in order to determine the interface cells. Before this step the input boundaries only have been taken into account for building the bounding box of the model.
4. Color the octree nodes (Section 4) and detect the forced interface nodes.
5. Refine the octree to fulfill the specific criteria for body-fitted meshes (defined in Section 5.3.3):
• RC 6: Forced nodes. This implies enter the forced nodes into the octree, which is assign an octree position to each one of them.
• RC 7: Tetrahedra distortion. (See page current).
• RC 8: Topology. (See page current).

This step implies the coloring of the appearing octree nodes. RC 4 (balance) has to be also taken into account during this refinement process, as it is mandatory for the tetrahedra generation following the given patterns.

From now on, the octree structure is frozen in the sense that its cells will not be refined any more.

6. Apply tetrahedra pattern (Section 3.3.3). In this step, for all the interface and inner cells (not the outer ones) the tetrahedra are created following the tetrahedra patterns.
7. Preserve geometrical features (Section 5.3.5). After this step all the forced edges are edges of the tetrahedra mesh, and the octree nodes linked to a forced node have been placed in the corresponding position.
8. Apply surface fitting algorithm (Section 5.3.6). Some nodes are placed onto the surface entities defining the contours of the volumes, and some tetrahedra may be split.
9. Color tetrahedra and create skin triangles of volumes (Section 5.3.7). After this step all the tetrahedra can be assigned to a specific volume or the outer part of the domain. Now the tetrahedra and the nodes of the outer part of the domain can be deleted. The skin triangles of the volumes is formed by the faces (triangles) between tetrahedra of different color, or interfacing a tetrahedron with the outer part of the domain.
10. Make-up and smoothing. This step can be considered out of the mesh generation algorithm itself, as it consists in the improvement of the quality of the elements generated, and it can be applied to any tetrahedra mesh. However, it is highly necessary considering the quality of the tetrahedra in the contours of the domain (Section 5.3.8).

Note that the refinement process involved in the fifth step is an iterative process. Every time an octree cell is subdivided, several aspects have to be considered in the new configuration of the octree:

• New octree nodes are created, so they have to be colored.
• Some of the new octree nodes may be forced interface nodes.
• Some previously existing forced interface nodes may become regular nodes (not forced ones).
• Some forced edge may be subdivided in order to accomplish the maximum Forced Edge Condition(Section 5.3.1).
• New forced nodes (interface or edge ones) may be created.

All this aspects may force the subdivision of other cells due to the other refinement criteria, so an iterative process is required in order to achieve an octree configuration where all the refinement criteria are accomplished. However, in our experience, few iterations are enough for satisfying all of them. The implementation of the algorithm is detailed in Section 6.5.

Another important characteristic of the algorithm is that the first six steps are based on the octree, while from the seventh step onwards, the operations are applied to the unstructured tetrahedra mesh.

A graphical example of the steps of the meshing process is shown in Table 5.3.4 using a 2D model to make the figures more understandable. In the first figure, the offset between ${\textstyle Bbox_{m}^{+}}$ and the octree root should be equal to ${\textstyle S_{max}}$, but it has been put smaller to make the following figures clearer.

It also has to be considered that some of the parts of the algorithm are intrinsically 3D (such as the process to preserve forced edges), so they cannot be illustrated with a 2D model. This example is formed by two surfaces and it has no mesh size information in the input data. It is important to note the presented algorithm is able to generate a body-fitted mesh with the only information of the input boundaries.

 1- Process input data and create the octree root. The red small squares represent forced points. The octree root is the black square, and the ${\textstyle Bbox_{m}^{+}}$ of the model is represented with dotted lines. In this example ${\textstyle S_{max}}$ coincides with the minimum side of the model bounding box (${\textstyle s_{box}}$), as there are no mesh size points. 2- Refine the octree accomplishing the size refinement criteria. As there are no mesh size point in the model, the octree is refined uniformly with ${\textstyle S_{max}}$. 3- Classify the input boundary entities into the octree. It is important to note that until now, the input boundaries only have been considered to build the bounding box of the model, but they have not been implied in the octree refinement process. 4- Color the octree nodes and set forced interface points. The nodes of the figure are painted with the corresponding color: orange or blue if they are inside the corresponding surface, black if they are close enough to the contour entities to be forced interface nodes. Nodes outside the domain are white. 5- Refine the octree according with the body-fitted refinement criteria. It can be appreciated that the octree refinement process leads to the creation of new octree nodes. Some of them become forced isolated nodes (red ones) or forced interface nodes (black ones). The arrows indicates the forced position of each forced node. 6- Apply the tetrahedra pattern. The tetrahedra are generated from the octree nodes. From this step on, the operations are performed to the corresponding tetrahedra mesh rather than the octree one. 7- Preserve geometric features. At this point forced nodes are moved into their forced positions. As a 2D example, only the moving of forced nodes can be appreciated in the figure (the preservation of edges has no sense in 2D). 8- Surface fitting process. This is the last step of the process which ensures the final mesh could represent the original topology of the input boundaries. 9- Tetrahedra coloring and identification of skin mesh. Elements owning to the orange or blue surface are painted accordingly in the figure. After this process the outer elements can be deleted. 10- Make-up and smoothing. In this step some mesh editing operations are done in order to improve the mesh quality.

Steps followed by the body-fitted meshing algorithm applied to a 2D example with two surfaces with no mesh size assigned.

### 5.3.5 Preserve geometric features

The preservation of geometrical features (corners and ridges) lies on the preservation of forced nodes (Section 5.3.2) and forced edges (Section 5.3.1).

The operations described in this section are performed after the creation of tetrahedra from the octree leaves following the tetrahedra pattern (Section 3.3.3), so at this point there are forced nodes, forced edges and a tetrahedra mesh got directly from the octree cells. This tetrahedra mesh covers the domain to be meshed, but does not preserve the geometrical features and does not fit the surface entities representing the interfaces between volumes (there are nodes outside the domain and tetrahedra with edges crossing volume interfaces). The aim of this process is that the tetrahedra mesh has nodes in the corresponding positions of the forced nodes, and edges corresponding to the forced edges.

Each forced node has an octree node associated to it. The process of preserving the forced nodes is reduced to move the corresponding octree nodes to the position of its forced node.

 (a) (b) Figure 48: Process of splitting a forced edge by inserting a node onto its base line. (a) Forced edge ${\displaystyle {\overline {AB}}}$ to be split by the node ${\displaystyle N}$; the dotted line is the base line of the forced edge. (b) Forced edges ${\displaystyle {\overline {AN}}}$ and ${\displaystyle {\overline {NB}}}$ result from splitting the forced edge ${\displaystyle {\overline {AB}}}$ by the node ${\displaystyle N}$; dotted lines are the corresponding base lines of the forced edges.

Concerning forced edges, the goal is to force the tetrahedra mesh to have edges coincident with them. As explained in Section 5.3.1, the forced edges are obtained as a linear mesh from the polyline entities. The portion of polyline entity enclosed between the extreme nodes of a forced edge can be defined as the base line of the forced edge. Note that this is a pure notation, as the polyline is made of generic line entities (in mesh or geometrical format). For notation purposes, a forced edge will be defined by its extreme nodes: the forced edge ${\textstyle {\overline {AB}}}$ is the one which extreme nodes are node ${\textstyle A}$ and node ${\textstyle B}$. The same notation is used for a generic edge of a tetrahedron.

From now on the concept of fitting the tetrahedra mesh to the forced edges is used to define the process of having an edge of the tetrahedra mesh for each forced edge. To achieve this goal a splitting process of forced edges and tetrahedra is proposed. This process involves three basic operations:

• Split a forced edge by inserting a node onto its base line. This process is shown in Figure 48 and it creates two forced edges from an original one connecting its extreme nodes with the splitting node. The base line of the forced edge is split by the node as well.
• Split an edge of a tetrahedron. This process involves all the tetrahedra surrounding the edge. Each of them is split into two tetrahedra as shown in Figure 49(b).
• Split a face of a tetrahedron. This process involves a tetrahedron and the neighbor to the face. Each of these two tetrahedra generate three tetrahedra as shown in Figure 49(d).
 (a) (b) (c) (d)
Figure 49: Process of splitting tetrahedra by a node. (a) Edge ${\displaystyle {\overline {AB}}}$ to be split by the node ${\displaystyle N}$, with its surrounding tetrahedra around (4 in this example). (b) 8 tetrahedra result from splitting the edge ${\displaystyle {\overline {AB}}}$ by the node ${\displaystyle N}$. (c) Face ${\displaystyle ABC}$ to be split by the node ${\displaystyle N}$, with the two tetrahedra sharing the face. (d) 6 tetrahedra result from splitting the face ${\displaystyle ABC}$ by the node ${\displaystyle N}$.

Hereafter the process needed for fitting the tetrahedra mesh to the forced edges is detailed. For each forced edge ${\textstyle {\overline {AB}}}$ we check whether nodes ${\textstyle A}$ and ${\textstyle B}$ are nodes of a same tetrahedra or not. In case they are, the tetrahedra mesh has an edge coincident with the forced edge, so the objective is already achieved. In case there is no tetrahedra containing nodes ${\textstyle A}$ and ${\textstyle B}$, the following procedure must be considered (a graphical example of the steps followed to fit the tetrahedra mesh with a forced edge is shown in Table 5.3.5):

• Find the tetrahedron surrounding node ${\textstyle A}$ which opposite face with respect to ${\textstyle A}$ (from now on ${\textstyle face_{A}}$) is intersected by the base line of the forced edge. Note that the base line can be represented by geometrical or mesh entities (NURBS curve or linear mesh). In case of a linear mesh, the intersection operation between the base line and the faces of tetrahedra is much more robust and computationally cheaper.
• Find the intersection point ${\textstyle P}$ between the base line and ${\textstyle face_{A}}$. Defining ${\textstyle min_{side}}$ as the length of the shortest edge of ${\textstyle face_{A}}$, a node (${\textstyle N}$) is created taking care on the following casuistic (a graphical view of the entities involved is shown in Figure 50):
 Figure 50: Entities involved in the process of making the tetrahedra mesh to have an edge coincident with a forced edge. Red line is the forced edge ${\displaystyle {\overline {AB}}}$ and the dotted red line is its base line. The tetrahedron ${\displaystyle AF_{1}F_{3}F_{2}}$ is the tetrahedron surrounding node ${\displaystyle A}$ which opposite face with respect to ${\displaystyle A}$ (${\displaystyle face_{A}}$) is intersected by the base line. Face ${\displaystyle F_{1}F_{2}F_{3}}$ is the ${\displaystyle face_{A}}$ (drawn in yellow). ${\displaystyle P}$ is the intersection point between the base line and ${\displaystyle face_{A}}$. The closest node of ${\displaystyle face_{A}}$ to ${\displaystyle P}$ is ${\displaystyle F_{1}}$, and the closest edge from ${\displaystyle face_{A}}$ to ${\displaystyle P}$ is the edge ${\displaystyle F_{1}F_{2}}$.
• If the distance between ${\textstyle P}$ and its closest node of ${\textstyle face_{A}}$ is lower than (${\textstyle \alpha _{vertex}\cdot min_{side}}$), then this node is moved to the position of ${\textstyle P}$ and it is considered as the node ${\textstyle N}$. The parameter ${\textstyle \alpha _{vertex}}$ is a real value between zero and one; the value used for it in the presented implementation of the algorithm is discussed in Section 6.7. This is the case shown in figures d) and e) of Table 5.3.5.
• If the previous case is not accomplished and the distance between ${\textstyle P}$ and its closest edge of ${\textstyle face_{A}}$ is lower than (${\textstyle \alpha _{side}\cdot min_{side}}$), node ${\textstyle N}$ is created in the position of ${\textstyle P}$, and the closest edge is split by the node ${\textstyle N}$. The parameter ${\textstyle \alpha _{side}}$ is a real value between zero and one; the value used for it in the presented algorithm is discussed in section 6.7. This is the case shown in figures f) and g) of Table 5.3.5.
• If the previous cases are not accomplished, node ${\textstyle N}$ is created in the position of ${\textstyle P}$ and ${\textstyle face_{A}}$ is split by the node ${\textstyle N}$. This is the case shown in figures b) and c) of Table 5.3.5.
• Split the forced edge ${\textstyle {\overline {AB}}}$ by the node ${\textstyle N}$. Now the forced edge ${\textstyle {\overline {AN}}}$ is an edge of the tetrahedra mesh.
• If the forced edge ${\textstyle {\overline {NB}}}$ is not an edge of the tetrahedra mesh, repeat the process with it.

Note that the nodes created in this process are not octree nodes, as they are not related with any octree position. Furthermore, these nodes do not need any coloring algorithm, as they lay on a interface entity.

 (a) Initial configuration of the forced edge ${\textstyle {\overline {AB}}}$ (red line) with its surrounding tetrahedra (in blue), and its base line (dotted red curved line). (b) Find the tetrahedra (${\textstyle AF_{1}F_{3}F_{2}}$) owning node ${\textstyle A}$ which opposite face to it (${\textstyle F_{1}F_{2}F_{3}}$) intersects the base line of forced edge ${\textstyle {\overline {AB}}}$. ${\textstyle P_{1}}$ is the intersection point. (c) As ${\textstyle P_{1}}$ is not close enough to ${\textstyle F_{1}}$, ${\textstyle F_{2}}$ or ${\textstyle F_{3}}$, and neither to the edges of face ${\textstyle F_{1}F_{2}F_{3}}$, create the node ${\textstyle N_{1}}$ in the position of ${\textstyle P_{1}}$ and split the face ${\textstyle F_{1}F_{2}F_{3}}$ and the forced edge ${\textstyle {\overline {AB}}}$ by the node ${\textstyle N_{1}}$. Now the forced edge ${\textstyle {\overline {AN_{1}}}}$ is already an edge of the tetrahedra mesh. (d) Proceed the treatment of forced edge ${\textstyle {\overline {N_{1}B}}}$. Find the tetrahedra (${\textstyle N_{1}F_{2}F_{1}F_{4}}$) owning node ${\textstyle N_{1}}$ which opposite face to it (${\textstyle F_{1}F_{2}F_{4}}$) intersects the base line of forced edge. ${\textstyle P_{2}}$ is the intersection point. (e) As ${\textstyle P_{2}}$ is close enough to ${\textstyle F_{4}}$ (closer than ${\textstyle \alpha _{vertex}}$ times the minimum edge of ${\textstyle F_{1}F_{2}F_{4}}$), move ${\textstyle F_{4}}$ to the position of ${\textstyle P_{2}}$ and split forced edge ${\textstyle {\overline {N_{1}B}}}$ by node ${\textstyle F_{4}}$. Now the forced edge ${\textstyle {\overline {N_{1}F_{4}}}}$ is already an edge of the mesh. (f) Proceed the treatment of forced edge ${\textstyle {\overline {F_{4}B}}}$. Find the tetrahedra (${\textstyle F_{4}F_{5}F_{7}F_{6}}$) owning node ${\textstyle F_{4}}$ which opposite face to it (${\textstyle F_{5}F_{6}F_{7}}$) intersects the base line of forced edge. ${\textstyle P_{3}}$ is the intersection point. (g) As ${\textstyle P_{3}}$ is close enough to edge ${\textstyle F_{5}F_{7}}$ (closer than ${\textstyle \alpha _{side}}$ times the minimum edge of ${\textstyle F_{5}F_{6}F_{7}}$), creation of node ${\textstyle N_{3}}$ in the position of ${\textstyle P_{3}}$ and split of edge ${\textstyle F_{5}F_{7}}$ and the forced edge ${\textstyle {\overline {F_{4}B}}}$ by the node ${\textstyle N_{3}}$. Now the forced edges ${\textstyle {\overline {F_{4}N_{3}}}}$ and ${\textstyle {\overline {N_{3}B}}}$ are edges of the mesh, so the process is finished. (h) Final configuration of the tetrahedra with the new forced edges created: ${\textstyle {\overline {AN_{1}}}}$, ${\textstyle {\overline {N_{1}F_{4}}}}$, ${\textstyle {\overline {F_{4}n_{3}}}}$ and ${\textstyle {\overline {N_{3}B}}}$.

Example of the process to fit a tetrahedral mesh with the forced edge ${\displaystyle {\overline {AB}}}$.

### 5.3.6 Surface fitting

This section describes the process of fitting the tetrahedra mesh into the volumes of the domain to be meshed, representing their interface surfaces accurately (from now on surface fitting process). The methodology presented is applied to the tetrahedra mesh (come from the tetrahedra pattern defined in Section 3.3.3), in which the process of preserving geometrical features (defined in Section 5.3.5) has been carried out. Some of the nodes of the mesh are forced nodes (fixed in a position in space) and some of its edges are coincident with the forced edges. The process defined from now on preserves the forced nodes as well as the edges corresponding to the forced edges.

The surface fitting process is based on the isosurface stuffing method published in [37] (from now on isostuffing method). However, it has some differences, as its objectives and restrictions are different:

• The isostuffing method is applied to an isolated volume. In this work all the volumes of the domain are meshed together, and they can be in contact one to each other.
• The isostuffing method do not preserve sharp edges nor corners. The present algorithm must respect the forced nodes and forced edges in order to preserve the topology of the model and its geometrical features.
• The isostuffing method does not preserve the topology of the domain, as the size of the cells used is chosen by the user a priori; if this size is much more greater than a thin part of the domain, this part will not be represented by the final mesh. The surface fitting process presented in this work ensures the topology of the mesh will be the same as the one of the input geometry. This is due to the topological refinement criteria applied to the octree the mesh is generated from.

An example of tetrahedra mesh generated using the isostuffing method is shown in Figure 51. It can be appreciated that the sharp edges relevant for the definition of the shape of the model are not preserved by the mesher.

 (a) (b) Figure 51: Example of a 3D model of a mechanical piece (a) and a tetrahedra mesh of it generated using the isostuffing method (b).

The surface fitting process is based on the edges of the tetrahedra mesh which are not forced edges: the iso-edges. It consists in the following steps:

1. Compute iso-edges intersections. Compute the intersection of the iso-edges with the input boundaries (intersection points). As described in Section 3.4, the following situations must be considered:
• If the intersection between an edge and the the surface entity is co-planar (${\textstyle P}$ intersection type), no intersection point is considered.
• Intersection points closer than ${\textstyle tol_{c}}$ are collapsed into a ${\textstyle MIP}$.
• As the extremes of an edge are already colored, if both belong to different volumes and no intersection point has been detected, the corresponding ${\textstyle GIP}$ must be considered as the intersection point.

Some of the edges can intersect more than one time the input boundaries, but not more than two due to the topology criterion (a). In these cases with two intersections only one of them is taken into account (no matter which one). From now on we only take care of the edges intersecting the input boundaries; the ones that do not intersect them are taken out from the iso-edges.

2. Move nodes. For each node of the iso-edges, we consider its closest intersection point (between the intersection points of the edges containing the node). It can be assigned to each of these nodes (${\textstyle nod}$) the closest intersection point (${\textstyle ip}$) and the edge the ${\textstyle ip}$ is into (${\textstyle e}$). Then, if ${\textstyle nod}$ is not a forced node, move it to the position of ${\textstyle ip}$ if
3.  ${\displaystyle d<\alpha _{iso}\cdot l_{edge}}$

where ${\textstyle d}$ the distance between ${\textstyle nod}$ and ${\textstyle ip}$, and ${\textstyle l_{edge}}$ the length of the edge ${\textstyle e}$. The parameter ${\textstyle \alpha _{iso}}$ is a real value between zero and one. The value used in the algorithm is explained in Section 6.7. This process of moving the node is only performed if the resulting configuration does not generate poor quality tetrahedra.

Moving a node implies to recompute the possible intersection points of all its connected edges, as one of their extremes is moved. This leads to an iterative process where new edges can be added to the iso-edges, and some existing one can be taken out (the ones with no intersection point). A 2D example of a situation where the movement of a node creates a new intersection point is shown in Figure 52.

 (a) (b) Figure 52: 2D example of creation of a new intersection point when moving a node. Black line represents the input boundaries, and part of the triangle mesh is shown in blue. (a) Initial configuration: node ${\displaystyle A}$ is close enough to intersection point ${\displaystyle ip_{1}}$, so it is moved. (b) Final configuration after moving ${\displaystyle A}$: the new intersection point ${\displaystyle ip_{3}}$ has appeared.
4. Split edges. Split the remaining iso-edges (the ones with intersection point after the node moving process) by their intersection point reconnecting the surrounding tetrahedra, as shown in Figures 49(a) and (b).
5. Repeat steps 1, 2 and 3 with a new set of iso-edges: the ones of the mesh having at least one if its nodes onto the input boundaries.
6. Repeat steps 1, 2 and 3 with a new set of iso-edges: the ones of the mesh fulfilling the following conditions (both of them):
• At least one of the nodes of the edge is onto the input boundaries.
• At least one of the nodes of the edge is connected (by an edge) to a forced node in edge.

 (a) Initial configuration of the mesh before surface fitting process. (b) Detection of the intersection points (red dots). Nodes ${\textstyle M_{1}}$ and ${\textstyle M_{2}}$ are close enough to its closest intersection points (red circumferences) to be moved. (c) Mesh configuration after moving points ${\textstyle M_{1}}$ and ${\textstyle M_{2}}$. Some intersection point have disappeared and a new one have appeared. Light gray lines represent the mesh in the previous step. (d) Splitting process by the intersection points. Dotted lines are the new edges created. There is an edge with two intersection points (${\textstyle A}$ and ${\textstyle B}$). Node ${\textstyle B}$ has been used for the splitting (${\textstyle A}$ could be used as well). Edges in gray represent outer edges with no intersections. (e) Second iteration. As there are no nodes to be moved, edges with intersections are split. There is an edge with three intersection points. Point ${\textstyle C}$ is chosen (arbitrarily) for the splitting. (f) Final configuration of the mesh.

The accomplishment of the topology criterion (a) allows the presence of edges intersecting twice the domain in the mesh coming from the tetrahedra pattern. For this reason the moving and splitting process has to be applied the second time to the iso-edges (step 4). The involved iso-edges this second time have at least one of their nodes onto the input boundaries (as they have already been processed the first time).

A 2D example applying two times the process of moving nodes and splitting edges is shown in Table 1.

 (a) (b) Figure 53: Example of mesh generated applying (a) or not applying (b) the process of moving nodes and splitting edges the third time. It can be appreciated in Figure (a) that the topology of the mesh near a forced edge does not represent well the topology of the domain.

For pathological situations, applying a third time the process of moving nodes and splitting edges is needed (step 5). These pathological situations happen when some of the new edges created in the previous steps has an intersection point relevant for the topology of the mesh. These cases occur near sharp edges, so the iso-edges involved are only the ones connected to a forced node in edge. An example is shown in Figure 53, where the difference between applying or not the moving and splitting process the third time can be appreciated. If there were not forced edges it will not be needed to repeat a third time the moving and splitting operations, because a set of tetrahedra would already represent topologically well each volume of the domain.

One can think about repeating the process of moving nodes and splitting edges as many times as needed until the mesh does not present any intersection point. This is not feasible as it could lead to almost infinite loops (in some configurations with curved input boundaries). Also, each time it is applied, the new elements will have worse quality.

As explained before, the process of moving nodes may create new intersection points. Also the splitting process can create more of them, as it implies the creation of new tetrahedra with the corresponding new edges. The new edges created may have more than two intersections with the input boundaries, but these intersection points are not relevant to preserve the topology of the volumes, so it is not needed to take them into account.

A 2D example of the described surface fitting process is depicted in Table 1. In this example, the process of moving and splitting could be repeated more times, as there are more intersection points (in edge ${\textstyle {\overline {FC}}}$ of the mesh shown in (f)). However, the topology of the input boundaries is already well represented, so no more iterations are needed.

### 5.3.7 Tetrahedra coloring

Once the surface fitting and preserving features operations have been performed, the mesh is in the following situation: it is ensured that the tetrahedra represent the topology of the input boundaries, and all the nodes of the mesh are colored (it is known if they are inside a volume or onto interfaces between volumes). Now the operation of tetrahedra coloring must be done: this is, assign a volume to each one of the tetrahedra. In this work the concept of assigning a color to an entity is used to determine inside which volume the entity is, so color and volume are used indistinctly.

The tetrahedra to be colored can present two basic situations:

• At least one node of the tetrahedra is inside a volume (inner node).
• All the nodes of the tetrahedra are onto interfaces between volumes (surface or line entities). From now on this kind of tetrahedra will be referred as interface tetrahedra.

Tetrahedra presenting the first situation are obvious to be colored. Because of the topological refinement criteria and the surface fitting properties, it is ensured that there is no edge of the the tetrahedra with nodes inside different volumes. This implies the tetrahedra belongs to the same volume of the inner node.

The second situation (interface tetrahedra) is much more complicated to solve. It has to be noted that each node of an interface tetrahedra is interfacing a number of volumes; they are the so called possible volumes of a node. Considering the forced edges lying on an interface between volumes, the concept of possible volumes of a forced edge can be defined analogously as for the nodes.

The possible volumes of each tetrahedron are clearly identified: they are the volumes which are interfaced by all the nodes of the tetrahedron. There are some configurations under which the color of an interface tetrahedron can be determined. This is the case where the tetrahedron has only one possible volume. The rest of interface tetrahedra are considered as undetermined ones; these are the cases that must be solved.

 Figure 54: 2D example of a triangle mesh of two surfaces: A (blue) and B (orange). Elements containing node ${\displaystyle M}$ are directly colored as ${\displaystyle A}$, as ${\displaystyle M}$ is a inner node to A. The other triangles are interface elements. Elements ${\displaystyle HGJ}$ and ${\displaystyle JPK}$ are undetermined, as they have two possible colors: ${\displaystyle B}$ and ${\displaystyle 0}$ (exterior).

A 2D example of these different types of elements (in terms of coloring process) is shown in Figure 54. Here the mesh for two surfaces (${\textstyle A}$ in blue and ${\textstyle B}$ in orange) is shown. All elements containing node ${\textstyle M}$ belong to surface ${\textstyle A}$, as ${\textstyle M}$ is an inner node to ${\textstyle A}$. The other elements are interface ones. The possible colors of each one of the interface nodes are listed in Table 2 (note that the outer part of the domain is considered as color zero):

 Node Possible colors E ${\displaystyle A}$, ${\textstyle 0}$ F ${\displaystyle A}$, ${\textstyle B}$ , ${\textstyle 0}$ G ${\displaystyle B}$, ${\textstyle 0}$ H ${\displaystyle B}$, ${\textstyle 0}$ J ${\displaystyle B}$, ${\textstyle 0}$ K ${\displaystyle B}$, ${\textstyle 0}$ L ${\displaystyle A}$, ${\textstyle B}$ N ${\displaystyle A}$, ${\textstyle 0}$ P ${\displaystyle A}$, ${\textstyle B}$ , ${\textstyle 0}$

Considering the interface elements, their possible volumes can be easily gotten analyzing the common possible volumes of each one of their nodes. They are depicted in Table 3.

 Element Possible colors ELF ${\displaystyle A}$ FLG ${\displaystyle B}$ HGJ ${\displaystyle B}$, ${\textstyle 0}$ JPK ${\displaystyle B}$, ${\textstyle 0}$ JLP ${\displaystyle B}$

It can be seen that the elements ${\textstyle HGJ}$ and ${\textstyle JPK}$ are undetermined in this example, as they can belong (topologically) to surface ${\textstyle B}$ or the outer part of the domain indistinctly.

As each undetermined tetrahedron has more than one possible volume, there are several possible configurations of colors for them. A configuration is understood as the assignment of a color to each one of the undetermined tetrahedra. It has to be noted that the triangles of the skin meshes of each volume are the faces of the tetrahedra shared by tetrahedra of different color, so each configuration of colors lead to a different tetrahedra skin mesh. Taking care about the topology of the input boundaries, the final configuration of colors has to accomplish the following conditions to be considered as valid:

1. Between the tetrahedra surrounding one node, it must be at least one tetrahedron for each possible volume of the node. This implies, for example, that the tetrahedra surrounding a node in interface cannot be all of them of the same color.
2. Between the tetrahedra surrounding one forced edge, it must be at least one tetrahedron for each possible volumes of the forced edge.
3. Considering the input boundaries in the near region from a given node in interface, if the input boundaries are manifold in this near region, the skin of the tetrahedra surrounding the node must be manifold. This condition can only be applied for watertight volumes and implies, for example, that if a node is in a manifold region of the input boundaries, there must be a closed loop of triangles around the node in the skin of volume's tetrahedra.
4. The mesh of a volume must present the same topology as the volume. This means, for example, that if the volume has ${\textstyle n}$ holes, its mesh must have also ${\textstyle n}$ holes.

The accomplishment of these conditions can reduce to one the possible colors of some tetrahedra. In this cases, these ones would not be undetermined. However, other tetrahedra may remain undetermined.

A good 2D example of undetermined elements considering only two different colors (interior and exterior) is shown in the (f) figure of Table 1 (see page current). The triangles ${\textstyle AGB}$, ${\textstyle LBG}$, ${\textstyle LGM}$, ${\textstyle CEF}$, ${\textstyle FED}$, ${\textstyle DHI}$ and ${\textstyle KLJ}$ have all their nodes in the boundaries. ${\textstyle AGB}$, ${\textstyle DHI}$ and ${\textstyle KLJ}$ are not undetermined elements, as they can only be colored as inside the surface taking into account the first condition defined above: if they were colored as outside, nodes ${\textstyle A}$, ${\textstyle H}$ and ${\textstyle K}$ respectively won't have any element inside the surface, and they are nodes interfacing the surface and the outer part. Then, in this example, the elements ${\textstyle LBG}$, ${\textstyle LGM}$, ${\textstyle CEF}$ and ${\textstyle FED}$ are the undetermined ones. These elements could be colored as interior or exterior each of them, and both situations may lead to topologically correct meshes. This example shows that the problem of coloring the undetermined tetrahedra has more than one solution, and it is not obvious (actually, in some cases it is impossible) to determine whether a solution is better or worse than another.

 (a) (b) Figure 55: 2D examples of pathological configurations for element coloring. Both cases are taken from the mesh shown in (f) figure of Table 1. (a) The element ${\displaystyle LBG}$ is a triangle of the surface, but its center (white dot) is outside it. (b) The element ${\displaystyle CEF}$ could be colored as inside or outside of the surface; the oriented normal vectors in nodes ${\displaystyle C}$, ${\displaystyle E}$ and ${\displaystyle F}$ does not point to the center of the element. Black arrows are the oriented normal vectors towards the surface, and the white dot is the center of the triangle.

A possible strategy for the tetrahedra coloring should be to assign the tetrahedra the color of their center using the node coloring technique explained in Chapter 4, but this can lead to wrong decisions as there are several pathological configurations where the center of an element of a volume is not inside it. This is the case of triangle ${\textstyle LBG}$ in (f) figure of Table 1. A zoom of this triangle is shown in Figure 55(a): ${\textstyle LBG}$ is topologically a triangle of the surface ${\textstyle A}$, but its center is outside it.

Another possible strategy for undetermined tetrahedra should be take care about the normal vectors of the interface surface entities in the nodes of the tetrahedra. This normal vector can be oriented in the sense that it can be determined which volume it is pointing to (between the volumes interfaced by the surface entity). The notation oriented normal to a volume V is used from now on to indicate the normal vector of a surface entity pointing towards the volume V.

The computation of the normal of a surface of the volume is not obvious taking into account that the contours of the domain can be non-watertight. We are leaving this aspect unsolved by now, and let us suppose the oriented normal vectors can be created. One could think about situations which the oriented normal vectors can determine whether the color of a tetrahedron is one or another depending on where they are pointing to. For example: the cases where the four normal vectors of an interface tetrahedron (one for each node) oriented to a given volume are pointing at the center of the tetrahedron could be considered as determined: the color of the tetrahedron should be the volume one. (In this context, point to the center of a tetrahedron means that the center is in the semi-space defined by the oriented normal vector).

Unfortunately, not all the tetrahedra accomplish with this condition, so there may still be some undetermined tetrahedron. This is the case presented in Figure 55(b): it is a 2D case where it can be seen that the oriented normal vectors to the surface in the nodes of the element ${\textstyle CEF}$ are not leaving the center of the element in the same semi-space.

Furthermore, the same topological situations can be represented by different volume boundaries which present a drastically different configuration of normal vectors in the nodes of the elements. This is the case shown in Figure 56, where different interfaces between surface ${\textstyle A}$ and ${\textstyle B}$ present the same element with the normal vectors in two of its nodes identical, and the third normal different for each shape of the boundaries. This 2D example demonstrates that the process of tetrahedra coloring cannot be based on the normal vectors of the interface entities in the nodes of the elements. It has to be taken into account that in 3D the possible configurations are much more complicated than in 2D.

 (a) (b) (c) (d) Figure 56: 2D example. The black curved line represent the interface between surfaces ${\displaystyle A}$ and ${\displaystyle B}$. (a), (b), (c) and (d) represent different shapes of the interface. In all the cases an undetermined triangle is depicted and the vectors are the normal vectors pointing to surface ${\displaystyle B}$. As it can be seen, in all the cases the element should be colored as inside surface ${\displaystyle A}$. The normal vectors in two of the nodes are identical for all the cases, and the third (the upper one) is different in all of them.

The reason why the two strategies defined above fail is they are based on geometrical conditions, and the problem to solve is more topological than geometrical. The color of an undetermined tetrahedron is independent on the percentage of the element inside or outside a volume: it is determined for topological conditions (the ones defined previously). To take care about the topology, a strategy is proposed to color the undetermined tetrahedra starting from the tetrahedra with known color. It is based in the following proposition:

Proposition 1: Considering an undetermined element ${\textstyle T}$ and an element ${\textstyle R}$ of color ${\textstyle A}$ which is neighbor of ${\textstyle T}$ by face ${\textstyle F}$. Being ${\textstyle N}$ the node of ${\textstyle T}$ opposite to face ${\textstyle F}$, if there is a point ${\textstyle P}$ inside volume ${\textstyle A}$ laying on face ${\textstyle F}$, and a continuous curve in space from ${\textstyle P}$ to ${\textstyle N}$ completely inside element ${\textstyle T}$ with no intersection with the boundaries of A, then the color of ${\textstyle T}$ is ${\textstyle A}$.

A graphical view of this proposition in a 2D case is shown in Figure 57. In this example the element ${\textstyle NML}$ is undetermined, as all its nodes are on the interface between surface ${\textstyle A}$ and the outer part of the domain. All the other elements are fully determined (they belong to surface ${\textstyle A}$), as they contain node ${\textstyle R}$, which is inner to ${\textstyle A}$. Point ${\textstyle P}$ is inner to ${\textstyle A}$, and it lays on the face shared by elements ${\textstyle NML}$ and ${\textstyle MRL}$. As it exists a continuous curve (shown in dotted red line) from ${\textstyle N}$ to ${\textstyle P}$, totally inner to element ${\textstyle NML}$ and with no intersections with the boundaries of ${\textstyle A}$, element ${\textstyle NML}$ can be considered inner to surface ${\textstyle A}$. It has to be noted that this process imply the coloring process of point ${\textstyle P}$, which is done following the algorithm described in Chapter 4.

 Figure 57: Graphical interpretation of Proposition 1 in a 2D example. The black line represents the contour of surface ${\displaystyle A}$, and its triangle mesh is represented by blue lines.

Proposition 1 is used to color the undetermined tetrahedra (the implementation of the algorithm is detailed in Section 6.5). However, it cannot be applied in some cases as the ones where an undetermined tetrahedron has no neighbor by face with a determined color, or the cases where the points in the common face are onto a surface entity (they are not inner to a volume), or the coloring process of these points fails. Another case where the proposition cannot be applied is when the volumes involved are non-watertight. For this reason, it has to be planned that after this coloring process some undetermined tetrahedra can still remain. From now on the methodology to color the remained undetermined tetrahedra is explained.

As it can be seen, the coloring of one tetrahedron can affect the coloring of its neighbors, as the skin of the tetrahedra changes and can affect the manifold condition of each node. However, each tetrahedron affects only to its neighbors, so there can be defined different clusters of undetermined tetrahedra which are independent one from each other; these clusters are made by the undetermined tetrahedra connected at least by one node. As an example, the clusters of undetermined tetrahedra of the mesh in figure (f) of Table 1 are identified in Figure 58. As these clusters are independent one from each other, the methodology presented can be applied on a cluster by cluster manner.

 Figure 58: Two clusters of undetermined elements (blue and orange) of the mesh shown in (f) figure of Table 1.

Considering all the undetermined tetrahedra of a cluster, each of them has its possible colors. All the possible color configurations can be obtained, and each one of these configurations implies a different skin mesh of each volume's tetrahedra. Checking the topological conditions that must be accomplished (detailed previously in this section) for the final mesh, some of the configurations are not valid, and some of them are valid. For the tetrahedra coloring purposes, any one of the valid configurations is set as the solution.

As an example, lets see the case of the orange cluster of undetermined elements shown in Figure 58. It is made of two elements: ${\textstyle LBG}$ and ${\textstyle LGM}$. Each of them has two possible colors: inside or outside the surface. All the different configurations are shown in Figure 59. The configuration shown in Figure 59(a) is not valid because the topology of the final mesh is not the same as the one of the surface it is representing (the mesh has two unconnected sets of elements). The configurations of Figure 59(b) and (c) are not valid because they have non-manifold nodes in regions where the contour of the surface is manifold (node ${\textstyle G}$ in case (b) and node ${\textstyle L}$ in case (c)). The only valid configuration is the one shown in Figure 59(d).

 (a) (b) (c) (d)
Figure 59: Zoom of the orange cluster of elements in Figure 58. All the possible configurations taking into account the different colors (inside or outside the surface) of elements ${\displaystyle LBG}$ and ${\displaystyle LGM}$.(a) Both elements are outside the surface. (b) ${\displaystyle LBG}$ is outside and ${\displaystyle LGM}$ inside. (c) ${\displaystyle LBG}$ is inside and ${\displaystyle LGM}$ outside. (d) Both elements are inside.

### 5.3.8 Make-up and smoothing

The make-up and smoothing operations are used to improve the quality of the mesh once it has been generated. Make-up operations are the ones changing the topology of the mesh, and smoothing process only implies movement of the nodes, maintaining the original connectivity of the elements.

The acceptable quality of the mesh elements is a relative parameter, as the mesher aims to be applied to different kinds of numerical simulations, and they are applicable in different ranges of qualities for the elements  [5]. A clear example of this is the boundary layer meshes, which elements should have an aspect ratio higher than 10000 in some cases; such distorted elements are not valid for standard FEM analysis. However, almost all the methods require non-inverted elements. This means that the Jacobian of the transformation of the element (from the parametric to the 3D space) should be positive. This lead to require, at least, a strictly positive Jacobian in the elements of the final mesh. A graphical interpretation of an inverted element is provided in Section 5.3.3.

In this section there are references to tetrahedra, triangles and edges. Edges are edges of the tetrahedra elements, and triangles are considered as the faces of the tetrahedra interfacing tetrahedra of different color. Thus, each volume tetrahedra mesh is surrounded by a triangle mesh.

Considering the mesh generated may have forced nodes and forced edges, there are some restrictions to be applied to the make-up and smoothing operations:

• Forced isolated nodes (Section 5.3.2) cannot be moved nor deleted.
• Forced nodes in edge which are extreme of a polyline entity (Section 5.3.1) cannot be moved nor deleted.
• Forced nodes in edge can only be moved along the polyline entity they belong to.
• Forced interface nodes can only be moved lying onto the surface entity they belong to.
• If both extreme nodes of a tetrahedron edge are forced nodes in interface and the edge is not a triangle edge, the edge cannot be collapsed.
• If both extreme nodes of a tetrahedron edge are forced nodes in edge and the edge is not a forced edge, the edge cannot be collapsed.
• If a tetrahedron edge has one forced node in edge and one forced node in interface, it cannot be collapsed.

From now on, the operations involving the collapse of an edge, or the movement of a node are only applied to the entities not violating these restrictions. It has to be taken into account that the collapse operation involves the deletion of one of the two nodes involved.

The elements obtained from the tetrahedra pattern defined in Section 3.3.3 have very good quality if the nodes are in the octree positions (see Section 3.3.2), so the make-up and smoothing operations are only needed for the tetrahedra containing a forced node, the ones coming from interface cells and the ones resulting from surface fitting or preserving features operations.

The parameters ${\textstyle \alpha _{vertex}}$ and ${\textstyle \alpha _{side}}$ (Section 5.3.5) try to guarantee a minimum quality in the tetrahedra resulting from the preserving features process, and ${\textstyle \alpha _{iso}}$ tries to do the same for surface fitting operations. However, some configurations of the mesh and the input boundaries may lead to low quality elements. This situation often causes the presence of small edges in the mesh. An example of these small edges is shown in Figure 60(a). To improve the quality of the elements surrounding these small edges, an edge collapsing step is performed (make-up operation). Taking into account that the octree cell containing each node gives an idea of the mesh size required for the mesh in that region (because of the user desired sizes or as a result of a topological refinement process), the edges smaller than a given portion of the size of the cell they are inside are collapsed. An edge is collapsed if it does not violate the restrictions defined at the beginning of this section and

 ${\displaystyle l_{edge}<\alpha _{c}\cdot s_{c}}$
(5.5)

where ${\textstyle s_{c}}$ is the size of the smaller octree cell between the ones where the extreme nodes of the edge are inside, ${\textstyle l_{edge}}$ the length of the edge and ${\textstyle \alpha _{c}}$ is a real value greater than zero. Its value is discussed in Section 6.7.

 (a) (b) Figure 60: Example of a mesh (a) before and (b) after the make-up and smoothing process.

The smoothing operation applied consists in a Laplacian-like smoothing that displaces a node into a position such that its surrounding elements improve their quality. The smoothing process is applied on a node by node basis, and follows these steps:

1. Compute the quality of all the tetrahedra surrounding the node, and identify the worst of them (${\textstyle q_{0}}$).
2. Move the node to a candidate position.
3. Compute the quality of the tetrahedra surrounding the node in the new configuration and identify the worst of them (${\textstyle q_{1}}$).
4. If ${\textstyle q_{1}}$ is better than ${\textstyle q_{0}}$, the node is left in the candidate position, otherwise, the node is returned to its original one.

This process implies the definition of a quality measure of the elements, as well as a procedure to obtain the candidate position. As a quality measure, the minimum dihedral angle of the tetrahedra has been chosen, considering that a tetrahedra is worse than another one if it has a lower minimum dihedral angle. The election of the candidate position is explained later on.

It has to be noted that the mesh can have three kinds of elements: tetrahedra, triangles and linear elements (corresponding to the forced edges). All the nodes have tetrahedra around them, but only some of them have triangle or linear (1D) elements. Taking into account the restrictions defined in the beginning of this section, the degrees of freedom for moving a node are restricted by its nature: if it is a forced node in edge, it can only move along the polyline entity corresponding to the forced edge, and if it is a forced interface node it must lay on the surface entity it belongs to. This aspect governs the candidate position for each node:

• If the node is a forced node in edge, the mean position ${\textstyle m}$ of the opposite nodes of the forced edges it belongs to is calculated. Considering the polyline entity those forced edges come from, the candidate position is the closest point to ${\textstyle m}$ laying on the polyline entity.
• If the node is a forced interface node (it has triangles around it), all the nodes of the triangles around it (except the node itself) are considered. Their mean position ${\textstyle m}$ is calculated, and the candidate position is the closest point to ${\textstyle m}$ laying on the surface entity the triangles belong to.
• If the node has no triangles nor forced edges, the candidate position is the mean position of the nodes of the surrounding tetrahedra of the node (except the node itself).

As the movement of a node affects the quality of all its surrounding elements, the smoothing operation is thought as an iterative process where several loops over all the nodes are performed in order to improve the quality of the mesh.

Apart from edge collapsing and nodes smoothing, the edge flipping operation (a make-up operation) is performed in the mesh [55]. The cases 2 to 3, 3 to 2, 4 to 4 and 5 to 6 are implemented. These operations imply the removal of a face (the first case) creating an edge, or the removal of an edge (the other cases) creating some extra faces. Apart from the pathological configurations described in [55] which determine the situations where the edge flipping cannot be made, some topological configurations must be considered in order not to invalidate the topology of the mesh generated. In particular, the following restrictions are considered:

• If an edge is a forced edge or the edge of a triangle, it cannot be deleted.
• If a tetrahedra face is a triangle of the mesh, it cannot be deleted.

The collapse, edge flipping and smoothing operations are applied in a sequential iterative manner to account for the updated configurations of the mesh each time. A maximum of four iterations is set in the present work. The result of applying them to the mesh shown in Figure 60(a) is depicted in Figure 60(b). In Chapter 7 several examples are shown indicating the mesh quality before and after the make-up and smoothing operations.

### 5.3.9 Mesh quality

This section focuses in the analysis of the quality of the final mesh generated by the new mesher.

For the embedded mesher, the quality of the mesh is totally guaranteed, as all the tetrahedra come directly from the predefined patterns. However, for the body-fitted mesher, it has to be considered that a body-fitted mesh is forced to respect the boundaries of the domain. If a part of the domain is bounded forming a very small dihedral angle, the tetrahedra representing it will have the same dihedral angle. This limits the scope of minimum quality guarantee to the inner parts of the domain, or the boundaries which are relatively smooth.

The tetrahedra of the inner parts of the domain come directly from the predefined patterns, so they have a very good quality. Their minimum dihedral angle is ${\textstyle 45}$ degrees (Section 3.3.3). For the tetrahedra near the boundaries there are several aspects that make impossible to ensure a minimum quality for the elements in the final mesh:

• Forced nodes. As the forced nodes are not placed in the octree positions, the tetrahedra pattern generated may yield tetrahedra that are distorted with respect to the ones coming from regular octree nodes. The tetrahedra distortion refinement criterion (RC 7) avoids inverted tetrahedra, which ensures a minimum dihedral angle equal to zero. A solution could be to set a higher minimum dihedral angle, thus ensuring a better quality for the elements. However, considering that the boundaries of the domain may have small dihedral angles, this could lead to infinite refinement in certain regions, because it may be impossible to represent a given shape from a tetrahedron coming from an octree cell.
• Preserving forced edges operations. These operations involve displacement of nodes and split of tetrahedra edges and faces. Parameters ${\textstyle \alpha _{vertex}}$ and ${\textstyle \alpha _{side}}$ govern the behavior of these operations. Depending on the position of the intersection point of the forced edge and the given tetrahedron face, these parameters establish whether a node must be displaced or a splitting operation must be performed. Theoretically, the tuning of these parameters could bound the quality of the final configuration of tetrahedra after these operations, if the minimum element quality of the initial configuration is known. This is not the case, because of the presence of forced nodes (as explained in the previous paragraph), so a minimum element quality cannot be ensured after the process of preserving the forced edges.
• Surface fitting process. The isostuffing method [37] ensures a minimum element quality by tuning two parameters that govern the decision of moving a node or splitting an edge in the surface fitting process. This is possible because the method does not preserve the geometrical features. The mesh where these operations are applied comes directly from a predefined tetrahedra pattern, which quality is known a priori. In our case, the quality of the mesh before the surface fitting process is not bounded, as explained in the previous two points. Hence, the tuning of the parameter ${\textstyle \alpha _{iso}}$ tries to get well shaped tetrahedra, but cannot guarantee any minimum quality for the final mesh.

For all these aspects, although the make-up and smoothing operations described in Section 5.3.8 reach an acceptable quality for the meshes generated, a given minimum element quality cannot be guaranteed theoretically.

### 5.3.10 Extension for surface meshing

As explained in Section 1.2.4, a secondary objective of the present work is to be able to apply the octree mesher to the meshing of surfaces not belonging to any volume. In this section the concept of inner surface is used to refer a surface entity not belonging topologically to the boundaries of a volume. It can be inside or outside a volume. From the mesher point of view this difference takes no sense as the outer part of the model is considered as volume zero.

Taking into account the new mesher can provide with the surface mesh corresponding to the boundaries of a volume, it could seem that it can be applied directly to an inner surface. Actually, after the surface fitting process there are tetrahedra at both sides of the inner surfaces that conform to them. However, the faces of the tetrahedra which form the triangle mesh of the inner surface cannot be detected in the same way as the surface entities boundary of a volume (regural surface entities).

Extracting the surface mesh of a regular surface entity is automatic after the process of tetrahedra coloring (described in Section 5.3.7): the triangles of the surface entities are the faces of the tetrahedra interfacing two tetrahedra of different color. In the case of inner surfaces, tetrahedra at both sides of the triangles have the same color, so the same strategy cannot be applied. Some modifications have to be performed in order to detect the triangle mesh of an inner surface.

At the time of extracting the surface meshes of the inner surfaces all the tetrahedra have been generated, the processes for preserving geometrical features and surface fitting have been carried out, and the tetrahedra have been colored. At this point, there are already tetrahedra faces corresponding to the triangles of the inner surface mesh, but they have to be detected. This situation is illustrated in Figure 61(a) using a 2D example (an inner line and the surrounding triangles is used instead of an inner surface and the surrounding tetrahedra).

 (a) (b) (c) (d) Figure 61: 2D example of the steps for detecting the line elements of an inner line. (a) Initial configuration: elements at both sizes of the inner line. Thick black line is the inner line and the mesh is in blue. Black dots are forced nodes in interface or in edge. (b) All the candidate 1D linear elements in dotted black lines. (c) Candidate 1D linear elements (in black) accomplishing the topological properties set to definitive. (d) Final mesh (in black) for the inner line once the remaining invalid candidate linear elements have been discarded by Proposition 2.

A first consideration must be done: all the triangles in the mesh of a surface entity have their nodes laying on it. This means that the nodes of those triangles are forced nodes in interface or in edge. The tetrahedra faces with this characteristics (all their nodes being forced nodes in interface or in edge) are called candidate triangles. In Figure 61(b) the candidate linear elements (in 2D there are candidate 1D linear elements instead of candidate triangles) of the configuration shown in Figure 61(a) are depicted in dotted lines.

Considering all the candidate triangles of an inner surface, some topological properties analogous to the ones made in the tetrahedra coloring strategy can be done. If a node is inside a surface entity, the surface mesh around it must be manifold (let us call it a manifold node). This implies that there should be a closed loop of triangles surrounding it in the mesh of the surface entity. If there are a closed loop of candidate triangles around a manifold node, and they are the only candidate triangles containing the node, they can be set as valid (triangles belonging to the inner surface final mesh). Also the requirement for every forced node in interface or in edge to belong at least to one triangle of the final mesh can be used: if one of these nodes belongs to one candidate triangle only, this one must be valid. For the case shown in Figure 61(a), the candidate elements set as valid considering these topological properties are depicted in Figure 61(c).

However, this topological properties may not solve the problem of detecting all the triangles of an inner surface mesh, as there are nodes which are not needed to be manifold (the ones in edges), and there may be triangle configurations around manifold nodes which are not manifold. This leads to the need for another strategy to detect the right triangles among all the candidate ones. The strategy proposed is based on disregarding the wrong candidate triangles, rather than detecting the right ones, using Proposition 2. It is based on Proposition 1 (defined in section 5.3.7 for tetrahedra coloring), but with slight modifications:

Proposition 2: Consider a candidate triangle ${\textstyle T}$ of an inner surface ${\textstyle S}$ and the two tetrahedra around it ${\textstyle R_{1}}$ and ${\textstyle R_{2}}$. Being ${\textstyle P_{1}}$ and ${\textstyle P_{2}}$ the nodes of ${\textstyle R_{1}}$ and ${\textstyle R_{2}}$ not belonging to ${\textstyle T}$, if there is a point ${\textstyle P}$ inside triangle ${\textstyle T}$ and a continuous curve in space from ${\textstyle P_{1}}$ to ${\textstyle P_{2}}$ passing through ${\textstyle P}$ completely inside the union of ${\textstyle R_{1}}$ and ${\textstyle R_{2}}$ with no intersection with ${\textstyle S}$, the candidate triangle is invalid.

A graphical view of this proposition(using a 2D case) is shown in Figure 62.

 (a) (b) Figure 62: Graphical interpretation of Proposition 2 in a 2D example. Black line represents the inner line (${\displaystyle S}$). Triangle mesh is represented by blue lines and the candidate linear element to be checked (${\displaystyle T}$) is the dotted line. (a) The candidate linear element is invalid, because there is a curve (in red) not intersecting the inner line. (b) The candidate 1D linear element cannot be set as invalid because all the possible curves intersect the inner line.

In Figure 62(a) the candidate 1D linear element is set as invalid because there is a curve (drawn in red) accomplishing Proposition 2. Figure 62(b) shows that there could not be any curve accomplishing the proposition, so the candidate linear element cannot be set as invalid.

Once all the invalid candidate triangles are disregarded using Proposition 2, the rest of the candidate triangles are the ones corresponding to the final mesh of the inner surface (Figure 10(d)).

Note that this proposition can set as invalid some correct triangle in some cases, specially when the the surface entities are curved. Some improvements should be applied to the algorithm to be more robust in some pathological configurations.

The implementation of the algorithm defined for extracting the mesh of an inner surface is detailed in Section 6.5.4.

# 6 Implementation aspects

This chapter details all the implementation aspects relevant for the meshing algorithm described in the previous chapters. The specific implementation of the ray casting technique used for the nodes coloring is explained in Section 4.3, as it can be considered as an independent process.

It is important to note that implementation matters do not affect the result of an algorithm, but they affect its performance and efficiency. In this sense, the implementation of a meshing algorithm is crucial if it has to be applied to industrial problems with complex geometries, or to generate really large meshes. A bad implementation can lead to unaffordable problems in terms of memory and computational time.

The algorithm has been implemented as a static library. Although the implementation of the mesher is done taking into account the GiD pre and post-processing system [56,57,58] for getting the input data and the visualization of the meshes, its connection to GiD has been done as an external library, by a general interface specially designed to make it accessible for other programs. It is really important to provide access to the mesher as a library, as some numerical simulations require interaction with the mesher during the simulation itself (i.e. for optimization loops or in remeshing processes).

## 6.1 General aspects

The implementation of the new algorithm has been carried out paying special attention to saving as much memory as possible, and trying to improve its performance taking into account shared memory parallel processing (Section 6.6).

One important characteristic common to many meshing algorithms is that they need an extra memory to generate the mesh compared to the memory needed to store the mesh once it has been generated. This leads to the need of memory saving strategies in the implementation of the mesher:

• The main part of operations involving octree cells are applied only to inner and interface cells (not to outer ones), so they are the ones from which the tetrahedra are generated. These operations include cells refinement and creation of octree nodes. It has to be considered that in the first stages of the meshing algorithm it is not obvious to know if a cell is outer or not, as the coloring operations are not done yet. For this reason, in these stages a simple check is made in order to detect some of the outer cells in an easy way: if the cell is outside the bounding box of the model, it is taken as an outer cell. Checking whether a cell is inside or outside a bounding box is a really inexpensive process.
 (a) (b) Figure 63: 2D example of a thin model (the solid surface) with its bounding box (dotted line) and the quadtree (black lines). Cells out of the bounding box are marked with points. (a) Quadtree refined considering only the cells colliding the model bounding box. It can be appreciated that the quadtree is not balanced out of it. (b) Quadree refined considering all its cells.

Depending on the aspect ratio of the bounding box of the model, this strategy can save a lot of memory. This is the case of very thin models in some direction (like the 2D example depicted in Figure 63) where, actually, the main part of the octree root is out of the bounding box of the model. Not considering the refinement criteria of the octree out of the model bounding box (like, for instance, the balance criterion) can reduce drastically the number of cells of it, as it can be seen in Figure 63. It can be appreciated that the quadtree refined considering all the cells (b) has much more cells than considering only the cells colliding the model bounding box (a).

• Tetrahedra from the inner part of the volumes are generated at the very last step of the algorithm. Tetrahedra from the inner leaves have really good quality, they are not affected by the operations involved in the body-fitting process (preserve geometric features and surface fitting) and their color is totally determined (all its nodes have the same color). For this reason, these tetrahedra can be generated in the last step of the algorithm without affecting the resulting mesh. This part is implemented as follows:
• Once the octree is refined following all the refinement criteria, the interface cells and their neighbors (one level of neighborhood) are considered. Then, their tetrahedra are generated.
• The following operations are performed on these tetrahedra: preserve geometrical features, surface fitting process, tetrahedra coloring and make-up and smoothing operations. These operations require typically an extra topological information of the mesh (like the neighbor elements of a node). Once these processes are finished, this extra information is deleted in order to save memory.
• Then, as a last step, the tetrahedra generated from the inner cells are created. In this step no additional information than the mesh itself must be stored.

The flowchart of the body-fitted algorithm illustrating this process can be seen in Figure 69.

Considering the input data, the meshing algorithm is designed in a way that the input boundaries can be mesh or geometrical entities. The only operations required for them to be used in the algorithm are:

• Get the bounding box of the entity.
• Compute the intersection between an edge and the entity.
• Get the closest point of the entity to a given position of the space.

The algorithm has been implemented considering only mesh entities for the definition of the input boundaries. This simplify much the geometrical operations defined above.

## 6.2 Octree implementation

The octree is the key structure of the presented algorithm. As mentioned in previous sections, it is not only used for generating the tetrahedra, but also for searching purposes. Actually, octrees were created originally for this topic, so it is natural to take advantage from on its characteristics in this field.

The main operation the octree is designed for is to find the leaf containing a coordinate in space. For this purpose, beginning from the octree root, the given coordinate is analyzed to evaluate which of its children contains it (considering the uniform subdivision of a cell, determining in which octant of a cell the coordinate is represents an inexpensive process). If the child is not an octree leaf, the process is repeated (with itself instead of the octree root) until an octree leaf is reached. This process can be seen as going downstream by the octree structure, and often implies the use of recursive functions (the same function applied to a cell is applied to its child). The use of recursive functions is not desirable because it decreases the performance of the algorithm due to the function call overhead.

In order to achieve a good performance, a very efficient implementation of the octree has been carried out based on