Face-based Volume-of-Fluid interface positioning in arbitrary polyhedra
FFace-based Volume-of-Fluid interface positioning in arbitrary polyhedra
Johannes Kromer and Dieter Bothe † Mathematical Modeling and Analysis, Technische Universit¨at DarmstadtAlarich-Weiss-Strasse 10, 64287 Darmstadt, Germany † Email for correspondence: [email protected]
Abstract
We introduce a fast and robust algorithm for finding a plane Γ with given normal n Γ , which truncates anarbitrary polyhedron P such that the remaining sub-polyhedron admits a given volume α ∗ |P| . In the literature, thisis commonly referred to as Volume-of-Fluid (VoF) interface positioning problem . The novelty of our work istwofold: firstly, by recursive application of the
Gaussian divergence theorem, the volume of a truncated polyhedroncan be computed at high efficiency, based on summation over quantities associated to the faces of the polyhedron.One obtains a very convenient piecewise parametrization (within so-called brackets ) in terms of the signed distance s to the plane Γ. As an implication, one can restrain from the costly necessity to establish topological connectivity,rendering the present approach most suitable for the application to unstructured computational meshes. Secondly,in the vicinity of the truncation position s , the volume can be expressed exactly, i.e. in terms of a cubic polynomialof the normal distance to the PLIC plane. The local knowledge of derivatives enables to construct a root-findingalgorithm that pairs bracketing and higher-order approximation. The performance is assessed by conducting anextensive set of numerical experiments, considering convex and non-convex polyhedra of genus (i.e., number ofholes) zero and one in combination with carefully selected volume fractions α ∗ (including α ∗ ≈ α ∗ ≈
1) andnormal orientations n Γ . For all configurations we obtain a significant reduction of the number of (computationallycostly) truncations required for the positioning: on average, our algorithm requires between one and two polyhedrontruncations to find the position of the plane Γ, outperforming existing methods. Keywords—
PLIC, interface positioning, Volume-of-Fluid method
In geometric Volume-of-Fluid (VoF) methods, the approximative reconstruction of the interface from volume fractionsconstitutes a central element, in terms of both accuracy of the method and consumption of computational ressources. Inits most extensive scope, the reconstruction of an approximative interface includes the computation of the normal fieldfrom the volume fractions, followed by positioning the planar interface patches such that they match the respectivevolume fraction in each cell. The present work focusses on the second component, assuming the discrete field ofapproximative normals to be given. As we shall see below, positioning the interface translates to finding the uniqueroot of a scalar monotonous function, namely the volume of the truncated computational cell, which is parametrizedby the signed distance s . However, the truncated volume is only known implicitly, in the sense that its evaluationat some s requires a (computationally costly) truncation of the polyhedron. Therefore, there are two conceptuallyindependent components: Volume computation of a truncated polyhedron:
The vast majority of literature contributions, in one formor another, resorts to schemes involving truncation, (i) establishing topological connectivity and (ii) volumecomputation. The combinatorial complexity of task (i) depends on the topological properties of the originalpolyhedron as well as on the possible truncation outcomes. In the present context, computing the volume of1 a r X i v : . [ c s . C G ] J a n n arbitrary polyhedron by application of the Gaussian divergence theorem can hardly be outplayed in termsof efficiency; cf. Lopez et al. [15]. While explicit formulae can be given for primitive polyhedra (e.g., Yang andJames [25], Gueyffier et al. [10]), tetrahedral decomposition is applied in general , being both costly and limitedto topologically simple or convex polyhedra in most cases (e.g., Skarysz et al. [21] and tables 3 and 4 in Lopezet al. [15]). Root-finding:
The choice of an appropriate method depends, among others, on the available information: while thesecant method gets along with function evaluations,
Newtons method or consecutive cubic spline interpolation(e.g., Maric [16]) additionally require derivatives. For a thorough overview, the reader is referred to the booksof Stoer and Bulirsch [22] and Atkinson [2].The fact that the interface needs to be reconstructed within each computational cell and time step thus highlightsthe necessity for both general applicability (in terms of cell topology) and computational efficiency. The applicationof the Gaussian divergence theorem allows to address both of the aforementioned issues. As stated above, this is acommon concept in the literature. However, the present approach introduces the exploitation of a previously neglectedscope for design: pairing recursive application of a reduction in dimension via the
Gaussian divergence theoremwith proper translations of the coordinate system avoids the costly and, thus, disadvantageous necessity to extracttopological information from the truncated polyhedron. While this property considerably reduces the complexity forconvex polyhedra already, its full potential develops for non-convex polyhedra or those with a large number of faces.Furthermore, our method produces the derivatives with respect to s of the volume function at negligible additionalcosts, enabling the application of a root-finding algortihm based on locally quadratic approximation in combinationwith implicit bracketing .The remainder of this paper is organized as follows: after reviewing the relevant literature in subsection 1.1, theformal notation and data structure are summarized in subsection 1.3. Section 2 provides the mathematical detailsof the volume computation, along with a novel rigorous analytical treatment of the parametrized volume function.Since this paper introduces a numerical method meant to be integrated into existing codes, section 3 offers a detailedflowchart, which is further detailed throughout this work. In section 4, we present numerical results a set of convexand non-convex polyhedra. Graphics
The figures and flowcharts in this manuscript are produced using the versatile and powerful library pstricks . Forimplementation details as well as an instructive collection of examples, the reader is referred to the book of Voß [24]and https://tug.org/PSTricks/main.cgi , respectively.
Note 1.1 (Problem formulation) . In a VoF context, the phases are encoded by the discrete volume fraction field α ,from which, in geometric methods, at some point the numerical interface has to be reconstructed. The aforementionedreconstruction involves (i) the computation of the normal field from the volume fractions, followed by (ii) a piece-wise(i.e. cell-based) positioning of the disconnected numerical interface patches, ensuring that the volume fractions α inthe respective cells are conserved to a prescribed tolerance. The present work focusses on the latter step. Hence, theproblem under consideration, also referred to as volume-matching (Diot and Fran¸cois [8]) or volume conservationenforcement (VCE) (Lopez et al. [15]) in the literature, can be cast as follows: given Figure 6 in Ahn and Shashkov [1] contains a nice illustration of the implied complexity. Strang [23] introduces a directionally split geometric transport, requiring three reconstructions per time step. The term bracketing refers to finding the bracket containing the sought position s ∗ . a polyhedron P , • a unit normal n Γ and • a volume ≤ V ∗ ≤ vol( P ) or, alternatively, a volume fraction ≤ α ∗ ≤ with V ∗ = α ∗ vol( P ) ,find some s ∗ ∈ R such that the volume of the intersection of the polyhedron P with the negative halfspace of the PLICplane Γ , i.e. { x ∈ R : (cid:104) x , n Γ (cid:105) − s ≤ } , equals α ∗ times the polyhedron volume. More formally, one may writefind s ∗ ∈ R such that (cid:12)(cid:12)(cid:12)(cid:8) x ∈ P : (cid:104) x , n Γ (cid:105) − s ∗ ≤ (cid:9)(cid:12)(cid:12)(cid:12) = V ∗ . Irrespective of the spatial dimension of the polyhedron P , the interface positioning problem is one-dimensional. Thisproperty reflects the fact that the (signed) distance of two planes with a common normal n Γ , say, Γ( s i ) and Γ( s j ) , canbe computed from two arbitrary points on the respective planes, i.e. (cid:104) x − y , n Γ (cid:105) = s i − s j for all ( x , y ) ∈ Γ( s i ) × Γ( s j ) . Rider and Kothe [18] find the position of the PLIC plane iteratively by application of the method of Brent [3]. Ananalytical relation of the volume of truncated cuboids is introduced in Scardovelli and Zaleski [19], resorting to apiecewise definition within the respective so-called brackets . A similar set of formulae was developed earlier byZemach [26] and implemented in Kothe et al. [11]. Ahn and Shashkov [1] augmented the work of Rider and Kothe [18]by secant and bisection methods, also extending the applicability to general polyhedral meshes. However, Ahn andShashkov [1] resort to a tetrahedral decomposition of the original polyhedron. Based on a small number of differentintersection configurations, Yang and James [25] and Gueyffier et al. [10] derive analytical expressions for the volumeof truncated tetrahedra and cuboids, respectively. Diot and Fran¸cois [8], extending the work of Diot et al. [7] from twoto three spatial dimensions, propose an analytical approach to obtain the PLIC position in general convex polyhedra.After aligning the normal with the z -axis by applying a rotation of the coordinate system, i.e. n Γ (cid:55)→ e z , the polyhedronis truncated at its vertices ν i by so-called s i -planes with s i := (cid:104) ν i , n Γ (cid:105) and s i ≤ s i + i . The intersection of P withthe negative halfspace of each of these s i -planes induces a volume V s,i , where V s,i ≤ V s,i +1 . Since 0 ≤ α ∗ ≤
1, thereis exactly one bracket B i := ( s k , s k +1 ] such that s ∗ ∈ ( s k , s k +1 ] and V s,k < V ∗ ≤ V s,k +1 . The corresponding planestruncate a prismatoid from the polyhedron; cf. figure 1 for an illustration. A polynomial of degree three in s , whose n Γ r n Γ r Figure 1: Rotation of cuboid s.t. n Γ → e z ; center: planes at s k / s k +1 (red, passing through cell vertices (cid:4) ) bracketingthe sought plane Γ (blue); right: prismatoid (obtained from two consecutive planes truncating the cuboid) whosevolume is computed analytically; cf. Diot and Fran¸cois [8].coefficients Diot and Fran¸cois [8] provide analytically, governs the enclosed volume, implying that the sought s ∗ canbe computed directly by means of root-finding. In order to obtain analytical expressions for the prismatoid volume,Diot and Fran¸cois [8] require considerable connectivity information of the underlying polyhedron a priori, along withestablishing the connectivity of the intersections and cell vertices for each truncation operation. Despite the ability Section 2 below contains a formal definition. Diot and Fran¸cois [8] employ the notation d i , which we have adapted to s i for reasons of notational consistency.
3o considerably reduce the runtime by application of a tailored data structure, establishing connectivity increases thecomputational effort for complex polyhedra. For convex polyhedra, Skarysz et al. [21] propose an algorithm combininga tetrahedral decomposition with three root finding methods: bisection as well as the algorithms of Muller [17] andBrent [3], where the latter, by design, makes no use of derivatives. By means of topological considerations of thetwo possible intersection configurations, the computation of the volumes of truncated tetrahedra can be carried outefficiently. A secondary, yet very important finding of Skarysz et al. [21] is the fact that the number of tetrahedraresulting from the decomposition crucially affects the overall algorithm efficiency (for the cases reported, the averagenumber of volume computations is 6 − Gaussian divergence theorem to express the volume of the truncated polyhedron as aweighted sum of the areas of its bounding faces, the latter being computed as a sum of outer vector products. Thevertices of the truncated faces stem from both the original polyhedron and the edge intersections. For the applicationof the formula of Goldman [9], they require to be arranged counterclock-wise with respect to the face normal. Lopezand Hernandez [12] extract the aforementioned arrangement from the polyhedron topology, which quickly becomescumbersome for polyhedra with a large number of faces. The position s ∗ is obtained by direct computation of theroot of a cubic polynomial within a bracket B i (cf. section 2), which requires sequentially intersecting the polyhedronwith various planes parallel to Γ passing through its vertices. Despite a carefully designed selection of the vertices,Lopez and Hernandez [12] state that ” this part of the algorithm may be the most time consuming, especially when thenumber of vertices of P is large “. Since Lopez and Hernandez [12] do not report the average number of truncations,a meaningful comparison to the present work cannot be made. Lopez et al. [14] develop a set of geometrical tools forconvex polyhedra, where the PLIC positioning resorts to an algorithm of Lopez and Hernandez [12]. In essence, acoupling of bracketing and analytical description of the volume function ( C oupled I nterpolation- Br acketed A nalytical V olume E nforcement, CIBRAVE for short) allows to obtain the sought s ∗ from spline interpolation. Each iteration s n , obtained from linear interpolation of { α left , α right } , commences by computing the index i of the parenting bracket(such that s i ≤ s n ≤ s i +1 ). Then, a truncation of the polyhedron at s i yields the four coefficents of the cubic spline S i ,allowing to evaluate α ( s i ) = S i ( s i ) and α ( s i +1 ) = S i ( s i +1 ). For α ( s i ) ≤ α ∗ ≤ α ( s i +1 ), the sought s ∗ corresponds tothe root of S i ( s ) − α ∗ . Otherwise, the volume fractions at the bracket boundary and { α left , α right } are used to update s n by linear interpolation; cf. note 2.4 and figure 2 for an illustration. Lopez et al. [14] provide numerical results fortetra-, hexa-, icosa- and dodecahedra, to which we compare this work within section 4. By introducing a conceptuallyextensive method to extract connectivtiy from topological information, Lopez et al. [13, 15] extend the applicability ofthe previously mentioned algorithms to non-convex polyhedra, thereby nicely illustrating the complexity of establishingvertex connectivity after truncation. Remark 1.1 (Truncation of non-convex polyhedra (Lopez et al. [15])) . The figure below resembles an adaption offig. 8 from Lopez et al. [15], illustrating the capability of handling non-simply connected faces. However, the properarrangement of the polyhedron vertices and intersections, indicated by the sets on the right, requires (a) connectivityon two hierarchical levels (vertex to face/face to face) and (b) interlacing with the intersections induced by the plane Γ . Lopez and Hernandez [12] employ the notation C Γ c , which we have adapted to s ∗ for reasons of notational consistency. sα ∗ l i n e a r i n t e r p . y i e l d s s n s n u r α left := α ( s i ) s i s i +1 α − α ∗ sα ∗ l i n e a r i n t e r p . y i e l d s s n + r α right := α ( s i +1 ) s i s i +1 u s n +1 sα ∗ s i s i +1 e x a c t s p l i n e y i e l d s s ∗ b s ∗ Figure 2: CIBRAVE method proposed by Lopez et al. [14]: the first iteration s n (left, (cid:78) ) lies not within the targetbracket B ∗ = [ s i , s i +1 ]. Due to α ( s n ) < α ∗ , α left is updated, such that a linear interpolation yields s n +1 (center, (cid:78) ).Since α ( s i ) ≤ α ( s n +1 ) ≤ α ( s i +1 ), the polyhedron is truncated at s i . Finally, due to α ( s i ) < α ∗ < α ( s i +1 ), a cubicspline interpolation yields the sought s ∗ (right, • ), resulting in two polyhedron truncations. n Γ bbb bbbb b
12 3 4 567 89 1011 1213141516 1718 { ν , ν , ν , ν } → { ν , ν , ν , ν }{ ν , ν , ν , ν } → { ν , ν , ν , ν }{ ν , ν , ν , ν } → { ν , ν , ν , ν }{ ν , ν , ν , ν } → { ν , ν , ν , ν }{ ν , ν , ν , ν , ν } → { ν , ν , ν } , { ν , ν , ν }{ ν , ν , ν , ν , ν } → { ν , ν , ν } , { ν , ν , ν }{ ν , ν , ν , ν }{ ν , ν , ν , ν }
12 1114 12139 1516 10 1817
Even for the comparatively simple case illustrated above, figs. 9-12 in Lopez et al. [15] indicate the computational effortbehind establishing the connectivity.
Unfortunately, Lopez et al. [15] only provide a performance assessment of their algorithm in terms of relative CPUtimes, but not in terms of truncations, stating in their fig. 16 that the evaluation of the spline coefficients accounts forroughly 80% of the total runtime. As a secondary finding, table 3 in Lopez et al. [15, p. 684] indicates the superiorityof the application of
Gaussian divergence theorem over convex decomposition. Investigating an extended versionof the problem under consideration here, Chen and Zhang [4] derive analytical formulae for the derivative of theparametrized volume function. The authors construct a ” modified Newton’s method which utilizes the monotonicity of V ( s ) [notation adapted] to ensure the convergence of the iteration “. In comparison to the bisection/secant method ofAhn and Shashkov [1], Chen and Zhang [4] report a reduction of iterations by 60 to 66%. Dai and Tong [6], followingan idea of Diot and Fran¸cois [8], apply a cell rotation to align the PLIC normal with e along with a decompositionof the convex cell into N − N is the number of unique signed distances of the rotated vertices.The volume of the respective sub-cells is computed by exploiting the Gaussian divergence theorem, which allows toidentify the subcell within which the sought PLIC plane must be located to obtain the prescribed volume fraction. Theenclosed volume is expressed as a cubic polynomial of the normalized distance of the bounding planes, whose root iscomputed by a
Newton-Raphson algorithm. The cell types under consideration in Dai and Tong [6] include variouspolyhedra with 12–48 vertices and 8–26 faces; e.g., hexagonal prisms, regular icosahedra and truncated cuboctrahedra(sic!); see Dai and Tong [6, tab. 1]. The instances of the normal n Γ and volume fractions α for the numericalexperiments are arranged in one and two sets, including the potentially troublesome cases where the PLIC plane isaligned with one of the cell faces and/or the cell is almost full or empty; cf. subsection 1.2. The authors report very5ccurate reconstruction results, with both L - and L ∞ -type deviations of magnitude 10 − . Unfortunately, Dai andTong [6] do not provide average iteration numbers to which the results in section 4 could be compared. Recently, Maric[16] proposed a consecutive cubic spline (CCS) interpolation, which approximates the sought root s ∗ by iterativelyconstricting s left < s ∗ < s right such that [ s left , s right ] maintains a sign change; cf. figure 3. Starting from s left := s − and s right := s + , each iteration starts by updating[ s left , s right ] := [ s n , s right ] D ( s n ) < s left , s n ] D ( s n ) > D ( s ) := α ( s ) − α ∗ . (1)Figure 3 illustrates the update s n +1 as a function of the sign of D ( s n ) D ( s n − ), where, without loss of generality, itcan be assumed that s n > s n − . If the current ( s n , (cid:78) ) and previous ( s n − , (cid:4) ) intersection produce deviations D withi) different signs (center), they comprise the sought root s ∗ , i.e. s ∗ ∈ [ s n − , s n ], such that a cubic spline interpo-lation based on D and its derivative at s n and s n − , respectively, provides the next iteration ( s n +1 , • ).ii) coinciding signs (left), the sought root s ∗ lies outside of [ s n − , s n ]. In this case, the next iteration ( s n +1 , • ) iscalculated from the root of a linear approximation of D around the current iteration s n .It is worth noting that, if s left and s right are located in the same bracket, i.e. [ s left , s right ] ⊂ [ s i , s i +1 ] (rightmost infigure 3), the cubic spline exactly interpolates the deviation D . The numerical results of Maric [16] highlight that thisis especially advantageous for α ∗ ≈ α ∗ ≈
1, rendering this approach superior to
Newton -based ones, whoseconvergence strongly suffers from diminishing derivatives.0 0 . . s − α ∗ − α ∗ li n . a p p r o x . ur D ( s n − ) < D ( s n ) < b s n +1 DD . . s − α ∗ − α ∗ c u b i c s p l i n e ru D ( s n − ) > D ( s n ) < b s n +1 DD s − α ∗ − α ∗ s i s i +1 s n − s n e x a c t s p l i n e r u Figure 3: CCS method proposed by Maric [16]: computation of next iteration ( s n +1 , • ) for coinciding (left) anddifferent (center) sign of deviation D at current ( s n , (cid:78) ) and previous iteration ( s n − , (cid:4) ). If the interval [ s left , s right ]is confined by a bracket, the spline interpolates D exactly (right).Beyond insightful comments on the involved floating point arithmetics, Maric [16] assesses the computationalperformance of his algorithm in detail, concluding that, in comparison to the polyhedron truncation, the computationalcosts of rootfinding are negligible. Section 4 contains a detailed comparison of the results. As a consequence of the literature review, we formulate two main objectives for the present work: O1 . Irrespective of the conceptual differences of the existing approaches, the main bottleneck for most of the publishedalgorithms is an efficient and conceptually simple computation of the volume of truncated polyhedra. In the vastmajority of the cited publications, extracting connectivity from topological information is responsible for most ofthe computational effort. Thus, the first objective of this work is to establish an algorithm for the computation ofthe volume of a truncated arbitrary polyhedron, which restrains from extracting topological information s much as possible . Note 1.3 further discusses the computational effort of extracting toplogical information. O2 . As has been pointed out by Skarysz et al. [21], Dai and Tong [6] and Maric [16], the performance of a recon-struction algorithm crucially depends on the input parameters { α ∗ , P , n Γ } . Hence, the numerical experimentsconducted to assess the performance shall resort to a dyadic combination of representative sets of geome-tries P , volume fractions α ∗ and normal orientations n Γ . While table 1 and appendix A gather the polyhedraunder consideration, the volume fractions and normals are adapted from Maric [16, eqs. (25) to (28)], i.e. S n Γ ( M n Γ ) := (cid:40) [cos ϕ sin θ, sin ϕ sin θ, cos θ ] T : ( ϕ, θ ) ∈ π M n Γ [1 , , . . . , M n Γ ] × πM n Γ [0 , , . . . , M n Γ ] (cid:41) , S α ∗ ( M α ∗ ) := { − k : 5 ≤ k ≤ } ∪ (cid:26) − + m − M α ∗ − (cid:16) − · − (cid:17) : 1 ≤ m ≤ M α ∗ (cid:27) ∪ { − − k : 5 ≤ k ≤ } . We set M n Γ = 80 and M α ∗ = 50, which corresponds to 2 M n Γ ( M n Γ +1)( M α ∗ +10) = 777 600 numerical instancesper polyhedron. Note that S n Γ corresponds to a discrete approximation of the boundary of the unit sphere,denoted S := { x ∈ R : (cid:107) x (cid:107) = 1 } . Taking into account volume fractions close to zero or one becomes relevantin cases where numerical whisps occur. Here, a reconstruction tolerance of 10 − is chosen. Geometric quantities
We consider an arbitrary polyhedron P bounded by N F planar polygonal faces F k withouter unit normal n F ,k . To ensure applicability of the Gaussian divergence theorem, ∂ F k and ∂ P = (cid:83) k F k arerespectively assumed to admit no self-intersection . For the desired application within finite-volume based methods,however, self-intersecting polyhedra have no relevance. The N ν k vertices { x F k,m } on each face are ordered counter-clockwise with respect to the normal n F ,k , implying that the m -th edge E k,m is spanned by ( x F k,m , x F k,m +1 ). Fornotational convenience, assume that the indices are continued periodically, i.e. let x F k,N ν k +1 = x F k, . To each edge E k,m on face F k we assign the outer unit normal N k,m such that (cid:104) N k,m , n F ,k (cid:105) = (cid:104) N k,m , x F k,m +1 − x F k,m (cid:105) = 0. Figure 4illustrates the setup and its relevant quantities. n Γ x Γ (a) intersection with PLIC plane Γ n F k N k, N k, N k, N k, b bbb x F k, x F k, x F k, x F k, (b) face F k with normal n F k and co-normals N k,m F − k E − k, E − k, (c) interior edges E − k,m Figure 4: Illustration of the relevant quantities for a cuboid. Two edges of a polygon are considered intersected if they have a common point which is not one of the respective vertices. Theanalogous holds for faces. Note that the number of vertices in a closed polygon coincides with the number of edges. ummation Throughout this manuscript, summation over faces F k and associated edges E k,m employ the subscripts k and m , where the respective limits always are k = 1 . . . N F ( m = 1 . . . N ν k ). For ease of notation, the summationlimits are omitted. Representation of polyhedra
Table 1 gathers the polyhedra employed for the numerical experiments in section 4.Table 1: Polyhedra considered for the numerical investigation within this work; cf. note 1.2. name parameters unit cuboid(cube) – 6 12 8 edge length 1(iso)dodecahedron – 12 30 20 edge length 1torus R , γ , N , N N N N N N N cf. appendix A.2letter A – 15 33 22 cf. appendix A.2unit tetrahedron – 4 6 4 – Note 1.2 (Parametrized polyhedra) . In principle, the classes of polyhedra given in table 1 admit further parametriza-tion: e.g., the normalized edge lengths of a general cuboid are { , ψ , ψ } with ( ψ , ψ ) ⊂ R . The implications ofconsidering the parametrization becomes apparent for the tetrahedron: without loss of generality, assume that one ofthe vertices coincides with the origin and one of the faces (containing the aforementioned vertex) lies in the xy -plane.Then, any instance {P , α ∗ , n Γ } with an arbitrary (non-degenerate) tetrahedron P can be transformed to the unit tetra-hedron by a linear transformation Q : R (cid:55)→ R . While such a transformation preserves the volume fraction α ∗ , thenormal experiences a non-linear transformation n Γ (cid:55)→ T ( n Γ ) with T ( x ) := Qx √ (cid:104) Qx , Qx (cid:105) . QQ − bb bbbb b bbb (1 , ,
0) (0 , , , , bb bbbb b bbb ( a, ,
0) ( b, c, d, e, f ) Note that the transformation T maps the boundary of the unit sphere onto itself, i.e. S = T ( S ) . Thus, if the set ofnormals S n Γ is dense in S , it is sufficient to numerically investigate the unit tetrahedron. Since S n Γ resembles a finiteapproximation of the boundary of the unit sphere, this is not the case in practice; cf. subsection 1.2. Furthermore,even moderate deviations from the unit tetrahderon induce strong distortions to the sample space for the normals S n Γ .Thus, the appropriate choice of S n Γ for general transformations Q poses a complex task in itself, which, however, goesbeyond the scope of this work. ϕθ ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ϕθ ××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××× × × × × × × ×××××××××××××××××××××××××××××××××××××××××××× × × × × × × × × ××××××××××××××××××××××××××××××××××××××××××× × × × × × × × × × ××××××××××××××××××××××××××××××××××××××××××× × × × × × × × × × ×××××××××××××××××××××××××××××××××××××××××××× × × × × × × × ××××××××××××××××××××××××××××××××××××××××××××××× × × × × × ××××××××××××××××××××××××××××××××××××××××××××××××× × × ××××××××××××××××××××××××××××××××××××××××××××××××××× × ××××××××××××××××××××××××××××××××××××××××××××××××××× × × ××× ××××××××××××××××××××××××××××××××××××××××××××××× × × ×× × ×××××××××××××××××××××××××××××××××××××××××××××× × × ×× × × × ××××××××××××××××××××××××××××××××××××××××××××× × ×× × × × × ××××××××××××××××××××××××××××××××××××××××××××× ×× × × × × × ×××××××××××××××××××××××××××××××××××××××××××× ×× × × × × × × ××××××××××××××××××××××××××××××××××××××××××××× × × × × × × × ××××××××××××××××××××××××××××××××××××××××××××× × × × × × × × ××××××××××××××××××××××××××××××××××××××××××××× × × × × × × × ××××××××××××××××××××××××××××××××××××××××××××× × × × × × × × ×××××××××××××××××××××××××××××××××××××××××××× × × × × × × × × × ××××××××××××××××××××××××××××××××××××××××××× × × × × × × × × × ××××××××××××××××××××××××××××××××××××××××××× × × × × × × × × × ××××××××××××××××××××××××××××××××××××××××××× × × × × × × × × ××××××××××××××××××××××××××××××××××××××××××××× × × × × × × ×××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××××× ×××××××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××××××××××××××××××××××××××××××××××× ××××××××××××××××××××××××× Distortion of uniformly discretized S for { / , / , / , − / , − / , / } , corresponding to the figure above.
8n algorithmic terms, the following two data sets represent a polyhedron:1. a list of N ν common vertices { ν i } . ( In theory, every face F k could hold its own set of vertices. )2. a list of N F faces {F k } , each defined by a periodic set of indices { i m } with 1 ≤ m ≤ N ν k . The correspondingvertices form a closed planar polygon, ordered counter-clockwise with respect to the face normal n F ,k . ( Here,by design, a hole in a face is realized by an additional coplanar face with inverted orientation; cf. figure 5. ) b bbbb bb bbbb
14 3256 78 { , , , , , , , , , } b bbbb bb bbb
14 3256 78 { , , , } { , , , } Figure 5: Different strategies to define a closed polygon with a hole: single polygon with overlapping edges ( { , } and { , } ) and decomposition with appropriate ordering. Note 1.3 (Computational effort vs. topology) . The above data sets comprise the minimum information required tocompute the volume of a polyhedron. In fact, most unstructured meshes employed in VOF-codes resort to this typeof description. In theory, one could deduce a list of edges, each respectively shared by two faces. However, thisconnectivity has to be either (i) provided a priori (increasing memory consumption) or (ii) established at runtime(increasing CPU consumption), where dynamically changing meshes cannot resort to (i). At this point, the authorswould like to emphasize that, obviously, the computational effort of volume computation could be reduced roughly byhalf, since neglecting the edge connectivity implies that the number of edge operations (intersections) doubles with thenumber of faces N F . On the other hand, the complexity of establishing the edge-connectivity between faces increases quadratically with N F , i.e., it becomes relatively expensive, even for a moderate number of faces N F . Thus, wepursue a face-based design, exclusively resorting to intra-facial connectivity. For a given polyhedron P and a unit normal n Γ , we wish to compute a base point x ∗ Γ for a plane Γ( x Γ ) with normal n Γ , such that the intersection of the negative halfspace of Γ with P , denoted P − ( x Γ ) := { x ∈ P : (cid:104) x − x Γ , n Γ (cid:105) ≤ } , (2)admits some given positive volume V ∗ ≤ |P| , where henceforth V ( x Γ ) := |P − ( x Γ ) | denotes the volume of the truncatedpolyhedron. Note that V ∗ = V ( x ∗ Γ ). With an arbitrary (but spatially fixed) reference x Γ , , the signed distance s isused to parametrize the base point in eq. (2) as x Γ = x Γ , + s n Γ , implyingΓ( x Γ ) = Γ( s ; x Γ , ) = { x ∈ R : φ Γ ( x ; s, x Γ , ) = 0 } with φ Γ ( x ; s, x Γ , ) := (cid:104) x − x Γ , , n Γ (cid:105) − s. (3)Then, we obtain the problem formulated in note 1.1, i.e.find s ∗ ∈ R such that (cid:12)(cid:12)(cid:12)(cid:8) x ∈ P : (cid:104) x − x Γ , , n Γ (cid:105) − s ∗ ≤ (cid:9)(cid:12)(cid:12)(cid:12) = V ∗ . (4)9n what follows, a prime denotes the derivative with respect to the signed distance s . Non-dimensional formulation.
Throughout this work, we shall resort to the normalized volume or volume fraction α ( s ) := V ( s ) |P| − . The signed distances of the verticesˆ S := (cid:91) k (cid:110) (cid:104) x F k,m − x Γ , , n Γ (cid:105) (cid:111) (5)correspond to the values of s for which the plane Γ passes through a vertex (single or multiple) of the polyhedron P .The elements s i of ˆ S , being unique ( s i (cid:54) = s j for i (cid:54) = j ) and arranged in ascending order ( s i < s i +1 ), define mutuallydisjoint so-called brackets B i := ( s i , s i +1 ] with B := [ s , s ]. The target bracket B ∗ contains the sought root s ∗ ;cf. figure 6 for an illustration in two spatial dimensions. Furthermore, let s s =: s + s s s =: s − x Γ , n Γ s b Figure 6: Polygonal face with signed distances of the vertices.˜ s := s − s − s + − s − with s − := min i s i , s + := max i s i and s ∆ := s + − s − , (6)where, henceforth, the tilde denotes the normalized pendant of a variable. Note that ∂ ˜ s = s ∆ ∂ s . Regularity.
Globally, V ( s ) ∈ C is a strictly increasing function, whose derivative corresponds to the area of thePLIC plane contained in the computational cell , i.e. V (cid:48) ( s ) = |P ∩ Γ( s ) | ≥
0. Locally, i.e. within a bracket B i , thevolume is a third-order polynomial in s . Under mild restrictions concerning the PLIC normal n Γ , cf. note 2.1, thetruncated volume is globally continuously differentiable, i.e. V ( s ) ∈ C . Note 2.1 (Faces parallel to Γ) . Consider a face F k which is parallel to the PLIC plane Γ , i.e. ± (cid:104) n F k , n Γ (cid:105) =: σ k .Hence, depending on the value of s , the face is either entirely contained in the PLIC plane or not at all, corresponding toa multivalued (in the sense of non-smooth analysis; see, e.g., Clarke [5]) derivative α (cid:48) ( s ) for the respective parameters. This identity directly results from application of the
Reynolds transport theorem to eq. (11). However, Chen and Zhang [4, eq. (30)]offer an intuitive geometric derivation. x σ = − σ = σ = − n Γ α = 0 α ′ ∈ [0 , case i n Γ α = / α ′ ∈ [ / , case ii
01 ˜ s c a s e i c a s e ii s − s + αα ′ Let H P denote the convex hull of the polyhedron and ∂ H P its boundary, with ∂ H P = ∂ P = (cid:83) k F k for convexpolyhedra. The regularity of α ( s ) depends on the topological properties of the respective face F k :i. If F k is an element of ∂ H P ( in the leftmost figure), the corresponding signed distance is s − ( σ k = − ) or s + ( σ k = 1 ). These boundary-discontinuities of α (cid:48) ( s ) can be resolved by assigning the respectively interior limit(right-sided for σ k = − and left-sided for σ k = 1 ). Thus, let α (cid:48) ( s − ) := lim s (cid:38) s − α (cid:48) ( s ) = |P| − (cid:88) k |F k ∩ Γ( s − ) | and α (cid:48) ( s + ) := lim s (cid:37) s + α (cid:48) ( s ) = |P| − (cid:88) k |F k ∩ Γ( s + ) | , (7) where the sum accounts for multiple coplanar faces.ii. A face F k which is not an element of ∂ H P ( in the leftmost figure) induces an interior jump discontinuity( ◦ in the rightmost figure) in α (cid:48) ( s ) , where α (cid:48)(cid:48) ( s ) may jump as well if Γ simultaneously passes through a vertexnot contained in F k . The third figure above illustrates the multi-valued derivative as α (cid:48) ∈ [ / , at α = / .Moreover, the absence of faces parallel to Γ implies vanishing volume derivatives at the boundaries, i.e. { σ k } (cid:54)(cid:51) ± ⇒ α (cid:48) ( s ± ) = 0 , (8) which can be assumed to be the case in the vast majority of applications. Hence, any convex polyhedron always admitsa globally continuously differentiable volume function, i.e. it holds that α ( s ) ∈ C for all s − < s < s + , irrespective ofthe alignment of n Γ with the normals of the faces. For non-convex polyhedra, the regularity depends on { σ k } . The desired numerical treatment of eq. (4) suggests a non-dimensional formulation. Thus, we define the volumefraction deviation D ( s ) := α ( s ) − α ∗ with D : [ s − , s + ] (cid:55)→ [ − α ∗ , − α ∗ ] . (9)Since the volume fraction α is (at least) continuous and strictly increasing, s ∗ corresponds to the unique root of thedeviation D in [ s − , s + ], whose approximation we compute iteratively by s n +1 = s n + ∆( s n ) . (10)The computation of the initial value s and step size ∆ s n , respectively, are deferred to subsections 2.3 and 2.4. Theupcoming section derives the expressions for the truncated volume V ( s ) and its derivatives.11 .1 Volume computation By application of the
Gaussian divergence theorem, the volume of the truncated polyhedron P − can be written as V = (cid:90) P − x = 13 (cid:90) ∂ P − (cid:104) x − x , n ∂ P − (cid:105) d x = 13 (cid:90) ∂ P − \ Γ (cid:104) x − x , n ∂ P − (cid:105) d x + (cid:90) Γ (cid:104) x − x , n Γ (cid:105) d o , (11)where the choice of the reference x is arbitrary . The fact that x Γ ∈ Γ suggests to choose x := x Γ , + s n Γ , inwhich case the integral over Γ vanishes and eq. (11) simplifies to V ( s ) = 13 (cid:90) ∂ P − \ Γ (cid:104) x − x Γ , − s n Γ , n ∂ P − (cid:105) d x . (12) Remark 2.1 (Invariance under translation/rotation) . Any affine linear transformation T : x (cid:55)→ Q T x + x with det Q = 1 and Q T Q = I preserves the volume of a bounded domain Ω ⊂ R , implying that the volume of Ω isinvariant under translation and rotation. The integrand in the rightmost expression of eq. (11) corresponds to the localsigned distance (which is invariant under rotation) of the boundary to the translated origin x . Since the boundary ∂ P − is composed of planar segments, the signed distance is constant within each segment and, hence, can be extractedfrom the integral, leaving to compute the area of the respective segment F − k . If the signed distance to the origin is zero,the respective segment does not contribute to the volume integral. Thus, if the origin is translated to be coplanar to Γ (note that x ∈ P ∩ Γ is not required), one obtains (cid:82) Γ (cid:104) x − x , n Γ (cid:105) d o ≡ . Note that the integration domain in the above expression is piecewise planar, with, say, ∂ P − \ Γ =: (cid:83) k F − k , where F − k = F − k ( s ) denotes the segment of F k located in the negative halfspace of Γ( s ); cf. figure 4. Equation 12 becomes V ( s ) = 13 (cid:88) k (cid:104) x F k, − x Γ , − s n Γ , n F ,k (cid:105) (cid:90) F − k o = 13 (cid:88) k ( γ V ,k + sγ V ,k ) (cid:90) F − k o = 13 (cid:88) k ( γ V ,k + sγ V ,k ) A k ( s ) , (13)where the coefficients γ V ,k := (cid:104) x F k, − x Γ , , n F ,k (cid:105) and γ V ,k := −(cid:104) n F ,k , n Γ (cid:105) (14)depend only on static geometric quantities , which is advantageous in terms of computational efficiency. Note thateq. (13) corresponds to eq. (11), where the spatial dimension is reduced by one. Applying the Gauss theorem oncemore, the area A k := |F − k | of a planar face can be written as A k ( s ) = (cid:90) F − k o = 12 (cid:90) ∂ F − k \ Γ (cid:104) x − x ,k ( s ) , N ∂ F k (cid:105) d s + (cid:90) ∂ F − k ∩ Γ (cid:104) x − x ,k ( s ) , N ∂ F k (cid:105) d s . (15) Remark 2.2 (Area computation) . In one form or another, most literature contributions, e.g., Lopez and Hernandez[12], apply the formula of Goldman [9] to compute the area of a planar polygon: In the literature, the common choice is x ≡ ; see, e.g., Lopez and Hernandez [12], Lopez et al. [14, 15]. Static in the sense of not depending on the position of Γ, i.e., on s . + + +Thereto, its vertices – potentially containing intersections of ∂ F k with Γ – must be ordered counter-clockwise withrespect to n F ,k , imposing additional computational effort. On the other hand, the application of the Gaussian diver-gence theorem requires no connectivity beyond the a priori assignment of vertices to edges:= + + + + +
In order to eliminate the last summand in eq. (15) as above, one needs to find a point x ,k ∈ Γ with zero normaldistance to F k , i.e. such that (cid:104) x ,k − x F k, , n F k (cid:105) = 0. Note that the latter does not imply that x ,k is an element of F − k (nor, in fact, of F k ). If the PLIC plane Γ and the face F k are not parallel, i.e. for |(cid:104) n F k , n Γ (cid:105)| <
1, their line ofintersection (cf. figure 7 for an illustration) contains x ,k ( s ) = (cid:104) x Γ , , n Γ (cid:105) + s − (cid:104) x F k, , n F k (cid:105)(cid:104) n F k , n Γ (cid:105) − (cid:104) n F k , n Γ (cid:105) n Γ + (cid:104) x F k, , n F k (cid:105) − ( (cid:104) x Γ , , n Γ (cid:105) + s ) (cid:104) n F k , n Γ (cid:105) − (cid:104) n F k , n Γ (cid:105) n F k . (16)As a beneficial side-effect, this choice of the reference point, compared to, say, the origin, reduces the numerical errors;cf. Shewchuk [20]. Decomposing the boundary as ∂ F − k \ Γ =: (cid:83) m E − k,m and inserting the base point from eq. (16),eq. (15) becomes A k ( s ) = 12 (cid:88) m (cid:104) x F k,m − x ,k ( s ) , N k,m (cid:105) (cid:90) E − k,m s . (17)Introducing l k,m ( s ) := (cid:82) E − k,m s for notational convenience, one obtains A k ( s ) = 12 (cid:88) m (cid:16) γ A ,k,m + sγ A ,k,m (cid:17) l k,m ( s ) (18)with the coefficients γ A ,k,m := (cid:104) x F k,m , N k,m (cid:105) + (cid:104) x F k, , n F k (cid:105)(cid:104) n F k , n Γ (cid:105) − (cid:104) x Γ , , n Γ (cid:105) − (cid:104) n F k , n Γ (cid:105) (cid:104) N k,m , n Γ (cid:105) and γ A ,k,m := − (cid:104) N k,m , n Γ (cid:105) − (cid:104) n F k , n Γ (cid:105) . (19)As in eq. (14), the coefficients depend only on static geometric quantities, implying that their evaluation is requiredonly once per cell. Note that l k,m , corresponding to the length of the segment of the edge E k,m located within thenegative halfspace of Γ, is a positive piecewise linear function of s , implying that the derivative l (cid:48) k,m is a positiveconstant, once again depending only on static geometric quantities.Finally, combining eqs. (18) and (13) yields the truncated volume V ( s ) = 16 (cid:88) k ( γ V ,k + sγ V ,k ) (cid:88) m (cid:16) γ A ,k,m + sγ A ,k,m (cid:17) l k,m ( s ) , (20)which is a continuous function composed of third-order polynomials in s ; cf. note 2.1. Conceptually, the above equationcan be interpreted as follows: by recursively applying the Gaussian divergence theorem, the volume integral reduces toa set of line integrals (corresponding to the computation of the interior edge lengths l k,m ), while the static coefficients13 b bbrrb x , k x F k, x F k, x F k, x F k, i n t e r s e c t i o n o f Γ a n d F k u br x , k x F k, (cid:10) x F k , − x , k , N k , (cid:11) > N k, l k, = |E − k, | u rr x , k x Γ1 ,k (cid:10) x Γ , k − x , k , N k , Γ (cid:11) = N k, Γ Figure 7: Illustration of the coefficients for the face area computation given in eq. (19). Note that, due to the choiceof x ,k in eq. (16), the boundary segment located on Γ (right) does not contribute to the area integral. { γ V ,k , γ V ,k , γ A ,k,m , γ A ,k,m } contain the required topological information. The derivatives of eq. (20) are V (cid:48) ( s ) = 13 (cid:88) k ( γ V ,k + sγ V ,k ) A (cid:48) k ( s ) + γ V ,k A k ( s ) ,V (cid:48)(cid:48) ( s ) = 13 (cid:88) k ( γ V ,k + sγ V ,k ) A (cid:48)(cid:48) k + 2 γ V ,k A (cid:48) k ( s ) and V (cid:48)(cid:48)(cid:48) ( s ) = (cid:88) k γ V ,k A (cid:48)(cid:48) k (21)with the face area derivatives A (cid:48) k ( s ) = 12 (cid:88) m (cid:16) γ A ,k,m + sγ A ,k,m (cid:17) l (cid:48) k,m + γ A ,k,m l k,m ( s ) and A (cid:48)(cid:48) k = (cid:88) m γ A ,k,m l (cid:48) k,m . (22) Remark 2.3 (Comparison to Lopez et al. [14]) . In their eq. (17), Lopez et al. [14] obtain a representation of thevolume similar to eq. (20), say V i ( s ) = β ,i s + β ,i s + β ,i s + β ,i for s ∈ [ s i , s i +1 ] by truncating the polyhedron at either s i or s i +1 . However, the complex geometric operations involved imply a con-siderable increment in computational effort. Noting that the highest non-trivial derivative of eq. (20) is V (cid:48)(cid:48)(cid:48) ( s ) , beingpiecewise constant within a bracket [ s i , s i +1 ] , yields the polynomial coefficients as a by-product of a single truncationat an arbitrary s ∈ [ s i , s i +1 ] ; cf. subsection 2.4.1. The robust numerical computation of the lengths l k,m and, consecutively, the areas A k , along with their derivatives,requires considering the topological properties of the involved entities, which shall be addressed in the upcomingsubsection. The decomposition of a non-degenerately (note 2.2 considers the degenerate case) intersected face F k requires ahierarchically consistent evaluation of the topological properties of its vertices { x F k } and edges {E k,m } . With respectto the orientable hypersurface Γ, described by the level-set function given in eq. (3), the logical status S of any of avertex is either interior ( S = − S = 0) or exterior ( S = 1). Implementation detail (Robustness) . For the purpose of numerical robustness, the status assignment employs a ubular neighborhood of thickness (cid:15) zero around Γ , corresponding to the interval ( − (cid:15) zero , (cid:15) zero ) , i.e. S ( x F k ) := if | φ Γ ( x F k ) | < (cid:15) zero sign (cid:16) φ Γ ( x F k ) (cid:17) if | φ Γ ( x F k ) | ≥ (cid:15) zero . (23) In other words, any point whose absolute distance to Γ falls below (cid:15) zero is considered to be on Γ . The choice of anappropriate tolerance strongly depends on the absolute value of the characterisic length scale h of the polyhedron,especially for h (cid:28) . Throughout this work, we have h ≈ and, thus, let (cid:15) zero := 10 − . In fact, all zero-comparisonsare implemented in this way. The hierarchically superior entity is an edge E , whose status is a surjective function of the status of its associatedvertices { ν i , ν j } . An edge with at least one of its vertices located on Γ becomes degenerate, where the type of degen-Table 2: Edge status as function of the stat¯us of the associated vertices; cf. note 2.2. degenerate exterior interior intersected exterior interior intersected S ( ν i ) 1 − ± − S ( ν j ) 1 − ∓ − S ( E ) 1 − − n Γ exterior interior b b b bbb x F k, x F k, x F k, x F k, x F k, x F k, eration depends on the status of the respectively complementary vertex; cf. table 2. The length of an interior/exterioredge E k,m (including the respective degenerate cases) with vertices { ν i , ν j } is trivially computed as l k,m = S ( E k,m ) ∈ { , } , (cid:107) ν i − ν j (cid:107) if S ( E k,m ) ∈ {− , − } . (24)For an intersected edge ( S ( E k,m ) = 0), one obtains the monotonous piecewise linear function l k,m ( s ) = (cid:104) ν i − x Γ , , n Γ (cid:105)− s (cid:104) ν i − ν j , n Γ (cid:105) (cid:107) ν i − ν j (cid:107) if S ( ν i ) = − , (cid:104) ν j − x Γ , , n Γ (cid:105)− s (cid:104) ν j − ν i , n Γ (cid:105) (cid:107) ν i − ν j (cid:107) if S ( ν j ) = − , (25)whose derivative with respect to the signed distance s , i.e. l (cid:48) k,m = S ( E k,m ) ∈ {± , } , (cid:107) ν i − ν j (cid:107)|(cid:104) ν i − ν j , n Γ (cid:105)| if S ( E k,m ) ∈ { , − } , (26)depends only on static geometric quantities, where the potential discontinuities of A (cid:48)(cid:48) k ( s ) (cid:12)(cid:12) s = s i can be removed bydifferent choices of the status assignment. Figure 8 illustrates the concept, whose main implications can be cast asfollows:1. The assignment of the edge length derivative based on the edge status corresponds to choosing the left- (eq. (26);15 Γ n Γ n Γ (a) degenerate face intersection configurations; cf. table 2 . . . . s A k |F k | . . .
25 ˜ s A ′ k |F k | . . − . .
15 ˜ s r r r A ′′ k |F k | r A ′′ k ( s i ) |F k | (b) normalized interior face area (eq. (18)) and derivatives (eq. (22)) as function of normalized position ˜ s , where the dashed vertical linescorrespond to the configurations in (a) (left to right) Figure 8: Degenerate intersection configurations and corresponding discontinuities of second derivative of face area(eq. (22)) for polygonal face F k with vertices { (1 , , (7 , , (11 , , (5 , , (3 , } , base point x Γ , = [1 , T and PLICnormal n Γ = [4 , T √ . (cid:4) in figure 8(b)) or right-sided limit, i.e. A (cid:48)(cid:48) k ( s i ) := lim s (cid:37) s i A (cid:48)(cid:48) k ( s ) or A (cid:48)(cid:48) k ( s i ) := lim s (cid:38) s i A (cid:48)(cid:48) k ( s ) . (27)In theory, each interior discontinuity s − < s i < s + may be removed by an individual linear combination of therespective left- and right-sided limit. However, within this work, we choose the left-sided limits; cf. eq. (26).2. While, from a conceptual point of view, the continuation of A (cid:48)(cid:48) k is important, the numerical results shown insection 4 are barely affected by the choice of the status assignment, cf. eq. (26). This is to be expected, sincethe probability of the iteration producing some s n ∈ { s i } is zero.3. For a degenerately intersected edge E k,m , the assignment of l k,m (eq. (24)) and l (cid:48) k,m (eq. (26)) depends on theintersection status of its parenting face (index k ). Note 2.2 (Degenerately intersected faces ( S ( F k ) = 3)) . The reasoning behind the removal of the discontinuities of A (cid:48)(cid:48) k ( s ) analogously applies to V (cid:48)(cid:48) ( s ) , i.e. the assignment of the status of topologically inferior entities (faces) correspondsto the left- or right-sided limits. As for the faces, we choose the left-sided limits to resolve the interior discontinuities s − < s i < s + ; cf. figure 8. The boundary discontinuities { s − , s + } do not require to be removed for the problem under consideration here. .3 On the choice of the initial value s Numerical experiments have shown that, even for seemingly well-posed configurations of {P , n Γ , α ∗ } , computing theinitial value for eq. (10) by linear interpolation, i.e. s := s − + α ∗ s ∆ , may inhibit convergence to the desired root s ∗ . Hence, following Maric [16], the root of a global cubic approximation of the deviation D (fulfilling the boundaryconditions discussed in note 2.1) provides a suitable initial value. Plugging the derivatives of the volume fraction atthe lower and upper bound, say, α (cid:48)± := α (cid:48) ( s ± ), into the standard cubic spline given in eq. (36) yields D (˜ s ; α ∗ ) = (cid:16) α (cid:48)− s ∆ + α (cid:48) + s ∆ − (cid:17) ˜ s + (cid:16) − α (cid:48)− s ∆ − α (cid:48) + s ∆ (cid:17) ˜ s + α (cid:48)− s ∆ ˜ s − α ∗ with ˜ s = s − s − s ∆ , (28)where Maric [16] assumes α (cid:48)± ≡ . . s r s − α ∗ − α ∗ r D D (a) deviation D and cubic spline approxi-mation D for unit cube with normal n Γ = [1 , − , T √ and α ∗ = . α − = 0 α − = 0 α − = 0 α − = 0 α + = 0 α + = 0 α + = 0 α + = 0cuboid 99.95 0.00 0.00 0.05dodecahedron 100.00 0.00 0.00 0.00torus 100.00 0.00 0.00 0.00letter A 99.97 0.00 0.00 0.03tetrahedron 99.95 0.02 0.02 0.00 (b) Number of instances (in %, rounded) where the PLIC normal n Γ ∈ S n Γ (see subsection 1.2 for details) and face normal n F ,k are aligned for faces F k ⊂ ∂ H P ; cf. note 2.1. Figure 9: Illustration of global cubic approximation (left; cf. eq. (28)) and number of instances in S n Γ with non-trivialderivatives at the boundary (right).Despite the fact that the results shown in figure 9(b) strongly depend on the choice of S n Γ , this numercial ob-servation supports the assumption of Maric [16] ( α ± ≡ Newton -type method forthe computation of s for reasons of numerical error control, we prefer to pursue an analytical approach: since thediscriminant ∆( D ) = / α ∗ (1 − α ∗ ) is positive for 0 < α ∗ < D exhibits three distinct real roots, one of which islocated in [ s − , s + ]. From the formula of Vi`eta , one obtains s ( α ∗ ) := s − + s ∆ − cos (cid:32) arccos(2 α ∗ − − π (cid:33) . (29)While it is theoretically possible to impose non-trivial Neumann -type boundary conditions for D at s − and s + , thereare two major issues to consider:1. For general α ± (cid:54) = 0, the evaluation of the analytic formulae for the roots becomes numerically challenging, dueto the complexity of the involved expressions. Figure 10 shows the sign of the discriminant ∆( D ) for α ∗ = / ,which is decisive for the class of the roots and, hence, the respective analytical expression. In general, theevaluation of ∆( D ) requires to compute sixth powers of potentially very small numbers, which, in terms ofnumerical robustness, poses a highly complex task in itself. Their numerical computation, on the other hand,implies additional computational effort. Especially for α ∗ ≈ α ∗ ≈
1, a robust computation requires carefulimplementation.2. There are configurations { α (cid:48)− , α (cid:48) + } degrading the monotonicity of D , hence violating one of the properties of the(target) function, namely α (cid:48) = | Γ ∩P||P| >
0. Computing the derivative of eq. (28) yields a quadratic polynomial,17 1 2 3 4 5012345 α ′− s ∆ α ′ + s ∆ ∆ > ∆ > ∆ < D ) as a function of normalized boundary derivatives α (cid:48)± s ∆ .whose roots correspond to the critical points of D . After some simple manipulations, one obtains (cid:26) α (cid:48)± ∈ [0 , ∞ ) : 3 (cid:16) α (cid:48)− + α (cid:48) + − s − (cid:17) + ( α (cid:48)− − α (cid:48) + ) ≤ s − or ( α (cid:48)− + α (cid:48) + ) s ∆ ≤ (cid:27) (30)as the set of admissible configurations, being a composition of an ellipse with semiaxes s − ( √ , √
6) centered at s − (2 , T and rotated by π and a triangle; cf. figure 11. Despite that none of the numerical test cases presented0 1 2 3 4 5012345 α ′− s ∆ α ′ + s ∆ u r a d m i ss i b l e non-admissible n e u t r a l l y a d m i s s i b l e ×××× ×× . . . . s N : (1 , D + α ∗ D ′ . . . . s (cid:4) : (cid:0) , (cid:1) D + α ∗ D ′ Figure 11: Admissibility chart of normalized boundary derivatives α (cid:48)± s ∆ with illustrations of admissible ( (cid:78) , center)and non-admissible ( (cid:4) , right) approximations D .in section 4 (indicated by × ) suffers from being non-admissible, the admissibility requires to be probed in anyrobust implementation.Furthermore, numerical experiments have shown that the results presented in section 4 do not benefit from considering α ± (cid:54) = 0, once more supporting the assumption of Maric [16] discussed above. Thus, in what follows, we employ eq. (29)to obtain the initial iteration. ∆ s The analytic expression of the volume fraction obtained in eq. (20) implies that its derivatives with respect to thesigned distance s can be easily evaluated, suggesting an exploitation within a higher-order rootfinding algorithm. Thus,the present work introduces an algorithm based on implicit bracketing with locally quadratic approximation. Theupcoming subsections provide the key ideas underlying the numerical algorithm proposed in section 3.18 .4.1 Implicit bracketing Recall from subsection 2.1 that, within a bracket B i := [ s i , s i +1 ], the volume fraction α ( s ) is an increasing cubicpolynomial, denoted by S i ( s ). For any given s n ∈ B i , the polynomial reads S i ( z ; s n ) = α (cid:48)(cid:48)(cid:48) ( s n )6 (cid:0) z − s n (cid:1) + α (cid:48)(cid:48) ( s n )2 (cid:0) z − s n (cid:1) + α (cid:48) ( s n ) (cid:0) z − s n (cid:1) + α ( s n ) . (31)Hence, the truncation of the polyhedron P at any s n (implicitly) provides the full information of α ( s ) within thecontaining bracket B i . Exploiting that α i = S i ( s i ; s n ) and α i +1 = S i ( s i +1 ; s n ) suggests the following strategy:1. if the current iteration s n is not located in the target bracket B ∗ ( α ∗ < α i or α i +1 < α ∗ ), the next iteration isobtained from locally quadratic approximation (figure 12, left).2. if the current iteration s n is located in the target bracket B ∗ ( α i ≤ α ∗ ≤ α i +1 ), the sought s ∗ corresponds to theroot of S i − α ∗ , requiring no further truncation (figure 12, center). sα ∗ q u a d . a p p r o x . b r s i s i +1 s s α D + α ∗ sα ∗ α i α i +1 s i s i +1 e x a c t s p l i n e u ur s α sα ∗ α i α i +1 s i s i +1 e x a c t s p l i n e u ur s α D + α ∗ Figure 12: Implicit bracketing and locally quadratic approximation of volume fraction α ( s ); cf. subsection 2.4.2.Figure 12 illustrates the components of the strategy outlined above: with s (cid:54)∈ B i (left), a quadratic approximationyields s , which lies within the target bracket B ∗ = B i (center). The rightmost configuration already starts with s ∈ B ∗ , such that the spline interpolation directly yields the sought position s ∗ , implying that only a single truncationis required. As we shall see in section 4, configurations of this type are often obtained for α ∗ ≈ α ∗ ≈ Note 2.3 (Comparison to Maric [16]) . Figure 3 shows that the CCS method of Maric [16] also allows to obtain anexact interpolation of the volume fraction α . However, there are two major differences to the present algorithm:1. In general, the CCS method requires two subsequent iterations within the same bracket, i.e. { s n − , s n } ⊂ [ s i , s i +1 ] (assuming that s n > s n − without loss of generality).2. If two iterations { s n − , s n } ⊂ [ s i , s i +1 ] are found, the CCS method only exploits the exact interpolation in [ s n − , s n ] . Assume that s ∗ ∈ [ s i , s i +1 ] \ [ s n − , s n ] , i.e. the spline interpolation governs the sought s ∗ without asign change between the iterations, i.e. D ( s n − ) D ( s n ) > . Then, a linear approximation is applied to obtain s n +1 , implying another truncation, despite the fact that the target bracket B ∗ = [ s i , s i +1 ] was already found.As will be shown in section 4, these properties, in combination with a sufficiently small tolerance (cid:15) zero , yield an averagemargin of one polyhedron truncation for the proposed method in comparison to CCS. ote 2.4 (Comparison to Lopez et al. [14]) . While the present work resorts to the same basic building blocks asCIBRAVE, the novel method of volume computation introduced in subsection 2.1 allows for considerable improvementand a local formulation. The main differences of the algorithm of figure 15 and CIBRAVE are:
CIBRAVE present workinitial iteration s linear interpolation cubic interpolation (eq. (28)) analytical representation additional computation by-product of truncation of α ( s ) ( only at { s i } ) (at any s n ) update s n +1 linear interpolation quadratic approximationbased on { α ( s i ) } around s n (see 2.4.2)Let us note in passing that CIBRAVE requires topological connectivity on a cell level, whereas the present work onlyrequires a list of faces. Locally, at some given s n , the deviation D can be approximated by T D ( z ; s n ) = α (cid:48)(cid:48) ( s n )2 ( z − s n ) + α (cid:48) ( s n )( z − s n ) + α ( s n ) − α ∗ , (32)from whose roots (assuming, for the moment, that α (cid:48)(cid:48) ( s n ) (cid:54) = 0) we obtain[∆ s Taylor ]( s n ) = (cid:34) √ D − α (cid:48) α (cid:48)(cid:48) (cid:35) ( s n ) with the discriminant D = (cid:104) α (cid:48) α (cid:48) − α − α ∗ ) α (cid:48)(cid:48) (cid:105) ( s n ) . (33)Exploiting α ( s ∗ ) = α ∗ and enforcing s ∗ to be a fixed point of the iteration, i.e. ∆ s ( s ∗ ) ! = 0, eliminates the quadraticambiguity, i.e. −√ D − α ∗ α (cid:48)(cid:48) is ruled out in favor of the expression given in eq. (33). For vanishing second derivative, the Taylor step degenerates to the one obtained by linear approximation, i.e.lim α (cid:48)(cid:48) → (cid:2) ∆ s Taylor (cid:3) ( s n ) = (cid:20) − α − α ∗ α (cid:48) (cid:21) ( s n ) = [∆ s Newton ] ( s n ) . (34)Figure 13 provides an illustration for a cuboid, where the flowchart in figure 14 gathers the details of the numericalimplementation. 0 0 . . s − α − α ∗ α ′ α ′′ . . s − α ∗ − α ∗ α − α ∗ ∆ s Newton ∆ s Taylor
Figure 13: Volume fraction deviation D and derivatives for a cuboid with n Γ = [1 , − , T √ and α ∗ = (left) andassociated Newton (eq. (34)) and
Taylor (eq. (33)) steps (right).20 nput s n { α, α ′ , α ′′ } ( s n ) discriminant D ( s n ) ≥ degenerate α ′′ = 0? Newton ∆ s := − α − α ∗ α ′ Taylor ∆ s := √ D − α ′ α ′′ admissiblity s n + ∆ s ∈ [ s − , s + ] ? output ∆ s yes nono yesyesno(fallback)Figure 14: Flowchart of Taylor step computation.
Note 2.5 (Locally third-order approximation) . While it is theoretically possible to apply a third-order approximation,we restrain from doing so for two reasons:1. Numerical experiments have shown that, in general, | α ( s ) − S i ( s ; s n ) | becomes prohibitively large for s (cid:54)∈ B i ,being attributable to the insufficient regularity of α ( s ) (right); cf. note 2.1. This disadvantageous effect prevailsprominently in the beginning of the iteration.2. In cases where the third order approximation proves sufficiently accurate, the second order approximation providessimiliar results, however without the requirement for numerical computation of the root (left); cf. subsection 2.4.2below. r α second orderthird order r α second orderthird order Figure 15 provides a flowchart of the algorithm. After reading the input parameters {P , α ∗ , n Γ } , representing thepolyhedron, volume fraction and normal, respectively, the static measures and coefficients (cf. eqs. 14, 19, 24-26) arecomputed and the signed distances of the vertices s i are sorted in ascending order (see section 2 for details). Startingwith the initial value s given in subsection 2.3, within each iteration (superscript n ), the truncation of the polyhedronyields the volume fraction and its derivatives. From the index i of the bracket containing the current iteration s n , thegoverning cubic polynomial S i is assembled, allowing to evaluate the interval [ α i , α i +1 ]. For α ∗ ∈ [ α i , α i +1 ], implying21hat s n ∈ B ∗ , the sought position s ∗ corresponds to the root of the cubic polynomial S i , which is obtained numericallyresorting to the algorithm described in appendix A.1. In the case where B i does not contain the sought position s ∗ ,the locally quadratic approximation of subsubsection 2.4.2 provides the next iteration s n +1 . By design, the quadraticapproximation may produce an update s n +1 located in a different bracket than the current iteration s n . Because theadmissibility ( s n +1 ! ∈ [ s − , s + ]) cannot be enforced a priori, an appropriate check is performed, resorting to a Newton step if s n +1 (cid:54)∈ [ s − , s + ]; cf. figure 14. It is possible that, after an update with ∆ s Taylor , the deviation D ( s n +1 ) becomes output s ∗ fallback(bisection)input P , α ∗ , n Γ static coefficients { γ V ,k , γ V ,k , γ A ,km , γ A ,km } sort s i measures { l ′ k,m , |F k |} initialize n := 0 s := s ( α ∗ ) intersect P − ( s n ) evaluate { α, α ′ , α ′′ } ( s n ) converged | α − α ∗ | < ǫ zero ?assemble S i α i := S i ( s i ; s n ) α i +1 := S i ( s i +1 ; s n ) find bracket index j s.t. s n ∈ [ s i , s i +1 ] final bracket α i ≤ α ∗ ≤ α i +1 ?compute step ∆ s Taylor interpolation s ∗ := root( S i ; s n ) n < n max ? update s n +1 := s n + ∆ s Taylor n := n + 1yesno no yes noyes s ∗ := s n Figure 15: Flowchart of the main algorithm. For all numerical testcases reported here, a tolerance of (cid:15) zero = 10 − was used.smaller than (cid:15) zero . These cases emerge for a plane Γ being almost parallel to one of the faces F k for α ∗ ≈ α ∗ ≈ { s i } of signed distances of the respective face vertices, implying a very smallsize of the target bracket B ∗ . In these cases, a quadratic approximation around some s n (cid:54)∈ B ∗ , but in the vicinity of B ∗ , is likely to be sufficiently accurate. However, their incidence barely reaches statistical significance. Should theiteration not yield the sought s ∗ with the prescribed number of iterations n max = 100, the algorithm resorts to abisection method. Remark 3.1 (Implementation complexity) . As can be seen from the flowchart in figure 15, the algorithm presented inthis paper heavily relies on the volume computation descibed in 2.1. Despite its face-based formulation, an integrationinto an existing legacy code, such as, e.g., openFOAM, requires some implementation effort. While the presentedmethod gets along with fairly basic geometric operations, the CIBRAVE algorithm of Lopez et al. [14] moreover re-quires establishing the topological connectivity of truncated polyhedral cells, which induces significant implementationcomplexity. In contrast, the CCS method of Maric [16] only resorts to the derivative of α (cid:48) , which is proportinal to | Γ ∩ P| . For all methods, the implementation complexity of the rootfinding routines is negligible in comparion thegeometric operations. In theory, the
Newton step may also fail to produce an admissible update. However, we have not encountered such a configuration. mplementation detail (Object orientation) . The algorithm presented in this paper was implemented in Fortran2008 (8-byte real/4-byte integer), employing object-orientation wherever possible. In figure 15, light blue boxes belongto the root-finding part and the dark blue boxes correspond to methods of the abstract class polyhedron ; cf. subsection 1.3.Since the computation of the measures and coefficients (eqs. 14, 19, 24-26) can be implemented in a generic way, theconcretization of the derived classes ( cuboid etc., cf. table 1) only requires to implement the polyhedron definition, i.e.the assignment of vertices and topology.
In what follows, we assess the algorithm presented in section 3 for a comprehensive set of instances {P , α ∗ , n Γ } , wheretable 1 gathers the classes of polyhedra P and subsection 1.2 defines the set of volume fractions α ∗ ∈ S α ∗ and normals n Γ ∈ S n Γ . Separate processing and aggregation of the results highlights the influence of α ∗ and n Γ in subsections 4.1and 4.2, respectively. Before discussing the numerical experiments in detail, however, it is instructive to illustrate thevolume fraction and its derivatives as a function of the plane position s for some selected instances in figure 16. Thestrong variance of the third derivatives in figure 16(c) vividly illustrates the argument of note 2.5, stating that, ingeneral, a third order approximation S i ( s ; s n ) around some s n ∈ B i produces large deviations | α ( s ) − S i ( s ; s n ) | outsidethe bracket containing s n , i.e. for s (cid:54)∈ B i ; cf. eq. (31). A suitable and robust measure for performance assessment.
While total absolute execution times constitutean intuitive measure for the computational performance of an algorithm, the problem at hand intrinsically featuresa more suitable indicator: most of the literature contributions discussed in subsection 1.1 state that the polyhedrontruncation accounts for the largest share in computation time; see, e.g., Lopez et al. [14, fig. 16] and Maric [16,fig. 8]. This property becomes especially prominent for non-convex polyhedra, where the truncation may producedisconnected sub-polyhedra, whose convex decomposition poses a challenging and computationally expensive taskon its own. Despite the fact that the method proposed in subsection 2.1 significantly reduces the computationaleffort by avoiding the usage of topological connectivity, the inherent complexity of the geometric operations remainspredominant. Thus, we assess the performance in terms of the numbers of truncation operations per instance,i.e. N trunc = N trunc ( P , α ∗ , n Γ ), corresponding to the final value of n in the flowchart of figure 15. Furthermore, weemploy the execution times t ( P , α ∗ , n Γ ) to establish a qualitative idea of the relative differences between the classesof polyhedra. The execution time is obtained by subtracting the results of MPI WTIME() , called respectively at thebegin and end of the subroutine implementing figure 15, i.e. it includes all secondary operations, e.g., sorting { s i } . Comparison to literature.
Subsection 1.2 and table 1 outline the structure of numerical experiments conductedwithin this section and define the sample set for the volume fraction α ∗ and normal orientation n Γ . Out of theliterature contributions discussed in subsection 1.1, only the manuscripts of Lopez et al. [14] and Maric [16] containdisaggregated data suitable for a meaningful comparison, where only Maric [16] contains data for non-convex polyhedra.For example, Lopez et al. [15] only report relative CPU times, normalized by the average runtime for a quadrangularcell. Also, they do not comment on the coverage of the measurements in terms of algorithmic elements. Thus, basedon the data published in their manuscript, no sensible comparison can be made. Furthermore, Lopez et al. [14] employrandomly distributed normals n Γ , while Maric [16] resorts to the sample set given in subsection 1.2. For the algorithm given in figure 15, the number of truncations corresponds to the number of iterations, i.e. N trunc = N iter . processor: Intel Core i7-7700HQ, compiler: gfortran 8.1.0 (with optimization -O3 ).
23 0 . . . . . −
101 ˜ s k α ′ k ∞ = 1 . k α ′′ k ∞ = 6 . k α ′′′ k ∞ = 183 . α α ′ α ′′ α ′′′ (a) dodecahedron . . . . . −
101 ˜ s k α ′ k ∞ = 1 . k α ′′ k ∞ = 12 . k α ′′′ k ∞ = 6486 . α α ′ α ′′ α ′′′ (b) torus . . . . . −
101 ˜ s k α ′ k ∞ = 1 . k α ′′ k ∞ = 45 . k α ′′′ k ∞ = 1733 . α α ′ α ′′ α ′′′ (c) letterA Figure 16: Illustration of intersected polyhedra with volume fraction as function of normalized parameter ˜ s . Thedashed vertical lines indicate the bracket boundaries s i , where the circles indicate the potential discontinuity at s i .For the reader’s convenience, the derivatives are scaled by the respective (cid:107)·(cid:107) ∞ .24 .1 Influence of the volume fraction α ∗ Figure 19 gathers the aggregated results, where the figures in the respective rows share the following properties:i. The volume fractions α ∗ on the abscissa are subject to a logarithmic stretch in the vicinity of zero and one,respectively. Between 10 − and 1 − − , the scaling remains linear.ii. The averages are defined as N av ,α ( P , α ∗ ) := 1 N n Γ N n Γ (cid:88) k =1 N trunc ( P , α ∗ , n Γ ,k ) and N av ( P ) := 1 N α ∗ N α ∗ (cid:88) m =1 N av ,α ( P , α ∗ m ) , (35)where N n Γ = size( S n Γ ) = 2 M n Γ ( M n Γ + 1) and N α ∗ = size( S α ∗ ) = M α ∗ + 10; cf. subsection 1.2 for details. Theexecution times t av ,α and t av are defined analogously. For the reader’s convenience, the left column containstheir normalized counterpart, with the absolute scale t av provided at the top of the respective figure.iii. The solid dots ( • ) represent the performance measures, given in eq. (35), obtained by the proposed algorithm,where the shaded regions indicate a tube with a width of twice the mean standard deviation. If available, thesquares ( (cid:4) ) and triangles ( (cid:78) ) show the results of Lopez et al. [14, fig. 17(d+e)] and Maric [16], respectively,for comparison. Lopez et al. [14] obtain the relative execution times (last row in table 3) from 10 randomcombinations of α ∗ and n Γ .The main implications of figure 19 can be cast as follows:1. For all polyhedra under consideration, the average number of truncations is roughly symmetric to α ∗ = / ,reflecting the symmetry of the positioning problem, since {P , α ∗ , n Γ } corresponds to {P , − α ∗ , − n Γ } . Ingeneral, however, the inverted order of { s i } spoils a corresponding symmetry of the initial iteration s , obtainedfrom the cubic interpolation D given in eq. (28).2. The average number of truncations differs only slighty for the different polyhedra, which indicates a robustperformance of the algorithm. The same holds for the relative execution times (left column in figure 19).3. For the cuboid, tetrahedron, dodecahedron and torus, there are no instances for which the fallback strategy hasto be applied. For the letter A, out of all instances in subsection 1.2, there are 8 (approximately 0 . ‰ ) fall-backcases. Figure 17 highlights the reasons: in essence, the applied second-order approximation bounces back andforth in a loop between two positions. Configurations of this type cannot be excluded a priori. However, thisfallback-quota can safely be considered acceptable. s s ∗ r brbr b D r quadratic(a) no convergence s s ∗ r brb D r linear(b) convergence after 3 iterations Figure 17: Comparison of convergence using (a)
Taylor and (b)
Newton step for polyhedron letter A: while the
Newton method converges to the sought root within 3 iterations, the quadratic approximation bounces back andforth without convergence (indicated by (cid:4) ). In fact, all fallback cases reported for letter A are of this type.4. For all polyhedra under consideration here, table 3 reports between 1 .
13 and 1 .
52 truncations to find s ∗ , wherethe relative execution times clearly reflect the geometric complexity of the polyhedra. The anticipation concern-ing the differences to CCS, formulated by note 2.3, find their substantiation in figure 19: on average and for all25able 3: Aggregated execution statistics for proposed method; cf. figure 19 and table 3 in Lopez et al. [14, p. 353]. cuboid dodecahedron torus letterA tetrahedron N av .
13 1 .
18 1 .
52 1 .
46 1 . t av ( µ s) 2 . . . . . .
00 2 .
27 11 .
58 2 .
84 0 . α ∗ under consideration, CCS requires one more truncation than the proposed algorithm. This can be attributedto the consecutive character of the spline interpolation, requiring, in general, two interior truncations, i.e. at s − < s n , s n +1 < s + . While figure 19 only compares CCS for convex polyhedra, the manuscript of Maric [16] alsocontains a parametrized set of non-convex polyhedra, for which CCS exhibits robust performance and qualita-tively equivalent results. In comparison, Lopez et al. [14] only contains disaggregated data for convex polyhedraand Lopez et al. [15] contains no disaggregated data at all. For the tetrahedron, constituting the simplest possiblepolyhedron, the average number of truncations produced by CIBRAVE essentially coincides with those of theproposed algorithm. However, for the cuboid, the proposed method outperforms CIBRAVE by approximately0 . n Γ in Lopez et al. [14] is randomly computed, implying that degenerate cases(the plane Γ is almost parallel to one of the faces F k , i.e. |(cid:104) n Γ , n F ,k (cid:105)| ≈ Newton -based method, similar to the algorithm ofChen and Zhang [4]. As can be seen in the prototypical illustration of figure 18(a), besides the significant26101418 α ∗ N av ,α cuboid − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − N av b b b b b bbbb b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b bbbb b b b b b α ∗ N av ,α dodecahedron − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − N av b b b b b bbbb b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b bbbb b b b b b (a) average number of truncations as a function of α ∗ cuboid dodecahedron torus letterA tetrahedron N av .
37 5 .
39 5 .
23 5 .
78 7 . t av ( µ s) 4 . . . . . .
00 2 .
18 10 .
12 2 .
76 0 . (b) aggregated execution statistics Figure 18: Average number of truncations for
Newton -based method of Chen and Zhang [4].overall increase in the number of truncations by a factor of 4 to 5, the
Newton -based method also requiresan increasing number of truncations for α ∗ ≈ α ∗ ≈
1. This numerical artifact can be attributed to thediminishing derivatives α (cid:48) , resembling a well-known drawback of Newton -based methods. In a prototypicalsense, this substantiates the superiority of spline-based methods over methods based on local first- or second-order approximation of α ( s ). Among the different concepts of algorithmic incorporation, e.g., explicit bracketing26as employed by Lopez et al. [14] in CIBRAVE) or consecutive interpolation (as employed by Maric [16] in CCS),the impact on the average number of truncations is small in comparison to approximation-based methods.012 α ∗ t av ,α / t av t av = 2 . × − s cuboid − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b α ∗ N av ,α − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − N av u u u u u uuuu u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u uuuu u u u u ur r r r r r r r r r rb b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b cuboid α ∗ t av ,α / t av t av = 5 . × − s dodecahedron − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b α ∗ N av ,α − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − N av u u u u u uuuu u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u uuuu u u u u ub b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b dodecahedron α ∗ t av ,α / t av t av = 29 . × − s torus − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b α ∗ N av ,α − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − N av b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b torus α ∗ t av ,α / t av t av = 7 . × − s letterA − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b α ∗ N av ,α − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − N av b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b letterA α ∗ t av ,α / t av t av = 1 . × − s tetrahedron − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b α ∗ N av ,α − − −
91 0 − − −
81 0 − − −
71 0 − − −
61 0 − − −
51 0 − − − N av u u u u u uuuu u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u uuuu u u u u ur r r r r r r r r r rb b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b tetrahedron Figure 19: Average positioning time in µ s and number of truncations for polyhedra from table 1 with reconstructiontolerance (cid:15) = 10 − obtained with proposed algorithm ( • ), CCS ( (cid:78) ) and CIBRAVE ( (cid:4) ). Figure 20 associates the average number of truncations of the proposed method to the respective orientation of thenormal n Γ , suggesting a parametrization in spherical coordinates. In general, the spatial distribution of the peaksin the average truncations resembles the geometry of the underlying polyhedron P , eventually inheriting some of itssymmetries. Figure 20(a) vividly illustrates this for the cuboid: an orientation with θ ∈ { , π / , π } corresponds to adegeneration to two dimensions for arbitrary ϕ , implying that the associated volume fraction α becomes quadratic.27ince the proposed method locally applies a quadratic approximation, the computation of the sought position s ∗ requires less truncations. On the contrary, for θ (cid:54)∈ { , π / , π } , the volume fraction is composed of non-degenerate cubicpolynomials in general, imposing an increased computational effort to find s ∗ .0 60 120 180 240 300 360060120180 ϕθ .
00 1 . (a) cuboid ϕθ .
00 1 . (b) dodecahedron ϕθ .
00 1 . (c) letterA ϕθ .
00 1 . (d) tetrahedron Figure 20: Average number of truncations over normal orientation in spherical coordinates.
We have introduced an algorithm for the positioning of a PLIC plane in an arbitrary polyhedron. The recursiveapplication of the
Gaussian divergence theorem in its respectively appropriate form allows to obtain a highly efficientface-based computation of the volume of the truncated polyhedron, implying that no connectivity information has to beestablished at runtime. In combination with the precomputability of the involved coefficients, the face-based characterrenders the presented approach most suitable for parallel computations on unstructured meshes. Furthermore, werigorously derive the mathematical foundations and examine the relevant limiting cases. The root-finding componentrelies on a combination of bracketing and local second-order approximation, allowing to compute the sought positionto a user-defined precision ( (cid:15) zero = 10 − ). For both the main algorithm and its secondary components, we haveprovided flowcharts along with a comprehensive description. We have assessed our algorithm on a set of convex andnon-convex polyhedra of genus zero and one (i.e., with a hole), which, to this extent, has not been done before. Forall numerical instances, comprising a wide range of volume fractions 10 − ≤ α ∗ ≤ − − and normal orientations n Γ ∈ S , one to two truncations are required on average. The variations with respect to α ∗ and n Γ are small,indicating the robustness of the algorithm. In terms of the average number of truncations, being accountable for themajor part of the computational effort, the proposed algorithm outperforms existing methods. Thus, we draw thefollowing conclusions:1. The recursive application of the Gaussian divergence theorem in appropriate form allows for an efficient compu-tation of the volume of a truncated arbitrary polyhedron , along with its derivatives. This face-based problemdecomposition allows to avoid extracting topological connectivity , which is advantageous in terms of bothimplementation complexity and computational effort.28. The local knowledge of the first to third derivatives substantially accelerates the iterative root-finding .For a comprehensive set of combinations of polyhedra, volume fractions and normal orientations, one to twotruncations are required on average to position the PLIC plane.
References [1] Hyung Taek Ahn and Mikhail Shashkov. Multi-material interface reconstruction on generalized polyhedral meshes.
Journal of Computational Physics , 226:2096–2131, 2007.[2] Kendall E. Atkinson.
An introduction to numerical analysis . John Wiley & Sons, 1988.[3] Richard P. Brent.
Algorithms for minimization without derivatves . Prentice-Hall, 2002.[4] Xiang Chen and Xiong Zhang. A predicted-newton’s method for solving the interface positioning equation inthe mof method on general polyhedrons.
Journal of Computational Physics , 384:60 – 76, 2019. ISSN 0021-9991.doi: https://doi.org/10.1016/j.jcp.2018.12.038. URL .[5] Frank H. Clarke.
Optimization and Nonsmooth Analysis . Society for Industrial and Applied Mathematics, 1990.doi: 10.1137/1.9781611971309. URL https://epubs.siam.org/doi/abs/10.1137/1.9781611971309 .[6] Dezhi Dai and Albert Y. Tong. Analytical interface reconstruction algorithms in the plic-vof method for 3dpolyhedral unstructured meshes.
International Journal for Numerical Methods in Fluids , 91(5):213–227, 2019.doi: 10.1002/fld.4750.[7] S. Diot, M.M. Francois, and E.D. Dendy. An interface reconstruction method based on analytical formulae for 2dplanar and axisymmetric arbitrary convex cells.
Journal of Computational Physics , 275:53 – 64, 2014. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2014.06.060. URL .[8] Steven Diot and Marianne M. Fran¸cois. An interface reconstruction method based on an analytical formula for 3darbitrary convex cells.
Journal of Computational Physics , 305:63 – 74, 2016. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2015.10.011. URL .[9] Ronald Goldman. Area of planar polygons and volume of polyhedra.
Graphics Gems II , pages 170–171, 2004.doi: 10.1016/B978-0-08-050754-5.50043-8.[10] Denis Gueyffier, Jie Li, Ali Nadim, Ruben Scardovelli, and St´ephane Zaleski. Volume-of-fluid interface trackingith smoothed surface stress methods for three-dimensional flows.
Journal of Computational Physics , 152:423–456,1999.[11] Douglas Kothe, W. Rider, Stewart Mosso, J. Brock, and John Hochstein. Volume tracking of interfaces havingsurface tension in two and three dimensions.
AIAA Paper , 96, 03 1996. doi: 10.2514/6.1996-859.[12] Joaquin Lopez and Julio Hernandez. Analytical and geometrical tools for 3D volume of fluid methods in generalgrids.
Journal of Computational Physics , 227(12):5939–5948, 2008. doi: https://doi.org/10.1016/j.jcp.2008.03.010. URL .[13] Joaquin Lopez, Pablo Gomez, Claudio Zanzi, Felix Faura, and Julio Hernandez. Application of non-convexanalytic and geometric tools to a PLIC-VoF method. In
International Mechanical Engineering Congress andExposition . ASME, November 2016. 2914] Joaquin Lopez, Julio Hernandez, Pablo Gomez, and Felix Faura. A new volume conservation enforcement methodfor PLIC reconstruction in general convex grids.
Journal of Computational Physics , 316:338–359, 2016.[15] Joaquin Lopez, Julio Hernandez, Pablo Gomez, and Felix Faura. Non-convex analytical and geometrical toolsfor volume truncation, initialization and conservation enforcement in vof methods.
Journal of ComputationalPhysics , 392:666–693, 2019.[16] Tomislav Maric. Iterative Volume-of-Fluid interface positioning in general polyhedrons with Consecutive CubicSpline interpolation. arXiv preprint , 2012.02227, 2020. URL https://arxiv.org/abs/2012.02227 .[17] David E. Muller. A method for solving algebraic equations using an automatic computer.
Mathematical Tablesand other aids to Computation , 10(56):208–215, 1956.[18] William J. Rider and Douglas B. Kothe. Reconstructing volume tracking.
Journal of Computational Physics ,141:112–152, 1998.[19] Ruben Scardovelli and Stephane Zaleski. Analytical relations connecting linear interfaces and volume fractionsin rectangular grids.
Journal of Computational Physics , 164:228–237, 2000.[20] Jonathan Richard Shewchuk. Lecture notes on geometric robustness. Technical report, University of Californiaat Berkeley, 2013.[21] Maciej Skarysz, Andrew Garmory, and Mehriar Dianat. An iterative interface reconstruction method for PLICin general convex grids as part of a coupled level set volume of fluid solver.
Journal of Computational Physics ,368:254–276, 2018.[22] Josef Stoer and R. Bulirsch.
Introduction to Numerical Analysis . Springer-Verlag New York, 2002.[23] G. Strang. On the construction and comparison of difference schemes.
SIAM Journal of Numerical Analysis , 5:506–517, 1968.[24] Herbert Voß.
PSTricks: Grafik mit PostScript f¨ur TEX und L A TEX , volume 7. Lehmanns Media, 2016. ISBN9783865412805.[25] Xiaofeng Yang and Ashley J. James. Analytic relations for reconstructing piecewise linear interfaces in triangularand tetrahedral grids.
Journal of Computational Physics , 214(1):41–54, 2006. doi: 10.1016/j.jcp.2005.09.002.[26] Chuck Zemach. Notes on the volume of a ruled hexahedron behind a truncating plane (unpublished). Los AlamosNational Laboratory, 1993.
Acknowledgment
The authors gratefully acknowledge financial support provided by the German Research Foundation (DFG) withinthe scope of SFB-TRR 75 (project number 84292822). Furthermore, the authors would like to thankDr.-Ing. Tomislav Mari´c for his advice concerning the design of numerical experiments.30
Supplementary information
A.1 Rootfinding for a cubic polynomial
The cubic spline tangentially interpolating the nodes { p , p } with p i = [ x i , f i , f (cid:48) i ] and ∆ x = x j − x i reads S ( x ; p i , p j ) = ∆ x ( f (cid:48) i + f (cid:48) j ) − f j − f i )∆ x ( x − x i ) + 3( f j − f i ) − ∆ x (2 f (cid:48) i + f (cid:48) j )∆ x ( x − x i ) + f (cid:48) i ( x − x i ) + f i . (36)In the non-gegenerate case (i.e., the cubic coefficient is non-zero), the root of eq. (36) is numerically computed usingthe algorithm given in figure 21, while in the (quadratic or linear) degenerate case the respective root is obtainedanalytically. Unless stated otherwise, the initial value is obtained by linear interpolation, i.e. x := x i − f i ∆ xf j − f i . Lopez output x n fallback(bisection)input x n := 0 evaluate {S , S ′ , S ′′ } ( x n ) |S | ≤ ǫ zero ? discriminant D ( x n ) ≥ n < n max ? n := n + 1 S ′′ = 0? Taylor x n := x n + √ D −S ′ S ′′ S ′ = 0? Newton x n := x n − S S ′ yesno yes yesyes no yesnono noFigure 21: Flowchart of the cubic rootfinding algorithm.and Hernandez [12] state that ” [...] explicit formulae for cubic roots generally consume less computational time thanNewton–Raphson iterations “. However, they (i) do not specify the gain in terms of computational time and (ii) reportseizing robustness for α ∗ < − and α ∗ > − − . Thus, we prefer to resort to the algorithm in figure 21. A.2 Non-standard polyhedra
Torus.
The surface of a torus is parametrized by F ( ϕ, θ ; R, γ ) = R (cid:2) (1 + γ cos θ ) cos ϕ, (1 + γ cos θ ) sin ϕ, γ sin θ (cid:3) T , (37)where R is the major radius and γ ∈ [0 ,
1] is the minor radius ratio. The number of major (minor) nodes is denoted N ( N ), such that the vertices of the torus are (cid:40) F (cid:18) iπN , jπN ; R, γ (cid:19)(cid:41) N j =1 N i =1 , (38)where, throughout his work, we employ γ = / , R = 1, N = 9 and N = 7; see figures 16(b) and 22(a) for anillustration. 31 etter A. The vertices of the letter A are / { (0 , , d ) , (4 , , d ) , (6 , , d ) , (8 , , d ) , (10 , , d ) , (14 , , d ) , (10 , , d ) , (8 , , d ) } and / { (6 , , d ) , (8 , , d ) , (7 , , d ) } , with d ∈ { , } corresponds to the front and back face, respectively. Figure 22(b) provides an illustration. The hole iscreated by an inverted arrangement of the blue and red set of vertices, resulting in two coplanar faces with opposingnormals; cf. figure 5. bbbbb b b brrrrrr r r r r N N (a) Torus with relevant quantities. b b b b b bbb bb b (b) Non-convex base faces ofpolyhedron letterA ..