Robust intersection of structured hexahedral meshes and degenerate triangle meshes with volume fraction applications

Two methods for calculating the volume and surface area of the intersection between a triangle mesh and a rectangular hexahedron are presented. The main result is an exact method that calculates the polyhedron of intersection and thereafter the volume and surface area of the fraction of the hexahedral cell inside the mesh. The second method is approximate, and estimates the intersection by a least squares plane. While most previous publications focus on non-degenerate triangle meshes, we here extend the methods to handle geometric degeneracies. In particular, we focus on large-scale triangle overlaps, or double surfaces. It is a geometric degeneracy that can be hard to solve with existing mesh repair algorithms. There could also be situations in which it is desirable to keep the original triangle mesh unmodified. Alternative methods that solve the problem without altering the mesh are therefore presented. This is a step towards a method that calculates the solid area and volume fractions of a degenerate triangle mesh including overlapping triangles, overlapping meshes, hanging nodes, and gaps. Such triangle meshes are common in industrial applications. The methods are validated against three industrial test cases. The validation shows that the exact method handles all addressed geometric degeneracies, including double surfaces, small self-intersections, and split hexahedra.


Introduction
The need for computing intersections between polyhedral objects is common in many applications, such as computer graphics, simulations, and robotics [12].Polyhedral objects are typically used as discrete representations of real objects, for example, when computational fluid dynamics (CFD) is used to solve multiphase flow problems.The surface of a solid object can then be represented by a triangle mesh and immersed in a background fluid grid consisting of rectangular hexahedra.The mesh needs to be connected to the background grid.In the connection step or later in the solution process, it can be necessary or useful to compute the intersection between the triangle mesh and each fluid cell, to find information about the cut cells.
Two well-known methods for connecting the mesh and the grid are cut cell methods [2,11,26] and immersed boundary (IB) methods [19,22].In cut cell methods, the fluid cells that are intersected by the mesh are refined to fit the mesh.In IB methods, body-fitted grids and cut cells are replaced by other approaches for connecting the mesh and the grid.The cut cell geometric information is clearly needed in the connection step of a cut cell method, but can be of interest also in an IB method.A typical application is in the calculation of fluxes.There is potential to improve the accuracy in the simulated fluxes if the fraction of each cell face inside the triangle mesh is known.This information can be extracted from the cut cell information.
Several methods for extracting the geometry of a cut cell are found in the literature.Some of these are reviewed in Section 2. However, few of the publications focus on degenerate triangle meshes.Degeneracies such as hanging nodes (T-vertices), gaps, cracks, overlapping meshes, or overlapping triangles are common in triangle meshes used for engineering applications.This aspect is important to consider when designing or using an algorithm that extracts the geometry of a cut cell.One option is to repair the mesh before it is given as input to the geometry extraction algorithm.Mesh repair is a broad field, partly surveyed in Section 2. Another, less common, option is to let the geometry extraction algorithm itself handle the degeneracies.This is done for example in [27].
In this article, we present a method similar to [1,2,11] for extracting the exact geometric information about a Cartesian cell cut by a triangle mesh.In contrast to earlier works, our method is geometrically robust in the sense that it can handle largescale triangle overlaps, what we will also call double surfaces, without the need for previously repairing the mesh.To some extent, our method also handles more general self-intersecting meshes.This is an improvement since large-scale triangle overlaps appear frequently in engineering applications and are hard to repair without undesirable side effects (see Section 2).Moreover, none of [1,2,11] go into details about the method and implementation, that is why there is still room for explanation of important concepts.We will consider some of these concepts here.
The method is in the following referred to as "the exact method."It is not completely exact since we use finite precision floating point arithmetics in our intersection routines.What we mean by exactness is that the exact polygons and polyhedrons of intersection are calculated given the intersection points.An approximate method where a least squares plane is fitted to the cell-triangle mesh intersection is also developed and compared to the exact method.The approximate method is intended for highly resolved hexahedral grids, for which it is reasonable to approximate the cell-triangle mesh intersection with a plane.It has some limitations and is a complement to the exact method.
What we essentially have solved is the purely geometric problem of intersecting a degenerate triangle mesh and a rectangular hexahedron.However, the proposed methods are specially designed to be used in an immersed boundary method [20].As such, we are interested in both the cell-triangle mesh intersection and the facetriangle mesh intersections.This has been taken into account in the development of the methods.

Related work
Different methods for extracting the geometry of a cut cell are found in the literature.In 1994, Quirk [23] introduces the Cartesian boundary method for complex two-dimensional geometries, where the cut cells are identified and classified into one of twelve types depending on how the solid and the cell boundary intersect.The intersection between the geometry and the cell is approximated by a line.This approximation is motivated by the fact that the existence of a cell with more than two intersection points indicates that the grid is not fine enough to resolve the solid properly.
Yang et al. [26] take this a bit further and describe a cut cell method in three dimensions where the part of the triangle mesh in the interior of the cell is approximated by a non-planar quadrilateral.They find the area of the quadrilateral and each cut face and use Gauss's divergence theorem to calculate the volume of the part of the cell located inside the triangle mesh.
In [1,2], Aftosmis et al. present an approach for completely resolving the geometry of a cut cell in three dimensions.They introduce triangle polygons, face polygons, and face segments to describe the polyhedron that represents a cut cell.To construct the triangle polygons, they use the Sutherland-Hodgman algorithm [25] for clipping each triangle against the cell boundary.They mention that face polygons are easily formed by connecting face segments with the edges of the cut cell.Cieslak et al. [11] also present a cut cell method that preserves the real geometry by finding the exact cut out polygons and polyhedron.They find the face polygons through connectivity criteria such as common faces or common triangles.
In [11,23,26], little or nothing is mentioned about geometrically degenerate triangle meshes, while Aftosmis et al. [2] adopt a volumetric approach to degenerate and overlapping triangle meshes.Volumetric approaches belong to one of two groups of mesh repair methods.The other group consists of surface-oriented methods.Surfaceoriented methods such as [7,15,18] operate directly on the degenerate geometry, while volumetric methods such as [9,17,21] create an intermediate volumetric representation of the geometry.The intermediate representation is used to create a new mesh without degeneracies.Most repair methods either have restrictions on the input mesh or drawbacks [4,10].Volumetric methods are typically robust, but destroy connectivity structures of the input geometry and could lead to loss of model features.Surface-oriented methods are often better at preserving details of the input mesh, but are typically not as robust as the volumetric methods.A survey and categorization of mesh repair approaches was done by Attene, Campen, and Kobbelt in 2013 [4].
A robust surface-oriented repair method that focuses particularly on selfintersecting meshes and double surfaces was proposed by Attene in 2014 [5].It computes the exact outer hull of the input triangle mesh by finding the self-intersections, subdividing the triangles along the self-intersections, and finally, removing excess triangles.Exact arithmetics is used to ensure correct intersections, but only when ordinary floating point operations are not sufficiently accurate.
Some approaches to a combination of mesh repair and geometric operations are found in the literature, where the method by Aftosmis et al. [2] is already mentioned.Another approach was recently proposed by Zhou et al. [27], who integrate mesh repair and set operations such as difference and union of an arbitrary number of input meshes.They assign a generalized winding number to different parts of space, defined by the intersections between the input meshes.The winding numbers are used to repair the meshes or perform set operations.
The repair method proposed in [5] could be used to remove the self-intersecting double surfaces, and [27] could be used to find the requested cut cell information if each hexahedral cell is first triangulated.In contrast to the three-dimensional set operations in [27], we present a method that operates on the two-dimensional faces of the intersection between the triangle mesh and the hexahedron.Thus, our method is adapted for finding both cell-mesh and face-mesh intersections.Only minor changes to the original method are required to handle double surfaces and some other types of self-intersections.Due to the small computational overhead of handling the degenerate cases, we claim that our method is a good alternative to the previous methods.

Problem formulation
Let T = {T i } n T i=1 be a triangle mesh, where T i ⊂ R 3 , i ∈ {1, . . ., n T }, are triangles.T is represented by a set of indexed vertices {v i } n v i=1 , and each triangle T i is defined by three vertices.Note that n v ≤ 3 n T , since vertices are unique and shared between triangles.The triangular mesh is oriented by the triangle normals {n i } n T i=1 pointing out of T .If T encloses a bounded volume, the interior of T is denoted T .
Let C ⊂ R 3 be a Cartesian cell, i.e., an axis-aligned cuboid consisting of its boundary ∂C and interior C .Then, ∂C = 6 i=1 F i , where {F i } 6 i=1 , are the rectangular cell faces.A cell C also has twelve edges {E i } 12 i=1 and eight vertices {V i } 8 i=1 .In this paper, we consider the problem of calculating the percentage or solid volume fraction of C inside T .The solid volume fraction α is given by where V (X) denotes the volume of a subset X of R 3 .We also consider the percentage or solid area fraction of each cell face F i inside T .The solid area fraction β i corresponding to face F i is given by where A(X) is the area of a two-dimensional object X.A cell C intersecting a triangular mesh T is seen in Fig. 1a, and T C is seen in Fig. 1b.More generally, we consider multiple meshes T 1 , . . ., T N , where α and β i are the fractions of the cell and faces inside N i=1 T i .Each mesh is allowed to include overlapping triangles or double surfaces.To define the concept of a double surface, T i=1 be proper triangle meshes (without geometric degeneracies).Let S 1 ⊂ T 1 and S 2 ⊂ T 2 be such that A typical double surface is seen in Fig. 2.

Algorithms
In this section, the two proposed methods to calculate the solid volume fraction α and solid area fractions β i are presented.Both methods find a polyhedron P that represents T C and calculate α and β i by and i ∈ {1, . . ., 6}.The first method exactly calculates T C .It is presented in Section 4.3.The second method is an approximate method based on approximating In this first phase, we require that the triangle meshes are non-degenerate.This restriction is removed in Section 5, where we consider double surfaces and multiple meshes.Pseudo code for the algorithms is presented, with some helper methods outlined in Appendix A. The algorithms in the pseudo code are sometimes simplified, but cover the key steps of a successful implementation.

Preliminaries
Before we proceed, we define some important concepts: face polygon, cell polygon, and intersection point.We use a notation similar to that of [1,2].Given a polyhedron P and cell faces F i , i ∈ {1, . . ., 6}, the sets P F i • , are called face polygons, and we say that ∂P C consists of a number of cell polygons (see Fig. 3).Here, P F i • is the closure of the interior of the two-dimensional set P F i and ∂P C is the closure of the three-dimensional set ∂P C .Now, let S be a subset of a two-dimensional space, and let ∂ 0 (S) be the extreme points of S. The set is called the intersection points between C and T (see Fig. 4 where, for a particular triangle T i , ∂ 0 T i F j is marked by spheres).The set are referred to as the triangle vertices of T inside C .The divergence theorem is used to calculate the volume of a polyhedron P .
Applying the theorem to the vector field F(x, y, z) = 1 3 (x, y, z), we get where c i is the centroid, ni the unit normal, and A i the area of the i:th face S i of the surface of P .To see this, note that ∇ • F = 1 and Fig. 4 Intersecting a triangle and a Cartesian cell results in intersection points, here marked by spheres The areas and centroids in (7) are calculated as in [8]: and where c = (c x , c y ) and the polygon has vertices (x 0 , y 0 ), . . ., (x n−1 , y n−1 ), (x n , y n ), with (x n , y n ) = (x 0 , y 0 ).A change of coordinates is performed so that (x i , y i ) are the coordinates of vertex i in a basis for the two-dimensional space defined by the polygon plane.

Triangle-cell intersections
To calculate the area and volume fractions, we first need to find the intersection points between the triangle mesh and the hexahedral cell according to (5).What we basically want to do is to intersect a number of triangles and an axis-aligned hexahedron in a robust way.By robustness, we mean that we want to locate the points of intersection, we do not want to miss an intersection due to numerical imprecision, and we want related intersection tests to give consistent results.The latter requirement could fail for example if an edge that is shared between two triangles is not represented consistently between the two triangles [12].We also want unique intersections, meaning that an intersection point on a triangle edge is shared between the two triangles meeting at the edge, and that an intersection point in a triangle vertex is shared between all triangles that meet in the vertex.
The above robustness criteria makes the triangle-box intersection problem hard to solve.Determining whether a triangle and an axis-aligned box intersect or not is a well-studied problem, discussed for example in [3].Robustness issues are discussed by among others [12,24].It gets more complicated when the locations of the intersection points are also required.One approach to this is to perform a sequence of simpler intersection tests, where the triangle edges are intersected with the rectangular faces of the box, and the box edges are intersected with the triangle.Such line segment-rectangle and line segment-triangle intersection tests are common in several applications, often build upon the parametric representation of the geometric objects, and summarized for example in [12].
To get unique intersections, we have introduced vertices and edges that are shared between triangles.Similarly, we have introduced box vertices and box edges that are shared between the rectangular faces of the box.These joint edges and vertices are then used in a sequence of simpler intersection tests, as described above.Before the intersection points are calculated numerically, we, as far as possible, discretely determine how many intersections there are.This is done using coordinate comparisons similar to the outcodes described in [1,2].The coordinate comparison approach is appealing since the cell is axis-aligned, and the x-, y-, and z-coordinates can be considered separately.Neither are there any numerical issues with the coordinate comparisons, since no floating point calculations have to be performed.

Exact method
In the exact method, we find the exact polyhedron of intersection T C (see Fig. 1b).We apply (1)-( 2) and ( 7) to calculate the exact area and volume fractions of C .The main steps in the method can be summarized as follows: 1. Find the intersection points I between the mesh T and the cell C 2. Find the triangle vertices I 0 of T inside the cell C 3. For i ∈ {1, . . ., 6}, connect intersection points in I F i to polygons ij , and calculate the area fraction β i as the total area of the polygons ij divided by the area of the face F i 4. For each triangle T i intersecting the cell C , connect I I 0 T i to a cell polygon T i 5.The face and cell polygons ij and T i in steps 3-4 define the polygonal faces of the polyhedron T C .Calculate volume fraction of T C according to ( 1) and ( 7).
The intersection points in step 1 are given by ( 5), and the triangle vertices I 0 in step 2 by (6).More details about the construction of polygons in steps 3-5 are found in Section 4.3.1.The main steps in the algorithm are illustrated in Fig. 5.
Pseudo code for an algorithm that computes the exact solid volume and area fractions is presented in Algorithm 1.In the algorithm, we introduce the notation T = T (C ) for the set of triangles T ⊆ T that are coplanar to some cell face for some k ∈ {1, . . ., 6}, and I 0 T i = ∅.This means that an edge or vertex of T i ∈ T ↑ touches a cell face, but the intersection of the triangle and the interior of the cell is empty.
The method normal = UnitNormal(polygon) returns a unit normal of the input polygon.The method UnitNormals(polygons) returns a list of normals, one for each polygon in polygons.The method Area(polygon) and Centroid(polygon) are implementations of ( 9) and (10), respectively.Centroids(polygons) calls Centroid(polygon) for each polygon in polygons.The method polygons = ConnectPolygons(intersection points, target facet) connects the points in points to polygons.It is described in Section 4.3.1.
In the algorithm, we do not create cell polygons from triangles in T or T ↑ , and we remove intersection points in T \ (T T ↑ ).We do not create cell polygons from coplanar triangles to avoid that they contribute with the same area both to cell and face polygons, which would result in an error in (7).Since we do not create cell polygons from coplanar triangles, the coplanar triangles have to be accounted for in the face polygons.Instead of resolving every coplanar triangle in the face polygon algorithm, we remove from I the set {x ∈ I : x / ∈ T \ T }; that is the intersection points that are only in coplanar triangles.If a cell face F k has no intersection points after the removal, we know that the area fraction of F k is 0.0 or 1.0.To find out which, we check if a point x ∈ F k is on the inside or outside of T .This is done using a method based on angle weighted pseudo normals, similar to that of Baerentzen and Aanaes [6].Since we do not consider triangles in T , we must also disregard triangles in T ↑ .To see why, consider Fig. 6a, where a cell intersected by a triangulation T including triangles T 1 ∈ T ↑ and T 2 ∈ T is seen in cross section.One of the intersection points marked by spheres belongs to both T 1 and T 2 , and the other one only belongs to T 2 .The intersection point that is only in T 2 is removed since T 2 ∈ T .Assuming symmetry, the area fraction of the top cell face is 1.0.Not removing the intersection points in T 1 would erroneously give an area fraction less than 0.5.Therefore, we also remove from I the set {x ∈ I : x / ∈ T \ (T T ↑ )}, the points that are only in coplanar triangles or triangles pointing out of the cell.
As seen in Fig. 6b, we can not remove all x ∈ I T .Again assuming symmetry, the top area fraction is now approximately 0.6, since there is a face polygon that in cross section is the line between the two intersection points.We need to keep the intersection point that belongs to both T 1 and T 2 to construct this polygon.The difference to Fig. 6a where we removed both intersection points is that T 1 now does not belong to T ↑ but points into the cell.

Connecting intersection points to polygons
In this section, we describe how the intersection points in I I 0 are connected to polygons.It is critical that the intersection points are processed in correct order so that ( 9)-( 10) can be used to calculate area and centroid.As in [11], the connectivity of the triangle mesh and the cell is used to achieve this.While [11] briefly states that Fig. 6 Cross sections of a cell intersected by a triangulation T including triangles T 1 and T 2 the face polygons are found through connectivity criteria such as common faces or common triangles, we describe how this is done.
There are two types of polygons that need to be created.The first type is the cell polygon, which is the intersection between a triangle T i and T .The second type is the face polygon that for each cell face F i , i ∈ {1, . . ., 6}, has vertices I F i .While [1,2] uses the Sutherland-Hodgman algorithm [25] to clip each triangle T i against the cell faces to create the cell polygons T i , we use the connectivity to connect all points in I I 0 T i .This will be convenient when detecting double surfaces (see Section 5.1.1).To create the face polygons { ij }, we use the connectivity to connect all points in I F i .
Both types of polygons can be constructed in a similar way, which can be described after introducing target facets (F T ) and search facets (F S ).A target facet is the triangle T i when the cell polygon T i is constructed, and the cell face F i when the face polygons { ij } are constructed.The goal is to connect all intersection points on the target facet to one or more polygons.To do this, we need search facets.A search facet is a triangle, a triangle edge, a cell face, or a cell edge in which to search for the next vertex of the polygon that is under construction.
Given a correct current search facet, there should always be a uniquely determined intersection point to take as next polygon vertex.The main idea is to find the next vertex on the current search facet, and then change search facet.The next search facet is determined by the location of the newly added vertex with respect to the triangle mesh and the cell.This procedure, taking the next vertex and updating the search facet, is repeated until the polygon is completed, and until all intersection points on the target facets have been assigned to a polygon.Examples are given in Figs.7 and 8, where a simple face polygon and a simple cell polygon are constructed.
In Algorithm 1, the recently described polygon connection step is performed by the method ConnectPolygons(intersection points, target facet).It is given in pseudo code in Algorithm 3 (see Appendix A.1).

Approximate method
In the approximate method, T C is approximated by a convex polyhedron P with convex polygonal faces.P is defined by a least squares plane fitted to the points in I I 0 and a normal of the plane.The plane splits C in two parts, P and C \ P .The main steps in the method can be summarized as follows: 1. Find the intersection points I between the mesh T and the cell C 2. Find the triangle vertices I 0 of T inside the cell C 3. Fit a least squares plane π to the intersection points I and the vertices found in steps 1-2 4. Identify P and C \ P 5.For i ∈ {1, . . ., 6}, calculate the area of P F i and the area fraction β i 6. Calculate the volume of P and the volume fraction α.
We now comment on the details of this procedure.The intersection points in step 1 are defined by (5), and the triangle vertices I 0 in step 2 by (6).In step 3, a total where p 0 is a point in the plane and n is the plane normal.We can take [13] Let A be the matrix Then, (11) can be written where || • || denotes the Euclidean norm.The minimum can be found by computing the singular value decomposition of A. One then finds the minimum n as the third column of V .To see this, note that where λ 1 ≥ λ 2 ≥ λ 3 are the squared singular values of A, and z = V T n.The second equality in (16) follows since U is an orthogonal matrix, and the last equality since S(i, i) = √ λ i and S(i, j ) = 0, i = j .The expression in ( 16) is thus minimized when z = e 3 = (0, 0, 1), or equivalently when the third column of V .
In step 4, the correct sign of the plane normal n π is critical to determine what is P and what is C \ P .One heuristic to determine the sign is to compare n π to the area weighted average n aw of all normals of triangles intersecting the cell, where ni is the unit normal of T i .If n aw •n π > 0, the sign of n is assumed correct; otherwise, it is changed.A second alternative is to replace A(T i ) by A(T i C ) in ( 18).This will have advantages when the mesh includes double surfaces (see Section 5), but is more expensive since T i C has to be found.In addition, T i C will in general have a more complex shape than T i , which makes the area calculation more expensive.
For steps 5-6, we need the cell and face polygons of P (see Fig. 3).The only cell polygon is convex and defined by It is constructed from the points in the intersection between the plane π and the edges of C .The face polygon on cell face F i is denoted by i , and defined by It is constructed by connecting the points in I π F i with the cell vertices on the correct side of the plane π.The correct side is determined by the plane normal n π .What remains is to calculate the volume of the polyhedron according to (7) and divide it by the cell volume to get the volume fraction α.The face area fraction β i is given by the area of i divided by the area of F i .The main steps in the algorithm are illustrated in Fig. 9. Pseudo code for an algorithm that computes P and hence the solid volume and area fractions of C is described in Algorithm 2. The method plane = LeastSquaresPlane(points) on line 3 is an implementation of the least squares plane algorithm described in ( 11)- (17).The method CorrectNormal(plane, triangle mesh, cell) on line 4 corrects the normal n π of the input plane according to the sign of n aw • n π as described above, where n aw is given by (18).The method FindCellVerticesInPolyhedron(cell, plane) returns a list of the vertices of the input cell on the correct side of the input plane, where the correct side is determined by the plane normal.The method normal = UnitNormal(polygon) returns a unit normal of the input polygon.The methods area = Area(polygon) and centroid = Centroid(polygon) are implementations of ( 9) and (10), respectively.Finally, polygon = ConvexPolgon(points, normal) takes as input a list of polygon vertices and the normal of the polygon, and returns a list of the vertices sorted in clockwise or counterclockwise order.Pseudo code for this method is presented in Algorithm 11 in Appendix A.2.

Geometric complications
The algorithms described in Section 4 calculate solid area and volume fractions for a non-degenerate triangle mesh.In practice, it is not unusual that the mesh is degenerate in some sense.It can include T-vertices (hanging nodes), gaps, and cracks, or consist of overlapping meshes.The latter leads to what we call double surfaces, the intersection of overlapping triangles.In this section, it is described how the algorithms in Section 4 are modified to handle double surfaces with normals of the overlapping triangles pointing in opposite directions.It is also shown that the algorithms handle multiple meshes and split hexahedra.These are steps towards a method that calculates the solid area and volume fractions of a degenerate triangle mesh including overlapping triangles, overlapping meshes, hanging nodes, and gaps.Such triangle meshes are common in industrial applications.

Double surfaces
Recall the definition of a double surface from Section 3. In Fig. 2, we saw a double surface including triangle normals with opposite directions.In Fig. 10, there is another, schematic, example of a mesh including a double surface.Two typical cases of intersection between a cell and the double surface are seen in Fig. 11.

Algorithm 2 {α, {β i }} = CalculateApproximateAreaAndVolumeFraction(T , C )
Data: A triangle mesh T and a hexahedral cell C .Result: α: approximate fraction of C inside T .{β i }: approximate fractions of F i inside C , where F i , i ∈ 1, . . ., 6, are the faces of C .When T 1 and T 2 are handled as separate meshes both algorithms correctly accounts for the double surface as explained in Section 5.2.To solve the degenerate case when the two meshes are treated as one, the algorithms are modified as described below.In the rest of the section, we consider double surfaces on a single mesh T = T 1 T 2 .

Exact method
The exact method in Section 4.3 is modified to handle both cases in Fig. 11.It is not necessary that the double surface is axis-aligned.The cell polygons are unproblematic since a cell polygon is created by connecting all intersection points in a particular triangle.The contribution to the volume fraction from the triangles in the double surface will cancel out in (7).Possible problems occur when the face polygons are overlapping triangles Fig. 10 Cross section of two triangle meshes T 1 and T 2 glued together.The arrows represent triangle normals pointing out of the interior of the respective mesh.The overlapping triangles give rise to a double surface Fig. 11 Cross section of typical cases of intersection between a cell and a double surface created.Then, it has to be made sure that intersection points from the underlying meshes T 1 and T 2 are connected to separate polygons.The solution is to use the connectivity of T to extract information about the underlying meshes T 1 and T 2 , so that the two parts of T can be handled separately.
Figure 12 shows a cell face F i and intersection points I F i from Fig. 11a.The intersection points are marked by spheres and boxes.The spheres represent I F i T 1 , and the boxes represent I F i T 2 .The whole cell is inside, so the face polygons { ij } should cover F i .In Fig. 12b, the points in I F i T 1 have been connected to one polygon i1 and the points in I F i T 2 to a second polygon i2 , where F i = i1 i2 as expected.When the search facet is a triangle, there is no ambiguity in the choice of the next polygon vertex.When the next intersection point is searched on a face side, as in Fig. 12a, the presence of the double surface becomes apparent.Then, it has to be made sure that only intersection points from the correct part of T can be taken.Normally, the next intersection point on the face side is the one closest to and on the correct side of the current intersection point (see Algorithm 7 in Appendix A.1).At a double surface, there are two intersection points at the same position.The solution Fig. 12 Face polygon connection at a double surface.Intersection points marked by spheres come from T 1 and intersection points marked by boxes come from T 2 is to make sure that an intersection point in a triangle with normal pointing against the current intersection point cannot be taken.This is a check that the polygon that is created is actually located on the inside of the underlying mesh T 1 or T 2 , since the triangle normals point of T 1 and T 2 (at a double surface the surface normal of T is undefined).To account for the presence of double surfaces, Algorithm 7 has to be modified according to this observation.The result is presented in Algorithm 10 in Appendix A.1.

Approximate method
To see that the approximate method in Section 4.4 is inaccurate for meshes with double surfaces, consider Fig. 11a.The cell is intersected by a double surface, which will make all intersection points lie in a plane in the middle of the cell.The unmodified algorithm takes this plane as the approximating plane and gives a volume fraction around 0.5, though the whole cell is inside.The solution is to check if all normals of triangles intersecting the cell point in the same or opposite direction, and if the triangles lie in the same plane.The volume fraction is then set to 1.0.
In Fig. 11b, a cell is intersected both by a double surface and a regular part of T = T 1 T 2 .The intersection points at the double surface will erroneously affect the plane.This problem is not handled as a special case in the approximate method.
As pointed out in Section 4.4, it is better to use the normal weighted with the area of the part of a triangle inside the cell than the area of the whole triangle when determining the sign of the least squares plane normal.To see an example of why, consider Fig. 13, where a cross section of a cell intersected by a triangle mesh including triangles T 1  1 , T 1 2 , and T 2 1 is seen.Parts of T 1 2 and T 2 1 overlap in a double surface.The dashed line and normal represent the approximating plane π in the approximate method.Now assume the area of T 1  2 is much larger than that of T 1 1 and T 2 1 .The comparison normal can then be approximated by that of T 1  2 , and we would erroneously take the wrong sign for the normal of π by the sign test of the dot product.If we weight with the area inside the cell, the double surface will be canceled out and will not affect the comparison normal.
Fig. 13 Cross section of cell intersected by a double surface created by T 1 2 and T 2 1 .The approximating plane π in the approximate method is marked by a dashed line

Several meshes
Both algorithms in Section 4 can be applied on multiple meshes T 1 , . . ., T N by considering each T j separately and summing the area fractions β ij , i ∈ {1, . . ., 6}, and volume fraction α j from each mesh according to and

Small self-intersections or overlaps
When floating points are used, numerical issues could make the triangles in a double surface slightly overlap or be slightly separated, as in Fig. 14.The proposed algorithms work even if the triangles in the double surface are slightly separated and not overlapping.However, an extra condition has to be added if small overlaps are to be handled.
A limit ε for the maximal overlap allowed can be introduced and used to indicate when a double surface is found.Two intersection points at a separation distance smaller than ε are assumed to lie on a double surface if they also belong to triangles with opposite normals.Due to numerical issues, it will be necessary to have a limit on the oppositeness of normals.In the following, we assume that two normals are opposite if π − θ < δ, where θ is the angle between the normals and δ is a prescribed limit.
In the approximate method, it is enough to use the limits to determine when two triangles with opposite normals are coplanar.In the exact method, the limits are used in the face polygon connection step to indicate whether a certain intersection point can be taken as the next vertex in the face polygon or not.In the following, we discuss the treatment of double surfaces in the exact algorithm.Fig. 14 Cross section of typical cases of intersection between a cell and a double surface.The double surface includes small overlaps or small separations between triangles When a double surface intersects a cell, as in Fig. 11, we wish to connect the intersection points on T 1 to one polyhedron of intersection, and the intersection points on T 2 to another polyhedron of intersection as described in Section 5.1.The result after a correct polygon connection step is visualized in Fig. 15, where the different striped patterns correspond to the different polyhedrons of intersection, one for T 1 and one for T 2 .
For an overlapping double surface, the goal is to get to the result in Fig. 16.The two meshes T 1 and T 2 still contribute with one polyhedron of intersection each, but there will now be an overlap between the polyhedrons.This overlap corresponds to the overlap of the double surface.The area and volume of the overlap is counted twice in the solid area and solid volume fraction calculations.To get to the result in Fig 16, the face polygon connection step has to be modified, or the result will.be as in Fig. 17 where only a fraction of the correct polyhedrons of intersection is found.
The problem with the cases in Fig. 17 is that the algorithm cannot distinguish the overlapping double surface from a thin or sharp-edged part of the mesh.Examples of thin and sharp-edged meshes are found in Fig. 18a and b, respectively.A sharp edge is in this case formed when two triangles meet at an angle less than δ.For the cases in Fig. 18, it is correct to connect the two close intersections at the left and right cell faces, while it is wrong to do the same thing for the cases in Fig. 17.
To avoid the problem demonstrated in Fig. 17, the maximal overlap limit ε is used as an upper limit on how close two nearby vertices of a face polygon can be.If the distance between the current vertex and a candidate for the next vertex is closer than ε, and the two vertices belong to triangles with opposite normals; the candidate can not be taken as the next vertex.With this modification of the polygon connection algorithm, the result will be as in Fig. 16.A consequence is that sharp-edged meshes and meshes thinner than ε cannot be handled.
The above modification of the exact algorithm is necessary only if the double surface belongs to a single mesh T = T 1 T 2 .If T 1 and T 2 are handled separately, the double surface is resolved by finding the face polygons and polyhedrons of intersection for one mesh at a time as described in Section 5.2.The area and volume fractions from the different meshes are then added.
For separate meshes T 1 and T 2 , the overlaps that can be handled are not restricted to double surfaces.It is possible to handle arbitrary overlaps, as long as it is reasonable to count the overlap volume and area twice.Fig. 15 Cross section of typical cases of intersection between a cell and a double surface.The intersections have been connected into the correct polyhedrons of intersection, one for T 1 and one for T 2 Fig. 16 Cross section of typical cases of intersection between a cell and an overlapping double surface.The intersections have been connected into the correct polyhedrons of intersection, one for T 1 and one for T 2

Split hexahedra
When the hexahedral grid is coarse or the geometry is thin or detailed, a cell can be cut in several disconnected parts by one or more polyhedrons of intersection.This leads to the formation of so called split cells, or split hexahedra.An example of a split cell in two dimensions is found in Fig. 19.
The exact method easily handles most cases of split hexahedra by construction.The cell polygons are treated triangle by triangle as described in Section 4.3.1 and in Algorithm 3, which is not complicated by split hexahedra.The face polygons are uniquely determined by the connectivity of the triangle mesh under the same assumptions of self-intersections as for double surfaces (see Section 5.3), and not either complicated by most forms of split hexahedra.The face polygon construction algorithm described in Section 4.3.1 and in Algorithm 3 finds as many polygons of intersection as needed until all intersection points have been assigned to a polygon.
One case that is not covered by the algorithm described this far is when the triangle mesh intersects a cell face and forms one or more polygonal "holes," as in Fig. 20.This happens when a connected part of the triangle mesh intersects the face without intersecting the boundary of the face.The algorithm cannot determine what is the inside and what is outside of the "hole," since the triangle normals are not used in the connection step.
To handle polygonal holes, signed polygon areas are used.The sign of the area is determined by the normals of the triangles that form the polygon.The area is positive Fig. 17 Cross section of typical cases of intersection between a cell and a double surface.The intersections have in both cases been connected erroneously to only one polyhedron of intersection Fig. 18 Meshes that cannot be distinguished from an overlapping double surface in the method based on a maximal limit of the allowed overlap if the normals point out of the polygon, and negative if the normals point into the polygon.To determine if the normals point into the polygon or not, a method similar to [6,14] is used.The closest polygon vertex v to an arbitrary face edge is first found, the vertex is projected onto the edge, and the vector from the projection v e to the vertex is formed.The sign of the dot product (v e − v) • n v determines if the triangle normals point into or out of the polygon according to where sgn(A) is the sign of the polygon area and n v is the average of the triangle normals of the triangles meeting in v.An illustration is found in Fig. 21.
If none of the face polygons intersect the boundary of the face, the area of the whole face has to be added to the total area of intersection if the "outermost" polygon is a hole.This is the case in Fig. 20b.The outermost polygon is found by locating the vertex closest to an arbitrary face edge, and identifying the corresponding polygon.

Numerical examples and results
Three test cases have been studied for evaluation of the algorithms.The cases were designed to test the algorithms on double surfaces, numerical or small overlaps, and Fig. 19 The triangle mesh T splits the cell C into four parts.Two of the four parts are disconnected polyhedrons of intersection Fig. 20 Two triangle meshes T 1 and T 2 intersecting a cell C in a way that splits the hexahedral cell in three parts, without intersecting the boundary of any cell face split hexahedra.In Section 6.1, double surfaces and split hexahedra is tested with a geometry representing a heat sink used for air cooling of CPU:s.In Section 6.2, we test the exact method for double surfaces and small overlaps on a geometry consisting of two cylinders whose surfaces should meet in a double surface, but the surfaces are not perfectly matching.In Section 6.3, the geometry consists of four cylinders placed inside each other.It is used to test the exact method for split hexahedra and holes.
The algorithms were implemented in the C++-based multiphase flow framework IBOFlow [16].All computations were carried out on a machine equipped with a 3.50-GHz Intel Core i7 processor (5930K) and 64 GB of RAM (1.066 GHz DDR4).

Heat sink with a double surface
The exact method is validated and compared to the approximate method on a geometry including a double surface.The geometry represents a heat sink and is taken from an industrial application.The whole mesh is seen in Fig. 22a, and a zoom in on the part of the mesh that is marked by a circle is shown in Fig. 22b.The triangle pattern    At each refinement level, the smallest cell size was halved and eight new cells were formed.The test setup and the solid volume fractions for the base grid is seen in Fig. 23a.The solid volume fractions for grids with one, three, and six refinements are seen in Fig. 23b, c, d, respectively.The exact method was used in all four cases.
Results of the grid study are presented in Fig. 24a.The total volume of the interior of the geometry was calculated by running the exact and approximate volume fraction algorithms for each cell.This was repeated for each refinement level.The calculated volume was compared to the real volume 1.45 • 10 −5 m 3 , which was found by applying (7) to the whole triangle mesh.
A plot of the CPU time against the number of fluid cells is seen in Fig. 24b.The times are the average results from ten runs.A Cartesian grid with initial cell size ( x, y, z) = (4.16•10−3 , 4.16•10 −3 , 4.16•10 −3 m) was used for the time study.It was refined up to five times around the triangular mesh, resulting in a finest grid with 8, 438, 270 cells of size ( x, y, z) = (1.30• 10 −4 m, 1.30 • 10 −4 m, 1.30 • 10 −4 m) or larger.In summary, the exact method is independent of cell size, while the approximate method is second-order accurate.The exact method is a constant factor slower than the approximate method.The constant factor is close to 2.0 when the number of cells is large.
The grid study in Fig. 24a and the test setup in Fig. 25 indicate that the exact method handles split hexahedra well.In Fig. 25a, the initial grid has been refined three times, and several cells are split in two or more parts by the triangle mesh.One of the split cells is marked green.A zoom in on this cell and the polyhedrons

Overlapping cylinders
Two identical co-located cylinders of radius 0.5 m were triangulated, joined to one, and intersected by an axis-aligned hexahedral cell of dimensions ( x, y, z) = (0.23 m, 0.208 m, 0.312 m).The co-located cylinders simulate a double surface.Different triangulation schemes were used for the two cylinders so that the resulting triangle meshes slightly overlapped instead of forming a perfect double surface.The resulting setup is seen in Fig. 26.In Fig. 26a, it is seen how the hexahedron intersects the triangle meshes, and in Fig. 26b, we have zoomed in on parts of a cross section of the triangle meshes to show that the triangles overlap.The algorithm handles the overlapping double surface according to Section 5.3, where ε = 0.01 x m and δ = 2 π 45 rad is chosen.The polygons and polyhedrons of intersection are correctly connected, as demonstrated in Fig. 27.One of the polyhedrons is colored red and the other black, and they belong to one cylinder each.Both overlaps and gaps had to be handled to construct the face polygons.This is seen from the face polygons for face 2 (y-top) in Fig. 28.For visibility, the overlaps and gaps have been exaggerated in the figure.At the bottommost face edge, the factual overlap between the triangles is 0.002 x m, which is within the tolerance ε.The angle θ between the triangle normals satisfies π − θ = 1.23 • 10 −1 rad, which is less than δ.In the first setup, the centerline of each cylinder was aligned with the y-axis and passing through the center of the xz-faces of the hexahedron (Fig. 29).In the second setup, the cylinders were rotated π 4 rad about the x-axis from their initial position (Fig. 30).For these setups, the analytic volume fraction is to three significant figures α 1 = 8.95 • 10 −2 , α 2 = 1.13 • 10 −1 , and the area fractions are

Split hexahedra
where the superscripts 1 and 2 correspond to the first and second test cases, respectively.The area fractions in ( 24) and ( 25) are calculated from the formulas for the area of a circle and an ellipse, respectively.The first volume fraction, α 1 , is simply calculated from the formula for the volume of a cylinder.The second volume fraction is  where is the volume of intersection of the cylinder with radius r i and the hexahedron.In ( 27), V i is the sum of the volume of three parts.The first part is the largest cylinder C that fits into the intersection between the rotated cylinder and the cell.The remaining parts are the two identical bodies that are formed as the difference of the whole volume of intersection and C. The volume of the two identical bodies was calculated from the definition of the volume integral.The area and volume fractions calculated with the proposed algorithm agree with the theoretical results.

Discussion and conclusion
Two methods for calculation of solid area and volume fractions of the intersection between degenerate triangle meshes and a hexahedral grid have been proposed.The algorithms have been implemented in the multiphase flow framework IBOFlow and will be used to improve the accuracy in the calculation of fluxes between fluids and solids.The solid area and volume fractions indicate how much of each fluid cell that is intersected by the solid.The exact algorithm is the main result that handles all geometric complications addressed in this paper.The approximate method is intended for highly resolved grids, for which it is reasonable to approximate the cell-mesh intersection with a plane.There are some comments on what could have been done differently in the two methods.In the exact method, since both hexahedra and triangles are convex, the cell polygons are also convex.This could be used in an alternative method to connect the intersection points to cell polygons, for example as described in Section 4.4 and Appendix A.2.In the approximate method, since it is known that the vertices to be sorted are the result of intersecting a plane with a hexahedral cell, an alternative is to use the connectivity of the cell to connect the polygons.A method similar to, but simpler than, the one used in the exact method in Section 4.3 and Appendix A.1 could also have been adopted.The advantages and drawbacks of these alternatives have not been investigated.
A limitation of the current implementation is the numerical computation of the triangle-cell intersection points, which could fail due to numerical round-off.In case this would happen, we currently use the approximate algorithm as a backup.Since the approximate algorithm cannot handle all cases, a better alternative would be to switch to exact arithmetics when the floating point precision algorithms fail.An even better approach would be to improve the discrete method for determining the number of intersection point and their location.If this could be done perfectly, all numerical problems would be eliminated.
A next step towards a geometrically robust method is to extend the exact algorithm to handle more general triangle overlaps and hanging nodes.Overlaps could be found by using the fact that each underlying triangle mesh is handled separately.Hanging nodes could be handled by introducing triangle subedges, and account for these in the polygon connection algorithms.

A.2 Approximate method
In the approximate method outlined in Algorithm 2 in Section 4.4, we need to sort the vertices of a convex polygon in clockwise or counterclockwise order.This is performed by the method ConvexPolygon in Algorithm 11, which sorts a list of vertices given the normal of the polygon.Pseudo code for SortVerticesRecursive, which performs the sorting and returns a list with the sorted input points, is presented in Algorithm 12. Pseudo code for the method PartitionVertices used in Algorithm 12 is presented in Algorithm 13.An overview of the methods and where they are called is presented in Table 2.

Fig. 1 Fig. 2 2 TC
Fig. 1 Intersection between a triangle mesh T and a Cartesian cell C

Fig. 3
Fig.3The polyhedron in Fig.1bseen as a composition of face polygons and cell polygons.A face polygon is the intersection between a cell face F i and the polyhedron P .A cell polygon is a face of P in the interior of the cell C

Fig. 5
Fig. 5 Description of how the polyhedron of intersection is found in the exact method

Fig. 7
Fig. 7 Description of how a face polygon is created according to Algorithm 3

Fig. 8
Fig. 8 Description of how a cell polygon is created according to Algorithm 3

Fig. 9
Fig. 9 Description of how the polyhedron of intersection is found in the approximate method (a) The area is positive since n v points out of the polygon.(b) The area is negative since n v points into the polygon.

Fig. 21
Fig.21 The sign of the dot product (v e − v) • n v is used to determine if the polygon is a hole or not (a) Test geometry representing a heat sink.(b) Zoom in on the part in the left figure marked by a circle.The overlapping triangles indicate that there is a double surface.

Fig. 22
Fig.22The test case (a) The test setup in intersection and the resulting solid volume fraction for the initial grid consisting of 8 hexahedral cells.(b) Resulting solid volume fraction for a grid that has been refined three times.(c) Resulting solid volume fraction for a grid that has been refined once.(d) Resulting solid volume fraction for a grid that has been refined six times.

Fig. 23
Fig. 23 Test setup and resulting solid volume fraction for different number of refinements of the initial grid (a) Grid convergence of volume fraction algorithms run on the heat sink case.(b) CPU time of volume fraction algorithms run on the heat sink case.

Fig. 24
Fig. 24 Results from grid and CPU time studies (a) The test setup for three refinements.One cell that is split in several parts by the triangle mesh is colored green.(b) The resulting polyhedrons of intersection that split the green cell in several parts.

Fig. 25
Fig.25 The test setup for three refinements results in a number of split hexahedra (a) Two overlapping cylinder triangle meshes are intersected by a hexahedral cell (blue).The arrows represent triangle normals pointing out of the interior of the meshes.Only half of the triangle meshes are shown.
(b) Part of a circular cross-section of the triangle meshes, which shows that the triangle meshes slightly overlap.The overlaps have been exaggerated in the figure.

Fig. 26
Fig.26 Setup that tests the algorithm for a double surface formed by slightly overlapping, non axisaligned triangles

Fig. 27
Fig.27 Polyhedrons of intersection between two overlapping cylinders and a cell.The polyhedron with red border corresponds to one cylinder, and the polyhedron with black border corresponds to the other cylinder.Together the polyhedrons cover the whole of the cell, as expected.The overlaps and gaps have been exaggerated in the figure

( a )
The test setup.(b) The cylinder triangle meshes in intersection, with triangle normals pointing out of the interior of each mesh.

Fig. 29
Fig. 29 Four axis-aligned cylinders with the same center line and different radii intersecting an axisaligned hexahedron

4 ← ∅ 5 v 0 ← 6 F 7 v i ← v 0 8 t 0 ← true 9 initialize d 10 while 14 v 15 mark v i+1 as taken 16 if target facet is a triangle then 17 F 19 {F 20 let v i+1 be the next vertex of 21 v i ← v i+1 22 t 0 ← false 23 i 2 if 4 else 5 next 2 if 3 /else 8 // current search facet is a cell face 9 if no untaken intersection point in the current search facet then 10 if current intersection point is in a face side then 11 next 13 / 15 next vertex ← the untaken intersection point in current search facet 16 else 17 / 8 / 6 d p ← |s p − s| 7 if d p < d then 8 d ← d p 9 next vertex ← p Algorithm 8 4 if 5 next search facet ← t 6 finished ← true 7 else 8 / 15 next search facet ← s 16 if 25 /
Data: A triangle or a cell face, and the intersection points located on that triangle or cell face.Result: The cell polygon or face polygons associated with the triangle or cell face. 1 polygons ← ∅ 2 i ← 0 3 while i < number of elements in intersection points do first untaken intersection point in target facet S ← GetSearchFacet(target facet, v 0 ) (v i = v 0 ) ∨ t 0 do 11 if target facet is a triangle then 12 v i+1 ← GetNextVertexInPolygon v i , F S 13 else i+1 ← GetNextVertexInPolygon v i , F S , d, t 0 S ← ChangeSearchFacet(v i+1 , F S ) 18 else S , d} ← ChangeSearchFacetAndUpdateSearchDirection(v i+1 , F S , d) ← i + 1 24 append to polygons Algorithm 3 polygons = ConnectPolygons(intersection points, target facet) Data: A triangle or a cell face, and the intersection point representing the first vertex in a cell polygon or face polygon associated with the triangle or cell face.Result: The search facet on which to find the next vertex in the cell polygon or face polygon. 1 if target facet is a triangle then intersection point is in a triangle edge or a triangle vertex then 3 next search facet ← the triangle edge or triangle vertex in which the intersection point is included search facet ← the face side or face corner in which the intersection point is included 6 else 7 next search facet ← the triangle in which the intersection point is included current search facet) Data: The current vertex in a cell polygon under construction, and the search facet in which to find the next vertex in the cell polygon.Result: The next vertex in the cell polygon. 1 if current search facet is a triangle edge then no untaken intersection point in the current search facet then / current intersection point is in a triangle vertex 4 next vertex ← the untaken intersection point in the other triangle edge 5 in which the current intersection point is included 6 else next vertex ← the untaken intersection point in the current search facet 7 vertex ← next untaken intersection point in the second cell face in which the current intersection point is included 12 else / current intersection point is in a face corner 14 next vertex ← next untaken intersection point in one of the other two faces in which the current intersection point is included 15 else 16 next vertex ← the untaken intersection point in the current search facet Algorithm 4 next search facet = GetSearchFacet(target facet, intersection point) Algorithm 5 next vertex = GetNextVertexInPolygon(current intersection point, next vertex ← GetNextVertexInPolygonOnFaceSide(current intersection point, current search facet, current search direction) ← next untaken intersection point in one of the other triangles in which the current intersection point is included 14 else / current search facet is a face side 18 next vertex ← GetNextVertexInPolygonOnFaceSide(current intersection point, current search facet, current search direction) Algorithm 6 next vertex = GetNextVertexInPolygon(current intersection point, current search facet) Data: The current vertex in a cell polygon under construction, and the search facet in which the current vertex was found.Result: The search facet in which to find the next vertex in the polygon 1 if current intersection point is in a triangle vertex v then 2 next search facet ← the triangle edge including v that is not the current search facet 3 else if current intersection point is in a triangle edge e and e is not the current search facet then 4 next search facet ← e 5 else if current intersection point is in a face side or face corner then 6 next search facet ← a cell face in which the side or corner is included, and which has not already been searched 7 else / current intersection point is in a cell face f 9 next search facet ← f Algorithm 7 next vertex = GetNextVertexInPolygonOnFaceSide(current intersection point, face side, search direction) Data: The current vertex in a face polygon under construction, the face side in which to find the next vertex in the face polygon, and a search direction for the face side.Result: The next vertex in the face polygon. 1 s ← relative position of the current intersection point on the face side 2 d ← a value larger than 1.0 3 for each intersection point p in the face side do 4 s p ← relative position of p on the face side 5 if p on correct side of the current intersection point given search direction then next search facet = ChangeSearchFacet(current intersection point, etAndUpdateSearchDirection(current intersection point, current search facet, current search direction) Data: The current vertex in a cell polygon under construction, the search facet in which the current vertex was found, and the current search direction.Result: The search facet in which to find the next vertex in the polygon, and the next search direction which needs update if the next search facet is a face side. 1 finished ← false 2 next search direction ← current search direction 3 if current intersection point is in a triangle edge or a triangle vertex then there is a a triangle t in which the current intersection point is included, and t has not already been searched then / current intersection point is in the interior of a triangle t 9 if current search facet is a face side then 10 next search facet ← t 11 finished ← true 12 if ¬ finished then 13 // next search facet is a face side, so we have to update the search direction 14 if current intersection point is in a face side s then current search facet is a triangle with normal n then 17 if the projection n s of n on s is not zero then 18 next search direction ← the direction corresponding to −n s 19 else 20 find a triangle in which the current intersection point is included, and which normal has a nonzero projection n s on s 21 next search direction ← the direction corresponding to −n s 22 else 23 next search direction ← current search direction 24 else / current intersection point is in a face corner c 26 if current search facet is a face side then 27 next search facet ← the face side in which c is included, and which is not the current search facet 28 next search direction ← current search direction 29 else 30 // current intersection point is in a triangle t, find face side to search in and search direction 31 next search facet ← the face side in which c is included, and which is located on the inside of t assuming that the normal of t points outwards 32 next search direction ← the direction of the next search facet pointing away from c Algorithm 9 {next search facet, next search direction} = ChangeSearchFacsection point, face side, search direction) Data: The current vertex in a face polygon under construction, the face side on which to find the next vertex in the face polygon, and a search direction for the face side.Result: The next vertex in the face polygon. 1 s ← relative position of the current intersection point on the face side 2 d ← a value larger than 1.0 3 for each intersection point p in the face side do 4 s p ← relative position of p on the face side 5 if p on correct side of the current intersection point given search direction then 6 if (p − current intersection point) • normal of triangle including p > 0.0 then 7 d p ← |s p − s| 8 if d p < d then 9 d ← d p 10 next vertex ← p

Table 1
Overview of algorithms in the exact method

Table 2
Overview of algorithms in the approximate method