Fast advancing layer method for viscous mesh generation

Hongfei YE, Jianjun CHEN,*, Yufei PANG, Yang LIU, Yao ZHENG

a School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China

b Center for Engineering and Scientific Computation, Zhejiang University, Hangzhou 310027, China

c China Aerodynamics Research and Development Center, Mianyang 621000, China

KEYWORDS

Abstract An efficient Advancing Layer Method (ALM) is presented to create semi-structured prisms on viscous walls, in which a procedure that checks possible front intersections is essential to its efficiency.This paper develops various novel schemes to improve the algorithm’s efficiency precisely while not sacrificing its robustness and the resulting mesh quality.First, it employs a set of new techniques, and data structures are developed to improve the efficiency of the frontcheck procedure.Then,within each octant,a new filter is developed to reduce the intersection computations in the searching process.In addition, data structures are well designed to store the contiguously accessed data in each computing-intensive loop in a contiguous space for a potentially better cache hit ratio.We built a geometry model library formed by examples of industrial complexity to demonstrate the practicability of the algorithm.All the efforts mentioned above enable us to reduce the percentage of computing time taken by intersection check to an acceptable level(approximately 26%), which make it no longer be the most time-consuming part.

1.Introduction

Automatic generation of high-quality meshes remains one of the major issues when conducting high-Reynolds viscous flow analysis.Various techniques have been developed to address this issue, and the performance of these techniques in terms of viscous accuracy and ease of use has been analyzed and compared at length.1–3It was suggested that, among different techniques, the generation of a prismatic hybrid mesh represents a good compromise between viscous accuracy and ease of use.4In this mesh, the near field of viscous walls (referred to as boundary layers hereafter) is configured with layered prismatic elements to resolve high flow gradients normal to the walls.In contrast, the remaining domain and the surface geometry are filled with unstructured meshes so that the mesh generation pipeline could be automated.

The generation of a prismatic hybrid mesh contains two individual steps, i.e., creating layered prisms in the near-field region and creating unstructured elements in the far-field region.The generation of unstructured elements, by either the Advancing Front Technique (AFT)5or the Delaunay Triangulation (DT) approach,6–8has been well studied, whereas the robust and efficient generation of high-quality boundary layer elements remains a challenging task.Presently, the most prevailing approach is the ALM.9–14In its simplest form, the ALM defines an initial front on the viscous walls and computes a marching direction at each front point.A layer of prismatic elements is then created by setting a marching distance along each direction and extending all front triangles to prisms.The entire boundary layer mesh is formed by repeating the above procedure a few times.15,16Nevertheless, this simplest form of the ALM cannot deal with configurations in which the surface geometry contains small gaps and/or complex corners (such as concave/convex corners or their combinations).Various approaches17–20are thus developed to enhance the performance of the ALM.

Among various issues of boundary layer mesh generation,our focus in this study is on the front intersection issue.Fig.1 presents two typical cases of intersected fronts,i.e.,local intersections that occur in regions close to concave corners and global intersections that occur in small gaps between two close viscous walls.A precautionary measure that can reduce the possibility of front intersections is to reduce marching distances at regions where front intersections possibly occur.11However, this technique cannot completely remove the possibility of front intersections.Special treatments on front intersections are still necessary to ensure an industry-level robustness of the algorithm.Note that the front intersections usually occur at local regions.A pre-processing procedure is suggested to classify the fronts into two groups.12,21It is assumed that the prisms attached to one group of the fronts(e.g., the fronts in the flat region) are free of intersections;therefore,these elements can be generated without intersection checks.Nevertheless, the elements attached to the rest of the fronts are(e.g.,the fronts at corners and/or small gaps)generated by using the standard ALM.This approach improves the efficiency by reducing the number of elements involved in intersection checks, whereas the cost is that element quality at complex corners is compromised, as observed in the mesh examples reported in Ref.12

It is worth mentioning that the approaches other than the ALM usually have to tackle the front intersection issue as well.For the approach that defines each layer of fronts on the isosurfaces of a solution field or its variants,22–25the fronts should not intersect each other if the iso-surfaces are accurately defined.However, numerical errors are likely to introduce intersections in complex cases.Besides, the efficiency of this approach is an issue since it is very time-consuming to obtain the solution field by computing partial differential equations.The whole-layer-inflation approach11can be more efficient than all the approaches mentioned above.It creates an outmost front first and ensures that this front has the same mesh topology as the initial front.The intermediate layers are created by interpolating appropriately between both fronts.Front checks are employed to ensure the validity of the outmost front and the techniques such as reducing the layer heights locally are necessary to remove the possible intersections.Difficulties exist in ensuring the robustness of these fixing techniques and avoiding the generation of badly-shaped elements.

Fig.1 Two typical front intersection cases.(a) Local intersections; (b) Global intersections.

By comparison with the above approaches, the layerby-layer ALM has its evident advantage.In this approach,the marching direction and distance can be particularly specified at each front node.As a result, the robustness of the algorithm and the quality of elements at complex corners can be ensured by employing those well-developed approaches.14,26–28However, the efficiency of intersection checks could be an issue in this approach since these computations involve a tremendous amount of geometric computations and are thus very time-consuming.If no cautions are taken,the procedure that treats front intersections can dominate the entire mesh generation workflow in terms of computing time.14For instance,it was reported that,for the baseline algorithm selected for comparison in Ref.,12more than 90%of the total time was consumed by this procedure, and for baseline algorithm in this paper, created by the author, this number is about 70%.

In this study, we propose a new approach that computes front intersections robustly and efficiently.It employs a pregeneration technique to create one whole layer at one step.The outermost facets of the new layer are classified into those belonging to the previous layers and those newly created in the present layer.Possible intersections introduced by the newly created fronts are then computed.After that, some of the pre-generated elements are removed at regions where intersections are identified, and a few pyramids are inserted to ensure that no quadrilateral facets are exposed.A set of new techniques and data structures are developed to improve the efficiency of the approach, which are summarized below.

(1) In octree intersection checking, it was suggested to change the tree topology at time when a new layer is created.29On the contrary, our suggestion is to create an octree that keeps its tree topology unchanged in the entire mesh generation workflow.This choice can minimize the computing cost taken by updating the tree structure.However, the imbalance of computing loads among different leaf octants is acceptable, since the boundary layer is a very thin structure attached to viscous walls.

(2) In the procedure of computing the intersection of two facets, a pre-check whether the outboxes of two facets overlap is necessary to speed up the computation.However, this check could be time-consuming because its calling times is huge.A new filter is thus developed to reduce this value.It projects all the facets belonging to a leaf octant onto an axis and then sorts these facets according to their projection results.After that, a small set of facets on which the intersection computations are really defined could be searched quickly.Meanwhile, a Principal Component Analysis(PCA)30based algorithm is suggested to compute the projection axis.

(3) In the modern computer system, a few levels of caches are configured apart from the main memory.Improving the hit ratio of caches is the key to improve the efficiency of an algorithm.The techniques mentioned above, i.e.,the adoption of an octree and the sorting procedure of the facets, are all helpful to implement a cache-friendly intersection check algorithm.Apart from these techniques,data structures are well designed so that the contiguously accessed data in each computing-intensive loop are stored in a contiguous space.

All the efforts mentioned above enable us to reduce the percentage of computing time taken by intersection checks to an acceptable level(below 30%for typical aerodynamics configurations).The developed approach creates a boundary layer mesh at a speed of 46 times faster than a commercial tool does.A few test cases with industry-level geometry complexity are selected to demonstrate the robustness of the developed approach.These results reveal that the classic layer-by-layer ALM can perform efficiently without sacrificing its performance in terms of robustness and element quality.

The remainder of this paper is organized as follows.Section 2 presents an outline of the proposed mesh generation approach,followed by a detailed description of the intersection check algorithm and the intersection removal algorithm in Section 3 and Section 4, respectively.Section 5 presents various numerical experiments that demonstrate the performance of the proposed approach.Section 6 concludes with the outcomes of the study and recommends some directions for future studies.

2.Overview of boundary layer mesh generation

Fig.2 presents the proposed workflow of boundary layer mesh generation.It inputs a surface mesh defining the problem domain and a few user parameters defining the preferred property of the output mesh.A typical set of these user parameters includes the height of the first layer, the ratio between heights of neighboring layers and the maximally allowed layer number.As illustrated in Fig.2, the workflow takes the following steps to create a qualified output:

Fig.2 Workflow of boundary layer mesh generation.

Step 1.Update active fronts.The input surface elements tagged as wall are initially set as the active fronts.In the middle of mesh generation, those facets defining the outskirt of the present boundary layer will be set as the active fronts.Accordingly, the nodes defining these active fronts are referred to as front nodes.

Step 2.Pre-generate one layer of elements.If the set of active fronts is empty, the workflow is ended by outputting the created boundary layer elements.Otherwise,by computing a marching direction and a marching distance at each front node, each front triangle will be stretched to a prism (see Fig.3(a)).When the fronts are smooth enough,we can luckily obtain a complete layer of prisms.Nevertheless, degenerated elements might be required in two circumstances.Firstly,multiple marching directions are required in complex corners in which one marching direction could not produce qualified or even valid prisms.In this circumstance, a few tetrahedrons are used to fill in the boundary layer (see Fig.3(b)).Secondly,a pre-generated prism has to be removed when its quality could not meet the quality standard or it induces a front intersection.In this circumstance, its neighboring prisms exposed their quadrilateral lateral faces to the next layer and these faces have to be transferred to triangular faces by filling in a few pyramids (see Fig.3(c)).

Step 3.Check front intersections.The pre-generated elements may intersect each other.These intersections usually occur at concave corners and small gaps.A check is thus necessary to investigate whether the intersections occur.If yes,the data such as which elements are involved in the intersections are collected; otherwise, the workflow switches to Step 1.

Step 4.Remove the intersections.In this step, some of those intersected elements are removed from the boundary layer so that the intersection induced by these elements could be removed.To maintain the difference of layer numbers between neighboring fronts under a predefined value, more elements might be removed in some cases.Accordingly, the set of pyramids, which use the lateral faces of prisms as bases and ensure all the new front faces are triangles, are updated.After intersection removal, the shape of the fronts are redefined and new intersections might be induced.Therefore,the workflow switches back to Step 3.The cycle between Steps 3 and 4 is ended when the new fronts are free of intersections.After that, the procedure of creating one layer of boundary layer elements is ended and the workflow switches back to Step 1.

Fig.4 shows the illustrating workflow of the steps.

Step 2 is a key to the success of the entire workflow, in which various techniques have been developed on how to compute marching directions and marching distances,how to configure multiple marching directions, to name a few.Interested readers can refer to Refs.18–20,28,31,32 for a comprehensive review on these techniques.Those techniques adopted in our implementation can be referred to in Ref.26.In summary,the prevailing techniques used in Step 2 are mainly based on local operations.As a result, the marching normal calculation part runs rather fast and occupies a small fraction of the total timing cost of mesh generation.In addition, Step 2 also includes some auxiliary substeps such as memory allocation and data preprocessing.In the experiment,this part will be collectively referred to as ‘‘Others”.

Fig.3 Three circumstances in which different types of elements are created in boundary layer.(a) A prism is created by extending a triangular front facet;(b)a tetrahedron is created at a convex corner where multiple marching directions are configured;(c)a pyramid is created to ensure no quadrilateral facets are exposed on the outermost front of boundary layer.

Fig.4 A 2D example illustrating the front intersection issue.(a)Step 1,activated fronts are marked as orange;(b)Step 2,dotted line and arrow compose a new layer of mesh element; (c) Step 3, intersecting mesh elements are marked as red; (d) Step 4, intersecting mesh elements have been removed.

A 2D illustrative example is given in Fig.5 Intersected front edges are marked as red.After 3 iterations, the mesh is intersection-free.Our focus in this study is on Steps 3 and 4.Step 3 involves tremendous geometric computations and its timing cost usually dominates the total timing cost of the presented workflow.Step 4 is a key to ensure the robustness of the presented workflow, in which detailed considerations must be contained to ensure the removal of intersections without sacrificing the completeness of boundary layers too much.In other words,as few as pre-generated elements should be removed in this step.

3.Front intersection check

3.1.Basic problem

The general front intersection check problem is defined as follows:

Problem 1.Given a set of triangular facets(denoted by F),find a set of facet pairs I= ｛｜｜f1,f2?F, and f1intersected f2}.If I is empty,we say F is free of self-intersections;otherwise,F is a self-intersected front.

Nevertheless, the set of front facets F could be subdivided into two subsets, namely Foldand Fnew.Here, Foldrefers to the group of facets that have been checked in the previous intersection computations, and thus, Foldshould be free of intersections.Fnewrefers to the set of facets that need further checks to judge whether it is free of self-intersections.Fig.6 presents two typical circumstances how to classify F into Foldand Fnewby the proposed approach.

(1) After pre-generating a new layer of boundary layer elements, the set front facets that belong to old fronts are classified as old, while the remaining ones are classified as new (see Fig.6(a)).

(2) If the pre-generated front is self-intersected, local updates are performed to remove these intersections.In this procedure, those facets that are added into the present front are classified as new, while the remaining ones are classified as old (see Fig.6(b)).

Fig.5 A 2D intersected front that needs three steps to remove its intersections.(a) Initial intersected front; (b) update front after removing two intersected front edges (marked as red); (c) update front after another removal process.(d) final front after three removal processes.

Fig.6 An example of classification of F in 2D.(a) Front pre-generation; (b) front after self-intersection removal.

Certainly F=Fnew∪Fold.Actually, if we consider IF, F as the intersection triangular facets in F, then IF,F=IFold,Fold∪IFold,Fnew∪IFnew,Fnew.We assume that the algorithm works and IFold,Fold=?because the intersection has been removed from previous check.By eliminating the double check, the general front intersection check problem is equivalent to a more specific problem, as defined below.

Problem 2.Given a set of front facets F that meet the assumption as follows:

(1) F could be subdivided into two non-overlapping subsets Foldand Fnew, and F =Fnew∪Fold

(2) Foldis free of self-intersection, find a set of facet pairs I =I1∪I2in which:

For real configurations,the majority of initial fronts cannot be extended completely.In other words, the generation of boundary layer elements usually stops earlier than expected.In this circumstance, the ratio between the magnitudes of Foldand Fnewincreases along with the layer-by-layer mesh generation procedure.It is thus a much cheaper choice to solve Problem 2 than to solve Problem 1.To be concise, we present the solution of Problem 1 in the following discussion,while the readers should keep in mind that the real implementation is targeted at solving Problem 2.The difference between the solutions of two problems is minor and will be mentioned only when it is necessary.

3.2.Overview of front intersection check scheme

To speed up the procedure of front intersection check, an octree is adopted in the proposed approach.Algorithm 1 illustrates the front intersection check scheme,which takes the following steps:

Step 1.Initialize the octree.If the front is the input surface mesh, the octree is initialized from a root box enclosing the mesh.At first, all the front facets are referred to by the root box.After that,the box is subdivided into 8 octants.For each facet referred to by the root box, if the outer box of the facet overlaps a child octant, the facet is referred to by the child octant as well.This subdivision is recursively performed on each child octant and stops when the number of facets referred to by the octant is fewer than 90 or the subdivision level exceeds 14.Those octants having no children are thus named leaf octants; accordingly, the ancestors of leaf octants are named branch octants.

Step 2.Update the octree.If the front is not the input surface mesh, i.e., it is the front after the generation of one or a few layers of boundary elements, the octree needs to be updated.Here, we have two choices.One choice is to rebuild the tree from scratch as Step 1,and the other choice is to keep the tree unchanged but only update the references of the front facets.Since the boundary layer is very thin and each layer of front facets scatter similarly as the input surface mesh to some extent, the second choice is adopted in the proposed scheme.The remaining issue is how to refer to the front facets.A fast and robust scheme will be proposed in Section 3.3.

Step 3.Compute intersected facets.In this step,each pair of facets that are referred to by the same leaf octant and possibly intersect each other are checked.If they do intersect each other, the pair is added into the set of intersected facet pairs,denoted by I.The returned value of the entire scheme is true when I is not empty, indicating that the front is selfintersected; otherwise, the returned value is false, indicating that the front is free of self-intersections.This step usually consumes the majority of time taken by Algorithm 1.In Section 3.4, a new technique will be developed to speed up this step.

Algorithm 1.Check Front Intersections (F,T,I).1.Input variables:2.F: the set of front facets 3 T: the octree tree 4.Output variable:5.I: the set of intersected facet pairs 6.If F is the input set of wall facets then 7.InitializeOctree(F, T)8.Else 9.UpdateOctree(F, T)10.End if 11.Return ComputeIntersectedFacets(F, T, I)

In both Steps 2 and 3,cache coherence is an important factor that determines efficiency.In Section 3.4, a few techniques that can improve cache coherence will be presented.

3.3.Tree update procedure

Assuming that the calling of after the procedure of pregenerating the ith layer of boundary elements, the input front is thus denoted by Fi.Similarly, the outmost front of the (i - 1)th layer of boundary elements is denoted by Fi-1.Denoting Fkeep=Fi∩Fi-1,Fdel=Fi-1｛Fkeep,Fadd=Fi｛Fkeepand the procedure of updating the octree thus includes two steps as shown in Algorithm 2:

(1) Remove the references of Fdelfrom the octree.

(2) Add the references of Faddinto the octree.

Alogrithm 2.UpdateOctree (F,T).1.Input variable:2.F: the set of front facets;3.Input & Output variable:4.T: the octree;5.Fold: the set of front facets presently referred by the octree.6.Fkeep =F ∩Fold 7.Fdel =Fold｛Fkeep 8.Fadd =F｛Fkeep 9.For f ?Fdel do 10.L : the set of leaf octants f refers to 11.For l ?L do 12.Rl: the set of facets referred to by l 13.Remove f from Rl 14.End for 15.End for 16.For f ?Fadd do 17.Let lguess be the store leaf octant that is highly likely to overlap f 18.SearchLeaves(T, f, lguess, L)19.For l ?L do 20.Rl: the set of facets referred to by l 21.Insert f from Rl 22.End for 23.End for

For each facet f ?Fdel, the leaf octants overlapping f are pre-known and explicitly stored in our implementation.Therefore, removing its reference from the octree is performed by iterating all such leaf octants and removing f from the reference list of these octants.For each facet f ?Fadd, the leaf octants overlapping f should be first searched and a reference is then added into the reference list of each of such octants.A brute-force implementation of the search procedure is performed by starting from the root octant and flagging its child octants that overlap f.This top-down procedure is recursively executed when all the leaf octants overlap f are flagged.Averagely, the time complexity of this searching algorithm is equal to the depth of the octree.It can slow down the tree-update procedure of updating the tree evidently because Faddis usually a large set in the real.A more complex but much faster searching algorithm is thus developed as follows:

(1) For each facet f ?Fadd, query its sister facet f′?Fi-1.Here, f and f′ are named sisters when they combines a prism of the present layer.In the case that f is not a top facet of a prism, e.g., it is the lateral facet of a prism or a pyramid, f′is the underside facets of mesh element.

(2) In our implementation, an octant that overlaps F is recorded in the previous computations.If this octant contains f′completely and the octant l′contains f completely (i.e., no overlapping exists between f and the boundary of l′, and here l′ means the leaf octant result of the lower layer of f),add a reference of f into the reference list of l;otherwise,perform Algorithm 3 to search all of octants overlapping f in a bottom-up fashion with l′ as the start point.This octant is recorded for use by the next layer.

Alogrithm 3.SearchLeaves (f, T, l, L).1.Input variables:2.T: the octree;3.f: a facet;4.l: a leaf octant that is highly likely to overlap f;5.Output variable:6.L: the set of leave octants overlapping f 7.t = l 8.While t does not contain f completely do 9.Let t be the father octant of itself 10.End while 11.While a leave octant c under the sub-tree rooted at t overlaps f do 12.Insert c into L 13.End while

3.4.Intersection computing procedure

The intersection computing procedure is performed by iterating all the leaf octants.For each pair of facets referred to by the same leaf octant, check whether the bounding boxes of the two facets are checked whether they overlap each other.If yes, a facet-facet intersection routine33is employed to confirm whether the two facets intersect each other in real.In this computing procedure, the bounding box check could be timeconsuming because its calling times is huge;as observed in real examples,this value is larger than the calling times of the facetfacet intersection routine by two or three orders.In this study,we recommend a new one-demensional filter to reduce this value.One-dimensional filter is widely used in mesh generation.The method was first proposed by Bonet and Peraire,29and one-dimensional filters can work well for spatial tree.34

This filter needs to project all the facets onto the longest axis of the leaf octant.To be concise, it is assumed that the longest axis of the leaf octant is the X axis in the following discussion.It is worth emphasizing that the presented algorithm,with minor changes, is applicable to the case that the longest axis is the Y or Z axis.

For each facet fi(i ?[1, n]) referred to by a leaf octant, its projection on the characteristic axis, for instance, X axis is a segment, is limited in the range between xi,minand xi,max(see Fig.8).Here, xi,minand xi,maxare the X coordinates of the left and right ending points of the segment.Thus, for a pair of facets (denoted by fiand fj), they possibly intersect each other when the following inequality holds:

To search all the pairs of facet that possibly intersect the facets sorted according to their projections along the axis, a two-step procedure is employed as shown below.

Based on the one-dimensional filter defined by inequality Eq.(1), Algorithm 4 presents a brute-force implementation to find all the pairs of facets that may intersect each other.In this study, a two-step sorting procedure is introduced to improve the efficiency of this one-dimensional filter.

Alogrithm 4.ComputeIntersected Facets (F,T).1.Input variable:2.T: the octree 3.Output variable:4.I: the set of intersected facet pairs 5.for all leaf octants modified by Algorithm 2 in T do 6.x : characteristic axis.7.Sort all the facet as Step 1 in Section 3.4.8.Modify the xi,max as Step 2 in section 3.4.9.for fi in leaf octant do 10.Find the lower/upper bound reference xlower,max ≥xi,min and xupper,max ≥xi,min by binary search.11.for j between lower and upper do 12.if BoxCheck(fi,fj) & TriangleTriangleOverlapTest(fi,fj) then 13.add to I.14.end if 15.end for 16.end for 17.end for

Fig.7(b) presents a graphic view of the sorted facets, in which the facet fiis represented by a point piwith coordinates xi,min,xi,max.After sorting,pi+1is always located at the top right of pi,if not considering special cases of xi+1,min=xi,minand/or xi+1,max=xi,max.Given a facet fi(i=1,2,???,n),we thus could define a shaded rectangle in Fig.7(b).All the facets that may intersect fi(including fiitself) are represented by the points located inside this rectangle.29Evidently, these facets combine a sub-array of the sorted facet array fi(i=1,2,???,n),and the start and end indices of this sub-array could be computed by a binary search.

A further improvement of this one-dimensional filter is to compute a better vector, instead of using the longest axis of the leaf octant directly,so that the number of facet pairs passing this filter could be reduced.Fig.8 presents an example illustrating why an appropriate selection of the projection vector could save intersection computations.As illustrated in Fig.8, if using vector x, the two facets pass the filter because their projections overlap each other.However, if using vector x′,the two facets fail to pass the filter because their projections are separated; as a result, a facet-facet intersection calling could be saved afterward.Since it takes a sizable amount of time to compute this projecting vector, we must achieve a tradeoff between the time taken by the computation and the time consumed by the facet-facet intersection callings that the algorithm saves.In this study, a Principal Component Analysis (PCA) based algorithm is suggested for its simplicity and efficiency,30which attempts to find the vector a that can minimize the objection function as shown below.

Here,aj(j=1,2,???,m=3n)refers to the projection of an ending point of vj(j = 1, 2, ???, m = 3n) on the vector a.The solution of this non-linear programming problem could be reduced to the maximum variance of aj, Algorithm 4 summarizes the computing procedure of intersections.Fig.9 shows the basic procedure of this calculation, and technical details that deserve more words are listed below.

(1) The PCA computation is performed after initializing the octree, and the computed vector a for each leaf octant will not be changed in the entire workflow of boundary layer mesh generation, again because the shape of the outermost mesh in each octant does not change greatly during the growth, especially the first few layers.We insist on this because the time cost of the PCA computation is sizable.

(2) After employing the one-dimensional filter, the classic bounding box filter is followed for those pairs of facets passing the one-dimensional filter.

3.5.Cache coherence

Algorithm 4 is the most time-consuming procedure, in which the data members of front facets are frequently accessed (see Lines 10–13 of Algorithm 4).The following guidelines are followed in our implementation to improve the cache hit ratio of these data access operations.

(1) All the front facets are stored in an array;

(2) The facets referred to by a leaf octant occupy a contiguous space in the array.

Meanwhile, the sorting procedure(see Line 7 of Algorithm 4)ensures the facets that possibly intersect a given facet occupy a contiguous space in the array, which could further improve the cache hit ratio of data access operations defined in Lines 10–15 of Algorithm 4.Algorithm 2 involves tremendous data access operations as well (see Lines 9–15 and Lines 16–23 of Algorithm 2).A simple sorting procedure is employed for Fdeland Fadd, which ensures the facets referred to by the same leaf octant are dealt with contiguously.In this way, the access on Rl(see Lines 13 and 20 of Algorithm 2) can achieve a very high cache hit ratio.In addition,a Hilbert sorting is introduced to the leaf octants so that two sub-arrays of facets referred to by two spatially close leaf octants are close enough in the array storing all the front facets.This sorting procedure could further improve the cache hit ratio of Lines 13 and 20 of Algorithm 2 in the case that more than one leaf octant overlap a facet.

Fig.7 Illustration of triangle projection.The red arrow presents the modification of Step 2.(a)Every triangle can be projected into an interval; (b) 2D representation of intervals.

Fig.8 Illustration of two triangles projected to two different axes.

4.Intersections removal

4.1.Pyramid conflict

Once we get from Algorithm 4,our task is to untangle all intersecting relationships in I by removing some of the current layer’s mesh elements.

Fig.9 Process of feature normal selection of DLR F6.

There may be conflicts in the removal process.As seen in Fig.10(a) in 3D, it is supposed that the green top facet exists in I, and the green prism is to be removed.After the removal,the two exposed sides quadrangle 1 and 2 (see Fig.10(b)) will pyramid with the underside triangle, but these two pyramids cannot coexist due to the topology conflict (see Fig.10(c)and (d)).Hence, each exposed triangle after element removal can only be adjacent to one prism at most.

4.2.Basic problem

For a particular intersecting pair , we call the intersecting relationship solved if either fior fjis removed.The whole intersection removal problem revolves around the removal of mesh elements.In general,our target is to preserve the mesh completeness while ensuring robustness.Here, mesh completeness means removing as few mesh elements as possible.However, this goal is difficult to achieve the global optimum.In addition, there is no need to design a very complex algorithm for this goal.Thus, in this section, the authors will seek a local optimum.

Mesh robustness can be summarized as follows:

(1) There should be no pyramid conflict when generating new pyramids.

(2) For ?Ik?I (0 ≤k

Fig.10 Illustration of pyramid conflict.(a)Initial mesh elements.Green facet exists in I,and prism below needs to be removed;(b)mesh after removal.Two quadrilaterals are exposed;(c)(d)two kinds of pyramids that can be generated on the exposed quadrilaterals.These two pyramids cannot coexist.

(1) For ?pi,pj?P′, there is no p ?P P′that satisfy p adjacent to both piand pj.

(2) For ??IP, not both of pk1,pk2?P′.

If we regard the base triangle of a prism as a graph node and adjacency as an edge and mark whether the prism is removed as true or false of the node, Problem 2 can be regarded as a global graph partition problem.What’s more,a facet may intersect multiple facets, i.e., the number of intersecting pairs in IPthat can be solved by removing different prisms varies, which makes this problem more difficult.In the author’s opinion,an optimal solution can only be obtained by brute force graph searching, but that can be timeconsuming.Actually,a suboptimal solution satisfying Criteria 1 and 2 is acceptable in the actual industry because we are more concerned about robustness.Therefore, here we utilize a local greedy algorithm to solve this problem for efficiency.

4.3.Intersection removal procedure

Our algorithm takes the removal of manifold around a mesh point as the smallest mesh removal unit.Here manifold means a set of facets sharing the same mesh point.Fig.11 shows the manifold removal process in 3D.If we want to remove the red point in Fig.11(a), all the mesh elements sharing the redpoint will be removed as shown in Fig.11(b), then fill pyramids on the exposed quadrilateral as shown in Fig.11(c).It is easy to prove that there must be no new pyramid conflict after manifold removal.If we take all the vertex of facet in I in the current layer as V= ｛v1,v2, ???,vvmax｝ into consideration,only part of V need to be removed to solve all the intersections in IP.The intersection removal procedure takes the following steps:

Step 1.Calculate object function.For each v ?V,an object value B(v) can be calculated from two aspects:

(1) The first index Q1(v ):Suppose A=the number of intersecting pairs in I is solved if v is removed, and B = the number of facet in the manifold of

(2) The second index Q2(v ): the largest equiangle skewness of prisms in manifold of v.

In our algorithm, B(v )=α*Q1(v )+Q2(v ), α is a userdefined scale, which is set to 100 by default.

Step 2.Removal manifold.Find viwith maximal B(vi),then remove the manifold of vi,and update B(vj)for vjneighbour to vi.Besides,solved intersecting facet pairs will be removed from I.Keep doing this removal and update process until I is empty.

Step 3.Fill the pyramids and update active fronts.The removal process possibly results in exposed quadrangles, so the necessary pyramids will be filled.Then update F with the calculation of Fnewand Foldfor further intersection checking,and only facts in Fnewdo intersection checking with F in the next loop in Fig.2.

It is worth noting that in the mesh with multi-normal, this removal operation in the first layer possibly results in the generation of zero volume elements, which need additional topology re-connections and node merge operations.The details will not be discussed in this paper.

Fig.11 Manifold removal process in 3D.(a) Initial outermost mesh.The red mesh point will be removed; (b) outermost mesh after removing all mesh elements that share this point; (c) mesh after filling pyramids.

5.Numerical experiments

The tests presented here were conducted on one personal computer (AMD R5 3600, 3.6 GHz, memory size 24 GB), all the following tests have been conducted in a single core of a CPU by a single thread, all breakdown of timing data were obtained by CPU sampling statistics (4000 sample/s), and small errors are possible.

5.1.Overall time performance

5.1.1.Overview of algorithm timing performance

To analyze the time performance of our algorithm, we adopt several typical surface mesh downloaded from the Internet:

(1) The DLR F6 aircraft model is utilized for illustrating the benefits in terms of computing time achieved by our algorithm,which is from the 2nd AIAA CFD Drag Prediction Workshop.

(2) The F16 aircraft model (referred to as F16 hereafter)obtained from GrabCAD35is introduced to prove the industrial-grade robustness of the external flow viscous-layer mesh generation.

(3) London tower models(from 23rd International Meshing Roundtable, IMR) are introduced to demonstrate the robustness of the proposed method for the large-scale surface mesh.

(4) The NASA Common Research Model(CRM)36DPW4,downloaded from the Fourth Drag Prediction Workshop (https://dpw.larc.nasa.gov/DPW4), is introduced to compare with our algorithm and Pointwise.

(5) The Delta and Widebody models, obtained through the software sumo (version 2.7.10, generated with default parameters),11are introduced to compare with our algorithm and Pointwise.

The surface mesh of the above models is enclosed in a box,but we do not show it for simplicity.Boundary layer mesh is generated from the surface mesh by default parameters,and then the mesh statistics, time statistics and peak memory during generation of the above models are presented in Table 1.

As shown in Table 1,the amount of time and memory consumed is basically positively correlated with the quantity of the viscous mesh.

5.1.2.Breakdown of meshing time

A breakdown of the proposed algorithm’s meshing timing data of viscous mesh generation of F6,F16,and London Tower are listed in Table 2.To illustrate the value of the algorithm in this paper in different parts of the procedure,we also compare the algorithm without any method proposed by this paper(baseline algorithm).

As we can see in Table 2, the whole Front Intersection Check process of intersection takes about 30%of viscous mesh generation process after our optimizations in the previous section,and the time-consuming Intersections Removal process is negligible.The vast majority of the program’s time is spent in Mesh Pre-generation, which includes normal calculation,memory allocation, and some of the start-up steps necessary for the ALM.

5.2.Analysis of timing performance

To prove our conclusion, we design several algorithms based on the original algorithm with small changes:

(1) Algorithm α: The algorithm will perform a tree reconstruction at each layer instead of performing Algorithm 3.

(2) Algorithm β: The algorithm will only perform a tree reconstruction once in the middle of the generation process.

(3) Algorithm γ: The tree will be dynamically updated at each layer, and the specific changes can be summarized as follows: (A) When the number of facets in a octant exceeds a certain amount,leaf octant is split to generate child octants, (B)When the number of facets of the sibling octants is too little, the sibling leaf node is merged into one octant.

(4) Algorithm δ: The 1D filter is disabled.

(5) Algorithm ε: The order of Fadd& Fdelis defined by Nature order, which means the order of Fadd& Fdeldefined by the input surface mesh F0, and it is partially ordered in geometry because surface mesh is generated patch to patch.This also means that the facets are simply spatially sorted according to face id of geometry.

(6) Algorithm ζ: The order of Fadd& Fdelis artificially randomized.

We chose the F16 models,which has plenty of complex corners and narrow gaps, as an example, and applied the above algorithms.Through the F16 mesh generated by the above algorithm, we can compare the detailed timing performance of each algorithm with our algorithm.

5.2.1.Static octree

The balance of octree will affect the efficiency of Algorithm 3 and Algorithm 4,but maintaining tree balance comes at a cost.We insist on keeping the tree topology unchanged in the entire mesh generation workflow because of the physical property of the boundary layer.

Table 1 Overview of our algorithm’s timing performance.

Table 3 shows the cost of maintaining dynamic octree under different algorithms and the time consumption of Front Intersection Check,and Fig.12 shows the specific time changes of each layer of Algorithm 4.It can be seen that the efficiency of static octree Front Intersection Check is slightly lower than that of other algorithms, but overall, the performance of each method on this function is not very different.In β’s curve, a small valley can be observed on Fig.12.This is because the whole tree is reconstructed there so that the tree is more balanced.However, this advantage cannot be maintained for a long time.After four more layers of growth, the gap between static octrees disappears.This proves that the dynamic octree in Algorithm α, β and γ is unnecessary, because it brings too little improvement.In addition,maintaining them would make the algorithm more complex and more expensive.

5.2.2.1D filter

The filter designed in this paper affects the efficiency of Algorithm 4.Though the PCA calculation is only at the first layer,it is also time-consuming because it involves eigenvector calculation.To prove the effectiveness, we calculate time consumed by Algorithm 4 at each layer (see Fig.13).

Fig.13 shows that the filter can filter out most candidate facet pairs to avoid the time-consuming BoxCheck function.The effect is very stable at each layer, and is even better when the tree is unbalanced(i >= 25).Compared with the filterless version,we can see in Table 4 that the 1D-filter can reduce the time consumption of Algorithm 4 by about half,while the cost is negligible.

Fig.12 Time consumption of Algorithm 4 under different algorithms at each layer.

5.2.3.Cache coherence

Timing data varying from different options are shown in Table 5.Although Hilbert sorting process does not improve the program performance much compared to origin order in this case, the order of the original surface mesh still faces uncertainties because it depends heavily on the surface mesh generation algorithm and geometric modeling design.

5.3.Comparison with other mesh generation engines

5.3.1.Local detail of viscous mesh

The surface mesh of the F6 model is presented in Fig.14(a),the mesh is enclosed in a box, but we do not show it for sim-plicity(the same for others).A complex corner node at the tail of the engine and several concave regions near the joint of the engine,and wing are presented in the close-up views in Fig.14(a).Fig.14(b) presents a cut-out view of the hybrid mesh of this model.A close-up view of the viscous-layer mesh around the above complex mesh point is presented in Fig.14(c) and Fig.14(d).

Table 3 Cost of maintaining dynamic octree under different algorithms and time consumption of front intersection check in F16.

Fig.13 Specific time consumption of Front Intersection Check by our algorithm and algorithm δ in F16 at each layer.

Table 4 Time consumption of front intersection check by our algorithm and algorithm δ in F16.(s).

Table 5 Comparison among different algorithms with different facets order of array in F16.(s).

Fig.15(a) presents the overall input surface mesh of F16, and Fig.15(b) shows cut-out views of the hybrid mesh(i.e., a front cut-out view and a side cut-out view).Fig.15(c)and Fig.15(d) show local details of the narrow gap between missile and airplane of our algorithm and Pointwise, mesh transition is smooth, trimmings are invoked around the proximate regions of the model to avoid potential collisions.Fig.15(e)and Fig.15(f)show the detail of landing gears of our algorithm and Pointwise.Fig.15(g) and Fig.15(h) show the boundary layer mesh near the narrow gaps.

Fig.16(a) shows the surface mesh of the London tower.The model has a large number of concave-convex corners and complex corners,and its surface mesh complexity is much higher than that of general aircraft, multiple normal methods are required on this model.Fig.16(b) shows the outcome outermost mesh and cut-view of tetrahedron mesh.

Fig.17 displays the model surface of the NASA Common Research Model (CRM),36Delta and Widebody models.

It is difficult to get a meaningful comparison to Pointwise in memory usage because Pointwise is a black box to the author and we cannot only execute T-Rex and see the difference.But by our observation of Pointwise by task management tool provided by Operation System,there are no significant differences between our algorithm and Pointwise.

Fig.14 Hybrid mesh of F6 model.(a)Surface mesh and complex corner;(b)cut-view of front perspective;(c)detailed hybrid mesh in the joint of hanging bins and fuselage; (d) another cut-view of joining corner.

Fig.15 Hybrid mesh of F16 model.(a) Front view and (b) Side view; (c) detail of boundary layer mesh in narrow gap generated by proposed algorithm(d)same perspective on boundary layer mesh generated by Pointwise;(e)detail of boundary layer mesh near landing gears generated by proposed algorithm; (f) same perspective on mesh generated by Pointwise; (e)(f) detail of boundary layer mesh near narrow gap generated by proposed algorithm.

We built a surface mesh module library with complex examples for realistic industrial relevance consideration, such as extremely sharp corners and concave corners,proximity surface mesh with narrow gaps, and anisotropic surface mesh.The library contains models that vary from simple geometric models,public airplanes,and missiles to internal flow engines,and the scale of its surface mesh ranges from hundreds to millions.After boundary layer mesh generation, the tetrahedral mesh is generated between the outermost facets of the boundary layer mesh and the bounding box by algorithm.37

5.3.2.Mesh quality

We compared our algorithm with the algorithm in Ref.12 and the sumo algorithm for hybrid mesh generation on typical aircraft geometries.Figs.18 and 19 show the comparison of our algorithm results with Ref 12 algorithm results for the CRM model and the comparison of our algorithm results with the sumo algorithm results for the Delta and Widebody models,respectively.It can be seen that our algorithm has visual advantages in terms of mesh completeness,mesh orthogonality when dealing with complex corner points.

A hybrid mesh of six models using Pointwise and two models using sumo with the same input surface mesh and parameters were generated, and Figs.20 and 21 present the equiangle skewness percentage comparison of each model respectively.For clarity, we use logarithmic axes to show the percentage.The equiangle skewness is computed as Pointwise(https://www.pointwise.com/doc/user-manual/examine/functions/equiangle-skewness.html),and equiangle skewness is represented as the maximum ratio of the cell’s included angle to the angle of an equilateral element.The angle skewness varies between 0 (good) and 1 (bad).It is recommended that this skewness measure be kept below 0.8 for a good mesh; values below 0.9 are acceptable, depending on the solver.In DLR F6, 0.280% of the prismatic elements generated by Pointwise had quality values over 0.8, whereas the ratio of the counterpart generated using the proposed method was only 0.232%.The data for the other three figures are different.The ≥0.9 part of our algorithm is slightly higher in F16, Tower, Delta,and Widebody.The two approaches are comparable in ≥0.8 part.Fig.21 presents the distribution of sumo and ours, and it is evident that prismatic mesh generated by sumo is lower compared to Pointwise and our algorithm in these two models.In addition, we stumbled upon the existence of negative volume prism in these two models in our experiments.Also, the average and worst skewness are provided in Table 6.

5.3.3.Time performance

Fig.19 A cut-view local details comparison of Delta and Widebody model mesh generated by our algorithm and sumo.

Fig.20 Mesh quality distributions of the prismatic elements of F6, F16, CRM, and tower compared to Pointwise.

Fig.21 Mesh quality distributions of the prismatic elements of Delta and Widebody compared to Pointwise and sumo.

To prove that robustness and completeness are not compromised, we made a comparison with the commercial software,and other research is presented in Table 7.11,12,19Then we generated the hybrid mesh of the models using Pointwise, a piece of prevalent commercial software for CFD mesh generation,with the same surface mesh.Pointwise uses T-Rex technology to generate the boundary-layer mesh: it first generated an anisotropic tetrahedral mesh around the model and then tried to form the boundary-layer mesh by combining three tetrahedral meshes into a prism.38In the T-Rex,insertion checking is processed by anisotropic tetrahedral mesh, and the algorithm utilizes a dynamic axis-aligned Binary Space Partition (BSP) tree for fast determination of a vertex’s geometrically neighboring faces and does intersection/proximity checking.39Here we compare our algorithm with the following three types of available boundary layer mesh generators:

(1) Pointwise V18.3, a widely used commercial software for CFD.All comparisons between the meshes generated by Pointwise and our algorithm in this paper were performed on the same machine with the same surface mesh and parameters.

Table 6 Equiangle skewness of viscous mesh of different mesh generation tools.

(2) The software sumo(version 2.7.10).Two models named as Delta and Widebody are obtained from the software for comparison,and then viscous mesh is generated with default parameters.11All comparisons of the two models were performed on the same machine with the same surface mesh and parameters.It is worth noting that the viscous mesh generation timing data of sumo is obtained by subtracting the tetrahedral mesh generation time from the overall time.

(3) The data from the published papers of other research.Compared with sumo of CRM, the mesh generated by sumo was run in a single core of an Intel Xeon X5650 CPU (Westmere architecture).The mesh generated by Wang Zhi12was run in a single core of Intel E5-1620.Besides, Wang Feng19′s mesh run in Intel E5-1630.

To show that this algorithm is stable and adaptable, a test model library containing 44 models was built.Public complex 3D models were selected considering realistic industrial rele-vance, such as missile, aircraft-8, bullet, and rocket.The scale of surface mesh in the library varies from thousands to millions, which can meet the applications of PC simulation.All the models have been generated with boundary-layer mesh until prisms are isotropic,and then isotropic Delaunay volume mesh generator is employed to fill the inner space.34Then we use Pointwise to generate hybrid mesh with the same input by T-Rex.Timing data of the above two methods are recorded,as shown in Fig.22.The horizontal axis represents the number of viscous units, and the vertical axis represents the time consumed to generate these viscous mesh (excluding Input/Output).Logarithmic coordinates are used for both axes for clarity.

Table 7 Mesh size and mesh generation time of different mesh generation tools.

Fig.22 Comparison of execution time between Pointwise and proposed algorithm (left) for the models(right) in test library.

An evident linearly related relationship is observed for the mesh scale and time cost, which means that our algorithm has strong commonality and expandability.Two linear fit lines are marked in Fig.22 by the least-squares method for comparison,presented as a blue and red solid line.The two lines show that algorithm performance is near the fitted line, and it takes about 44 s to generate ten million viscous-layer mesh, and the same number of Pointwise is 225.Therefore, our algorithm is faster than Pointwise by more than four times for the cases in the library.As presented,the numerical results provide overall satisfactory results of timing data.

In addition, the author also collected the worst quality of the grid generated by the proposed algorithm and the grid generated by Pointwise.The author found that all inputs can generate a positive volume grid,and in the example,the algorithm proposed by the author has the maximum equiangle skewness of 0.9997, and the value of the grid generated by Pointwise is 0.9994, which indicates that the quality of the grid generated by the algorithm in this paper is basically equivalent to that of Pointwise.

6.Conclusions

This paper describes the implementation of an intersection check scheme for viscous mesh generation in a threedimensional complex configuration.This proposed scheme uses layerwise intersection check to reduce internal intersection checking cost with optimized cache coherence and onedimensional filter.Then an intersection removal strategy is used to eliminate the intersection mesh.We built a model library formed by a public CAD model of industrial complexity to test the algorithm’s performance.For the cases in the library, the proposed scheme can generate ten million viscous meshes in one minute and is over four times faster than Pointwise.This technique significantly reduces the mesh generation time while maintaining the same robustness, mesh completeness,and quality as the classic ALM algorithm.Using the proposed scheme, we can say that intersection check in the ALM will no longer be the bottleneck of its efficiency.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This study was co-supported by the Zhejiang Provincial Science and Technology Program, China (No.2021C01108),the Innovative Research Foundation of Ship General Performance,China(No.14022105)and the Science and Technology on Scramjet Laboratory Fund, China (No.2022-JCJQ-LB-020-05).