A Feature-aware SPH for Isotropic Unstructured Mesh Generation
AA Feature-aware SPH for Isotropic Unstructured MeshGeneration
Zhe Ji a , Lin Fu b , Xiangyu Hu a , Nikolaus Adams a a Chair of Aerodynamics and Fluid Mechanics, Department of Mechanical Engineering,Technical University of Munich, 85748 Garching, Germany b Center for Turbulence Research, Stanford University, Stanford, CA 94305, USA
Abstract
In this paper, we present a feature-aware SPH method for the concurrentand automated isotropic unstructured mesh generation. Two additional ob-jectives are achieved with the proposed method compared to the originalSPH-based mesh generator (Fu et al., 2019). First, a feature boundary cor-rection term is introduced to address the issue of incomplete kernel supportat the boundary vicinity. The mesh generation of feature curves, feature sur-faces and volumes can be handled concurrently without explicitly followinga dimensional sequence. Second, a two-phase model is proposed to charac-terize the mesh-generation procedure by a feature-size-adaptation phase anda mesh-quality-optimization phase. By proposing a new error measurementcriterion and an adaptive control system with two sets of simulation param-eters, the objectives of faster feature-size adaptation and local mesh-qualityimprovement are merged into a consistent framework. The proposed methodis validated with a set of 2D and 3D numerical tests with different com-
Email addresses: [email protected] (Zhe Ji), [email protected] (Lin Fu), [email protected] (Xiangyu Hu), [email protected] (Nikolaus Adams)
Preprint submitted to XXX March 3, 2020 a r X i v : . [ c s . G R ] F e b lexities and scales. The results demonstrate that high-quality meshes aregenerated with a significant speedup of convergence. Keywords:
Unstructured Mesh, Smoothing Particle Hydrodynamics, DelaunayTriangulation, Particle Method
1. Introduction
Automated mesh generation is a critical and challenging topic for a widerange of scientific problems. In recent years, tremendous advancements havebeen made especially in unstructured triangular and tetrahedral mesh gener-ation methods [1]. The most well-established unstructured mesh generationmethods can be classified into five categories, i.e. advancing front methods(AFT) [2][3], Octree refinement-based methods [4][5], Delaunay refinement-based methods [6][7], Delaunay variational-based methods [8][9] and Particle-based methods [10][11][12]. In addition to the aforementioned methods, lotsof hybrids exist too especially in the context of parallel mesh generation [13].For a complete review of mesh generation methods, we refer to [13][14][15].There are generally two steps for most mesh generators to obtain the finalmesh, i.e. the initial mesh generation step and the mesh-quality improvingstep, providing the input geometry and the target feature-size function. Itis uncommon that any mesh generation package can achieve both the objec-tives of target feature-size distribution and desired mesh quality without someform of mesh-quality improvement procedure [14]. The mesh-improvementalgorithms mostly are different from the algorithms used for obtaining theinitial mesh. For example, the code CGALmesh [16] employs the restricted2elaunay refinement approach for the initial mesh generating. Then differ-ent types of optimization methods, e.g. Optimal Delaunay Triangulation[17] and vertex perturbation [18], can be applied or combined sequentially toimprove the mesh quality. Moreover, the computational time consumptionsof the two steps are different too. Due to the rapid development of parallelalgorithms, the current state-of-the-art mesh generators [19][20] can generatebillions of meshes within minutes. However, the mesh-improvement step canbe several orders of magnitude more expensive than the sequential meshingtime if a high threshold is set for mesh quality and if the geometry is compli-cated [21][22]. Therefore, the two-step procedure hinders the automation ofhigh-quality mesh generation. Meanwhile, developing the parallel version ofdifferent mesh-improving algorithms increases the workload and difficultiesfor software maintenance.Moreover, in order to generate an initial mesh, a certain sequence is re-quired for most mesh generators. For AFT, the initial front has to be meshedfirst before being advanced, e.g. in order to generate a volumetric mesh, thesurfaces of input geometry or sub-domain boundaries have to be meshed first[3]. For Delaunay refinement-based methods, a surface mesh is required tooas the input for conforming boundary tetrahedralization. Subsequently, thevolumetric mesh is generated with interior Delaunay refinement [23]. Thesame process applies to generating a surface mesh, i.e. the boundaries ofthe surfaces have to be meshed first [24]. For the particle-based method[25][12][10], the sampling and optimization of particle positions are also pro-ceeded in a hierarchical way. The quality of the mesh at feature boundarieshas significant influence on the resulting interior mesh [26]. For complex3D geometries, obtaining a high-quality surface mesh is non-trivial, whichresults in even more time for generating the volumetric mesh. In this sense,concurrent mesh generation independent of feature type and dimension forthe above-mentioned methods is not well developed. The concurrency is-sue is mitigated in a parallel environment, where the geometry can be di-vided into sub-domains and meshed separately [13][27]. In [20], a hybridtwo-level Locality-Aware Parallel Delaunay imaging-to-mesh conversion al-gorithm (LAPD) is developed to exploit the parallelism from both coarse-and fine-grained perspectives. Although partitioning the geometry into sub-domains allows for better concurrency, inside each local mesh generator thesame steps have to be executed in sequential.To tackle the aforementioned issues, several critical properties, e.g. in-herently suitable for parallel computing, independent of geometry features,consistent for initial mesh generation and mesh-quality optimization, guar-anteed convergence and etc., are required for the mesh generator. In thispaper, the particle-based mesh generation methods are focused to achievethe automated and concurrent unstructured mesh generation due to theirunique characteristics.The particle-based mesh generation methods share high similarity withthe Delaunay variational-based methods. Both approaches require a tar-get density/energy function [8][23][10]. The target mesh-vertex distribu-tion is calculated directly from the target density/energy function. Themesh quality is iteratively improved by applying different numerical schemesand by minimizing the energy or interpolation error. The key differencebetween the two approaches is whether the connectivity information is re-4uired during the optimization procedure [28]. For particle-based methods[10][11][29], the triangulation/tetrahedralization of the mesh vertices is onlyapplied once the mesh is generated, and the mesh quality is improved im-plicitly by pair-wise particle forces. While for Delaunay variational-basedmethods, re-triangulation/re-tetrahedralization is required every iteration todetermine a new position for each vertex [8]. Another appealing character-istic of the particle-based methods is that owing to the Lagrangian nature,the method is inherently suitable for large-scale parallel computing and canachieve scalable performance [30]. Moreover, due to the simplicity of themethod, it is easy to program and maintain with various parallel techniques,e.g. Message Passing Interface (MPI) [31], OpenMP [32] and CUDA [33].Several well-established particle-based codes are already available for variousarchitectures [34][35][36].Previously, we have developed an unstructured mesh generator based onLagrangian-particle fluid relaxation strategy [10], and further extended toadaptive anisotropic mesh generation [37]. The target-density function is de-fined on a multi-resolution Cartesian background mesh considering the effectsof distance to the geometry surface, curvature and singularity points. A set ofphysics-motivated governing equations has been then proposed and solved byadaptive-smoothing-length Smoothed Particle Hydrodynamics (SPH) [38].By introducing a tailored equation-of-state (EOS), the relative discrepancyof particle density and target density is characterized as pseudo pressure. Thepressure gradient results in pair-wise particle interaction force and drives par-ticles towards target density distribution while maintaining a regularized andisotropic distribution [28]. Later, the method is extended to parallel environ-5ent utilizing both MPI and Thread Building Blocks (TBB) [39][40] tech-niques. By introducing a repulsive surface tension model between distinctpartitioning sub-domains, the targets of domain decomposition, communica-tion volume optimization and high-quality unstructured mesh generation areachieved simultaneously within the same framework.However, the method still suffers from several problems. First, due to thefriction model, particle velocities are nullified every timestep to maintain nu-merical stability. Therefore, a large number of iterations is needed to achievea convergence [10]. Although this issue can be mitigated by starting froman initial particle configuration generated from another mesh generator, themethod then functions only as a mesh-improving algorithm. Second, sinceparticles at geometry boundaries are treated as boundary conditions of innerparticles and no special treatments are applied to remedy the lack of kernelsupport at the boundary vicinity, a well distributed particle configuration isrequired in advance to provide sufficient support. Therefore, a serial sequenceis required to generate the final mesh [10].In this paper, we propose an improved particle-based mesh generatorby developing a feature-aware SPH method to achieve two additional ob-jectives in the same framework. First, by introducing a feature boundarycorrection term to address the lack of kernel support, particle interactions atthe boundary vicinity are more consistent and accurate. If the target den-sity function is smooth, particle evolution on feature curves, feature surfacesand volumes can be proceeded simultaneously without explicitly followinga dimensional sequence. Second, a two-phase model is proposed to accel-erate convergence of the particle evolution. Due to the feature boundary6orrection term, the constraint of zero velocities at each timestep can besignificantly relaxed. Therefore the target of feature-size-adaptation can beachieved with much less iterations compared to the original algorithm. How-ever, the mesh quality can be affected by high-frequency acoustic waves ifmomentum accumulation is allowed. The original simulation setup is pre-ferred for final global mesh-quality improving. In this scenario, we proposeto characterize our mesh-generation procedure by a feature-size-adaptationphase and a mesh-quality-optimization phase. Different simulation param-eters are utilized for each phase. Consequently, the initial mesh generationand mesh-improvement steps are merged into a monolithic formulation.The paper is arranged as follows: In Section 2, the SPH-based isotropicmesh generator developed in [10] is first briefly reviewed. Then the pro-posed feature-aware SPH method is introduced in detail in Section 3. Thefeature definition, the correction term for feature boundaries and the two-phase mesh generation model are elaborated respectively. Initial particlesampling, particle stabilization strategy, flowchart and etc. are presentedin this section too. In Section 4, a set of validation tests are carried outconsidering both triangular surface-mesh and tetrahedral volumetric-meshgeneration with the presence of various sharp features (creases, sharp edgesand singularity points). Conclusion remarks in terms of the performance ofthe proposed method are given in the last section.
2. The SPH-based isotropic unstructured mesh generator
In this section, the SPH-based isotropic mesh generator developed in [10]is briefly reviewed. 7 .1. Geometry definition
The surface of the underlying geometry is represented by a zero level-setfunction [41] Γ = { ( x, y ) | φ ( x, y, t ) = 0 } . (1)The positive phase, i.e. Γ + = { ( x, y ) | φ ( x, y, t ) > } , is defined as the volu-metric mesh-generation region. The zero level-set, i.e. Γ = { ( x, y ) | φ ( x, y, t ) =0 } , is defined as the surface mesh-generation region. A multi-resolution block-structured Cartesian background mesh is employed to discretize the level-setfunction. Since the underlying mesh generation method falls into the category ofvariation-based approaches, a target feature-size function is required for theoptimization of both mesh-vertices distribution and mesh quality. In gen-eral, the target feature-size function h t can be defined considering arbitrarycharacteristic fields. In [10], the effects of distance to the geometry surface φ , the minimum distance from the geometry singularities ψ and the diffusedcurvature field in the domain κ are considered, h t = f ( φ, ψ, κ ) ,ρ t = g ( φ, ψ, κ ) , (2)where ρ t is the target density function defined for the relaxation procedure. ρ t can be obtained from h t following ρ t = h d , where d is the spatial dimension.The three characteristic fields contributing to the feature-size function arecalculated by solving three modeling equations, utilizing the same Cartesian8ackground mesh.Based on ρ t , the total mass for generating a volume mesh defined indomain Ω and the total mass for generating a surface mesh can be calculatedby integrating ρ t over the positive-phase of the geometry and the geometrysurface respectively. For details, we refer to [10]. The evolution of mesh vertices is computed based on a fluid relaxationprocess and is modeled as an isothermal compressible flow. The Lagrangianformation of the governing equations is dρdt = − ρ (cid:53) · v , (3) d v dt = − F p + F v , (4) d x dt = v , (5)where ρ denotes the density, v the velocity vector, x the position. F p and F v denote the pressure force and the viscous force respectively.To close the system, an EOS is required p = P ( ρρ t ) γ , (6)where P is a constant pressure, and γ > γ = 2 in the EOS, the momentumequation can be discretized as d v dt = − (cid:88) j m j (cid:16) p ρ t,i + p ρ t,j (cid:17) ∂W ( r ij , h ij ) ∂r ij e ij + (cid:88) j m j η i η j η i + η j (cid:16) ρ t,i + 1 ρ t,j (cid:17) ∂W ( r ij , h ij ) ∂r ij v ij r ij , (7)where h is the smoothing length and characterizes the interaction range of thekernel function, W ( r ij , h i ) is the kernel function, ∂W ( r ij ,h ij ) ∂r ij the derivative ofkernel, r ij = r i − r j the connecting vector between particle i and j , e ij = r ij r ij the unit vector of r ij , v ij = v i − v j , h ij = h i + h j the averaged smoothinglength of particle i and j , and η = ρν the dynamic viscosity. By setting ν ∼ . r c | v | , (8)where r c is the cut-off radius of particle interaction range, the local Reynoldsnumber of the simulation is always on the order of O (10). Meanwhile, inorder to maintain the stability of the simulation, the velocities of particlesare nullified to damp the kinetic energy after every timestep.The system is advanced in time using a simplified Velocity-Verlet scheme.The timestep size is calculated considering the CFL criterion, the viscouscriterion, and the body force criterion,∆ t = min (cid:16) . (cid:114) r c | a | , r c | v | , . r c ν (cid:17) , (9)where the artificial speed of sound is assumed as c s ∼ | v | max .10 . The feature-aware SPH method Following [10], the zero level-set function is utilized for geometry def-inition due to the flexibility of defining complex geometries. A Cartesianbackground mesh is utilized to discretize the level-set function. The size ofthe background mesh is defined according to the minimum and maximumtarget-feature-size. In the proposed method, we assume that the minimumand maximum mesh size, i.e. target-feature-size, is known before generatinga mesh. In order to maintain accuracy, we set the minimum grid size slightlysmaller than the minimum target-mesh size. In this paper, we generally en-sure that the minimum mesh size is 1.5 to 2 times of the background gridsize.For complex geometries, multiple sharp edges and singularity points mayexist, which cannot be resolved by a single level-set function. Additionalinputs are employed in this paper to define these features to facilitate thegeometry recovery. Following [10][28], the positions of singularity pointsare directly imported and singularity particles are generated accordingly.Moreover, sharp edges are defined with piecewise-linear B-splines, and eachsharp edge is assigned with a unique index. To cooperate with the feature-definition system, each segment of the curves is mapped onto the backgroundgrid and all the cells containing feature-curves are characterized as feature-cells. 11 .2. Feature definition
The same Cartesian background mesh used for the definition of the level-set function introduced in Section 2.1 is utilized. To characterize the features,each cell C i is assigned with a unique type and five categories of feature cellsare defined accordingly, i.e. positive cell ( C + ), negative cell ( C − ), feature-surface cell ( C s ), feature-curve cell ( C c ) and singularity cell ( C si ). Figure1 shows a simple example of the tag system with a 2D Zalesak’s disk. Inthis case, four different types of cells are defined. For 3D geometries, thefeature-curve cell can be defined if the geometry also contains sharp edges.For complex input geometries consisting of multiple sharp edges and sin-gularity points, the geometry surfaces/volumes may be split into several re-gions. Therefore, these features are further characterized by a unique index.For the 2D Zalesak’s disk (Figure 1), a total number of 9 features are defined,i.e. 4 feature-surfaces, 4 singularities and 1 positive volume. All indices arerepresented with a unique color. In the following, each type of feature cellis distinguished by a subscript k , i.e. C ∗ ,k . A lookup table is constructed toassociate the feature index with the feature type.12 igure 1: (a) Contour of φ = 0 for the Zalesak’s disk. (b) Tag system built for characterizing the geometryfeatures. Based on the above defined tag system, we can associate SPH particleswith feature types. Four types of particles, i.e. singularity particles ( P si ),feature-curve particles ( P c ), feature-surface particles ( P s ) and positive parti-cles ( P + ), are defined. In addition, based on the feature cell SPH particleslie in, the same feature index can be assigned too to each particle accordingto the cell feature index, and the same subscript k is used.To evaluate the total number of mesh vertices, the target density is inte-grated first in each feature region to obtain the total mass of each feature.For volume mesh, the total mass M v,k of feature volume V k is calculated by M v,k = (cid:90) V k ρ t dv, (10)The total mass for generating a surface mesh ( M s,k ) or line mesh ( M c,k ) can13e calculated similarly with M s,k = (cid:82) S k ρ t ds,M c,k = (cid:82) C k ρ t ds. (11)where S k and C k denote the feature surface and feature curve with index k respectively. The calculation procedure is completely parallelized since wecan evaluate the mass in each cell independently. Only a reduction operationis required at the end for cells of the same feature index. The total numberof mesh vertices, i.e. SPH particles, can then be calculated assuming eachparticle possesses unit mass following [10][28]. During the mesh-generation process, particles belonging to different fea-tures are not mutually independent. In the current paper, the interaction ofparticles between different features can be specified once the above-mentionedtag system is constructed. Similarly to [10][28], the interaction relationshipis characterized in a hierarchical way following the dimensional sequence, e.g. P c provide repulsive force for P s and P + to prevent penetration and P si areused as boundary conditions for all the other types of particles. Moreover,for particles with the same feature type but assigned with different featureindices, no interactions are applied and they are evolved independently.Fig. 2 shows the interaction relationship between P s (red circles) and P + (blue circles) in a 2D scenario. In this case, P s is used to impose the boundaryconditions for P + . For particle i (the highlighted blue circle), all red particleswithin the cutoff range (the dotted circle) are treated as boundary particlesto enforce the impermeability condition and to complete kernel support.14here are different approaches to model the interactions between “inte-rior” particles and “boundary” particles. Due to the meshless and Lagrangiannature of SPH, imposing of robust and accurate boundary conditions is dif-ficult [42]. In our previous work [10], no special treatment has been applied.In [37], additional ghost particles are generated and evolved dynamically atthe boundaries to prevent particles from penetrating the domain boundary inorder to generate adaptive anisotropic triangular meshes. Since more parti-cles are needed to achieve the desired mesh quality, additional computationalcosts and memory are required, which is critical in 3D scenarios. Figure 2: Interaction with boundary particles and illustration of inconsistency in kernel support region.
In this paper, we propose to address this issue with a new approach,utilizing the geometry information from the level-set function. First, werecall that with the SPH method, the gradient of a field quantity ψ as thefunction of the coordinate x can be calculated through a convolution integral15ver the cutoff range Ω of a kernel function W as (cid:53) ψ ( x ) = (cid:90) Ω (cid:53) ψ ( y ) W ( r , h ) d y . (12)A discrete approximation of Eq. 12 for a set of particles at positions r i is (cid:53) ψ i = (cid:88) j ( ψ i + ψ j ) (cid:53) i W ( r ij , h ij ) V j , (13)where the summation is performed over all the neighboring particle j withinthe cutoff range of particle i with weight W ( r ij , h ij ), and V j is the volume ofparticle j .Eq. 12 and Eq. 13 are consistent only when the cutoff range of the kernelfunction is fully inside the domain [43][42][44]. For the situation illustratedin Fig. 2, this full kernel support condition is no longer valid for P + near ge-ometry boundaries. To reproduce the kernel approximation at the boundaryvicinity, a renormalization term γ ( x ) can be introduced to Eq. 12 [42] [45] (cid:53) ψ ( x ) = 1 γ ( x ) (cid:90) Ω (cid:53) ψ ( y ) W ( r , h ) d y , (14)where γ ( x ) = (cid:90) Ω W ( r , h ) d y . (15)Upon integration by parts, one can obtain (cid:53) ψ ( x ) = − γ ( x ) (cid:90) Ω ψ ( y ) (cid:53) W ( r , h ) d y + 1 γ ( x ) (cid:90) ∂ Ω ψ ( y ) W ( r , h ) n dy, (16)where ∂ Ω is the boundary of Ω and n is the normal direction of ∂ Ω [42]16see Fig. 2). Several approaches can be applied to discretize the equationbased on how the geometry boundary is characterized, i.e. by particles ordiscrete segments [45]. In this paper, since we already use feature-surface,feature-curve and singularity particles to discretize the geometry boundaries,the fully discretized form of the gradient operator is employed, (cid:53) ψ i = − γ i (cid:88) j ( ψ i + ψ j ) (cid:53) i W ( r ij , h ij ) V j + 1 γ i (cid:88) b ( ψ i + ψ b ) W ( r ib , h ib ) n b A b , (17)where the first summation is calculated for all particles j with the sametype and feature index, while the second summation is performed for all“boundary” particles b from the perspective of particle i . n b and A b denotethe normal direction and “area” of particle b , respectively.The renormalization term γ ( x ) of particle i can be directly calculatedsimilar to the discrete Shepard coefficient [46] with the target volume ofparticle j , i.e. γ i = (cid:88) j W ( r ij , h i ) V t,j , (18)where the summation is carried out over all the neighboring particles of “i”,and V t,j = h dt,j .Finally, we can rewrite the pressure force term in Eq. 7 as F p,i = 1 γ i (cid:88) j m j (cid:16) p ρ t,i + p ρ t,j (cid:17) ∂W ( r ij , h ij ) ∂r ij e ij + 1 γ i (cid:88) b (cid:16) p ρ t,i + p ρ t,b (cid:17) W ( r ib , h ib ) n b A t,b , (19)where the target information A t,b = h d − t,j is used for simplicity.From a physical point of view, the extra terms introduced in the pressure17orce can be interpreted as a “repulsive force” from the boundary particles,which can prevent positive-phase particles from penetrating to negative do-main. Therefore, the system is more stable and robust. A regularized andconverged particle distribution at the geometry boundary is not necessarybefore evolving particles inside the domain. The extra term in the momen-tum equation dynamically adjusts particle positions at the boundary vicinityand maintains the validity of impermeability boundary conditions. Thus theconstraints of evolving the system following a hierarchical sequence can beeliminated, and less iterations are required for generating a mesh.Moreover, the implementation of the extra terms is straightforward. Allthe additional information can be directly calculated based on the level-setfunction and target feature-size function. The computational overhead of theproposed method is through the additional summation over all the neighbor-ing particles before calculating the forces. A considerably small increasein computation-time is introduced while additional memory requirement isnegligible. The performance will be assessed in Section 4. Based on the proposed method in the previous subsection, the mesh-vertex placement procedure can be simplified so that particles of all featuretypes can be evolved concurrently. In this section, a two-phase mesh gener-ation model is proposed to further improve convergence.The key concept is to merge the initial mesh generation phase and meshquality optimization phase into the same formulation with different set ofsimulation parameters. After initialization, particle velocities evolve in time18o allow for accumulation of momentum. During this initial phase, particlesare driven by relaxation forces and relax rapidly towards optimum positions(minimum residual). Owing to the proposed correction term in Section 3.3,the system can maintain stability towards equilibrium state. Following thisinitial phase, particles are re-adjusted locally to improve final mesh quality.Re-adjustment is achieved by switching back to the original method [10][28],in order to suppress high-frequency acoustic noises and deterioration of re-sulting mesh quality.It is worth noting that, because an equilibrium state has already beenachieved in the initial phase, particles are close to the target-feature-sizedistribution. Re-adjusting particle velocity at every iteration only locallyaffect the positions of mesh vertices. From the experiments in Section 4, thetarget of mesh-quality improvement can be achieved within a relative smallnumber of iterations.With the proposed model, the same governing equations are solved butwith different set of parameters. No intermediate triangulation/tetrahedralizationstep is required, and the implementation is trivial. In the following subsec-tions, technical details of the method are presented.
In order to find a suitable criterion to switch from the initial phase (PhaseOne) to the re-adjustment phase (Phase Two), a proper measurement of con-vergence to the targeted feature-size distribution is required. One possibilityis to evaluate convergence by monitoring the total kinetic energy. However,in complex geometries with singularities, the gradient of the target functionmay not be smooth, which contributes to inter particle forces. Consequently,19article velocities never vanish in those cases and a proper criterion basedon the kinetic energy is impractical.According to the definition of the EOS, upon reaching equilibrium, thesmoothing length of each particle equals the target feature size, i.e. h i = h t,i ,and remains constant. In [47], the authors proposed to approximate theparticle volume as (cid:101) V i = 1 (cid:80) j W ( r ij , h i ) , (20)i.e. the inverse of the particle specific volume. Ideally, when the particlefeature-size approaches the target, we can obtain (cid:101) V i → (cid:101) V t,i = 1 (cid:80) j W ( r ij , h t,j ) . (21)The calculation of (cid:101) V i is less sensitive to the gradient of the target function andthus is more suitable as the measurement of convergence. The normalizedform of the total “volume” in each feature is { V + ,k , V s,k , V c,k } = 1 V ref,k { i ∈ P + ,k (cid:88) i (cid:101) V i , i ∈ P s,k (cid:88) i (cid:101) V i , i ∈ P c,k (cid:88) i (cid:101) V i } , (22)where the notation “+”, “ s ”, “ c ” denotes the feature type, k the featureindex, V ref,k the reference volume of feature k , witch is an approximatedvalue of the total volume/area/length of the feature. We employ 1D/2D/3DSPH kernels for P c / P s / P + to calculate (cid:101) V i respectively.Based on Eq. 22, we can obtain the time history of V ∗ ,k . If the relativeerror ( E t ) between V ∗ ,k and the last time-averaged (over a certain period) V ∗ ,k is smaller then a threshold and the relative error ( E avg ) between two20onsecutive time-averaged V ∗ ,k also is smaller then a threshold, then we con-sider that the system has achieved the target of Phase One. The maximumvalue between E t and E avg is used as the measurements of the convergenceerror and is denoted as E sys . For all the cases in Section 4, we use an intervalof 200 timesteps for time averaging and V ∗ ,k is measured every 20 timesteps,i.e. 10 sampling points each time-averaging period. The threshold is set to5 × − . Although a correction term is introduced to address the issue of incom-plete kernel support and to relax the need to setting particle velocity to zeroeach timestep, we set particle velocities to zero periodically to guarantee thatparticle distributions are regularized and that timestep sizes do not becomevery small. Extensive tests suggest that a period of 50-200 timesteps gen-erally achieves a good balance between stability and convergence. In thispaper, we use a period of 100 for all the cases below. In addition, we monitortimestep-size variations, when we observe a decrease by an order of magni-tude compared to the previous step, particle velocities are set to zero. Insome extreme situations, e.g. ill-posed initial particle sampling, multiple sin-gularities, particle clustering during simulation, an additional damping termis introduced to the momentum equation during Phase One d v dt = − F p + F v − ε v , (23)where ε is the damping factor. In this paper, a range of ε ∈ [0 , .
2] issuggested. The damping term is not needed in Phase Two.21 .4.4. Remarks
In order to achieve the targeted feature-size adaptation and mesh-qualityoptimization with the same framework, a two-phase mesh generation modelis proposed in this section. First, the convergence of the feature-size adap-tation phase, i.e. Phase One, is monitored by introducing the normalizedconvolutional particle volume estimation term, i.e. V ∗ ,k . Once the conver-gence error of the simulation drops below a certain threshold, Phase One iscompleted. In Phase One, accumulation of particle momentum is allowed,which facilitates faster convergence to optimized mesh-vertex positions. InPhase Two, the original mesh generation scheme is employed, which sup-presses high-frequency acoustic waves and achieves improved mesh-quality.The key differences between Phase One and Phase Two are two param-eters, i.e. the period of velocity nullification and the damping factor. Thesuggested values/value ranges can be found in previous sections. Moreover,between Phase One and Phase Two, we use an intermediate period and anattenuation function to ramp down the parameters to ensure a smooth tran-sition of the computation. Particles are characterized according to feature types and feature indices.Forces between particles of different feature types are only one-sided. Forparticles of the same type but with different feature index, interactions aremutually exclusive. Consequently, we can evolve particles in groups, andeach group possesses particles of identical feature index. The timestep sizeof each group can be different. Based on this observation, we modify Eq. 922s ∆ t ∗ ,k = min i ∈ P ∗ ,k (cid:16) . (cid:114) r c,i | a i | , r c,i | v i | , . r c,i ν (cid:17) , (24)where the notation “ ∗ ” represents different feature types.Based on the smooth target feature-size function assumption, in simplegeometries the timestep sizes among all features should be similar. However,for complex geometries and large resolution variations, timestep size distri-bution across all feature types may deteriorate the convergence of the fullsystem. By using feature-aware timestep sizes, asynchronous convergence ofeach feature improves robustness. Moreover, evolving particles in groups alsofacilitates pipeline parallelization of the code. In the previous particle-based mesh generation methods [10][28], the ini-tial mesh vertices, i.e. SPH particles, are generated randomly with constantprobability inside the computational domain. This approach demonstratesthat the proposed methods can generate a mesh from extreme poor initialconditions. However, the iteration number required for obtaining a high-quality meshes from this initial condition is very high, thus making it im-practical in realistic applications and much more expensive than other meshgenerators [15][28].Instead of using constant probability function, the initial particle distri-bution is generated randomly with probabilities proportional to the targetdensity field. We use the Mersenne Twister algorithm [48] to calculate ran-dom numbers, which is of O (1) time complexity. Therefore, no significantextra computational overhead is introduced. The initial sampling is closer to23he target feature-size distribution, and the number of iterations to achievean optimized mesh quality can be reduced. A detailed flowchart of the proposed method is given as Algorithm 1.The method is implemented based on the parallel framework for adaptive-resolution SPH developed in [30][49]. Since the equations of the proposedalgorithm exhibit high similarities with standard SPH for gas dynamics, itis straightforward to implement the mesh generator from an existing SPHcode. For more technical details we refer to [30][50].24 lgorithm 1
Flowchart of the feature-aware mesh generation method
1: Initialize the background Cartesian mesh and the level-set function (Eq. 1);2: Load user-defined feature curves and singularity points;3: Calculate the target density function (Eq. 2) based on the background mesh;4: Define the feature-based tag system for each cell ( C i ) and calculate the target infor-mation for mesh generation (Eq. 10 and Eq. 11);5: Allocate particles based on the target information and initial particle sampling;6: while t < t end do
7: Update the multi-level data structure; (cid:46)
See [30] for detailed description for k = 0 → K − do (cid:46) K is the total number of features if V ∗ ,k is not converged then
10: Define ρ t , h t , n and etc. for P ∗ ,k from background mesh;11: Construct Verlet neighbor list; (cid:46) See [30] for detailed description
12: Reset particle forces;13: Reset particle velocities; (cid:46)
See Section 3.4.3 for detailed description onwhen to do the velocity reinitialization
14: Calculate the renormalization term γ i (Eq. 18)15: Calculate pressure force F p (Eq. 19);16: Map F p to feature k for P ∗ ,k ∈ ( P s ∪ P c ) (Eq. 30 and Eq. 31 in [28]);17: Calculate ∆ t ∗ ,k (Eq. 24);18: Update the mid-point velocity (cid:101) v n + (Eq. 26 in [28]);19: Accumulate viscous force F v (Eq. 7);20: Map F v to feature k for P ∗ ,k ∈ ( P s ∪ P c ) (Eq. 30 and Eq. 31 in [28]);21: Calculate ∆ t ∗ ,k (Eq. 24);22: Update predicted velocity v n + (Eq. 27 in [28]);23: Update particle position r n +1 (Eq. 28 in [28]);24: Map r n +1 to feature k for P ∗ ,k ∈ ( P s ∪ P c ); (cid:46) See Section 3.5 in [28]
25: Calculate V ∗ ,k and check convergence;26: end if end for if Do post-processing then
29: Generate the corresponding mesh and calculate mesh quality;30: end if end while . Validation and demonstration In this section, a set of numerical experiments varying from generating2D triangular to 3D surface/volume meshes is considered. The performanceof our method is assessed, and the results are compared with previous workand other established methods. In all the experiments below, we refer to theresults obtained by the proposed method as “ improved ”. All cases in thissection are computed on a desktop workstation with Intel R (cid:13) Xeon R (cid:13) CPU(E5-2630V2, 2.60GHz, 12 threads).For 2D triangulation, the mesh is constructed following [10] utilizing alocal Voronoi tessellation. For 3D tetrahedralization, the open-source codeTetGen [51] is used. Flip operations included in TetGen (2-3 flip, 3-2 flipand 4-4 flip) are invoked to improve connectivity. It is worth noting thatthe mesh quality is only measured for comparison purposes. In the pro-posed method, the mesh quality is optimized implicitly and the triangula-tion/tetrahedralization step is not invoked in the mesh-generation procedure.The quality of isotropic triangular meshes is quantified by G = 2 √ SP H and the angle θ , where S is the triangle area, P the half-perimeter and H the length of the longest edge. θ max , G avg , G min / θ min and θ min are themaximum, average, minimum values and the averaged minimum angle ineach triangle respectively. θ < is the number of triangle with angle smallerthan 30 ◦ . The quality of isotropic tetrahedral meshes is quantified by thedihedral angle θ and radius ratio γ = 3 r in r circ , where r in is the inradius and r circ the circumradius of a tetrahedron. θ max , γ avg , γ min / θ min and θ min arethe maximum, average, minimum values and the averaged minimum dihedralangle in each tetrahedron respectively. The number of slivers are measured26y counting the number of tetrahedra with different smallest dihedral angles,i.e. 10 ◦ , 20 ◦ , 30 ◦ and 40 ◦ . First we consider a square with feature-size function h t ( x, y ) = h max − h min √ (cid:112) ( x − + ( y − + h min , (25)where ( x, y ) ∈ [0 , h max and h min ) are 4 .
88 and 0 .
244 respectively. The total number of particles is 2524.The results are compared with the previously proposed SPH-based meshgenerator [10]. A variant of the method [10] is considered (referred to as[10.v2]), where particle velocities are set to zero with same period, i.e. every100 timestep, as with the current method. All computations start with thesame initial particle distribution (see Fig. 3 (a)).The convergence history comparisons of E sys , θ min , G avg and G min areshown vs. the total wall-clock runtime in Fig. 3. The angle distribution ismeasured and compared in Fig. 3. The dashed lines in Fig. 3 (b) indicatethat the system has achieved the prescribed error threshold of Phase One.First it is observed that [10.v2] fails to converge in the plotted time range,and the mesh-quality is the worst. Due to the absence of the feature-awarecorrection term proposed in this paper, the system is not able to maintainstability at the vicinity of the boundary when the accumulation of particlemomentum is allowed. Successive particle losses at the domain boundary areobserved. Remapping particle into the positive phase are required to handlethe particle lost issue. Second, the “ improved ” method achieves the target-27eature distribution significantly faster than the previous method [10]. Theruntime for both methods to achieve the target of Phase One are approxi-mately 4 . s and 34 s , which implies a speedup factor of 8 .
1. This can also bedemonstrated by the particle distributions at different iterations (see Fig. 4).The quality of the generated meshes exhibits no obvious difference (see Fig.3 (c-f) and Table 1), which can also be seen in the Delaunay triangulationsof the final mesh in Fig. 5.
Runtime (s) E rr o r
10 20 30 40 5010 -8 -7 -6 -5 -4 -3 -2 -1 Fu et al.ImprovedFu et al. v2
Runtime (s) θ m i n Fu et al.ImprovedFu et al. v2 (c)(d) (e) (f)(a) (b)
Runtime (s) G av g
10 20 30 40 500.60.70.80.91
Fu et al.ImprovedFu et al. v2
Runtime (s) G m i n Fu et al.ImprovedFu et al. v2
Angle P e r ce n t a g e ( % )
50 100 150051015202530354045
Fu et al.ImprovedFu et al. v2
Figure 3: Square: (a) Initial particle distribution. The convergence histories of (b) E sys , (c) θ min , (d) G avg and (e) G min . (f) Histogram of the angle distribution. igure 4: Square: Particle distributions at (a) iteration 3000, (b) iteration 6000 and (c) iteration 12000.Upper row: results calculated by the algorithm developed by Fu et al. [10]. Bottom row: results calculatedby the proposed algorithm.Figure 5: Square: Delaunay triangulation of the final mesh and the zoomed-in view at iteration 20000.(a) Mesh generated by the algorithm developed by Fu et al. [10]. (b) Mesh generated by the proposedalgorithm. able 1: Mesh quality of the Square case G avg G min θ max θ min θ min θ < N tri runtime N iter Fu et al.
Improved
Next, we consider the 2D Zalesak’s disk following [10]. The geometryconsists of a slotted circle of radius 15. The width and length of the slot are5 and 25 respectively. The mesh generation region is inside the slotted circle.The definition of feature-size function and h min / h max are identical to [10].Instead of generating the initial sampling of particles following the targetdensity function ρ t ( x ), we use a uniform probability function to sample theparticles (see Fig. 6 (a)), which is the same with [10].Based on the convergence histories of E sys , θ min , G avg and G min in Fig6 (b-c), a high-quality mesh is obtained and both the targets of Phase Oneand Phase Two are achieved. The system remains stable during the mesh-generation procedure, which is demonstrated by the particle distributionsat different iterations (see Fig. 7). The mesh quality calculated by theproposed method is close to the result from [10], and is compared in Table2. The main differences of the two results are the iterations required toachieve the convergence. Since in [10], there is no explicit measurement ofconvergence, we consider the instant as converged when all the mesh-qualitycurves reach the plateau. The method [10] achieves convergence after ca. N tri the total number of triangles runtime denotes the total wall-clock time in seconds N iter denotes the number of iterations for the measured runtime All the mesh quality in this table are measured at N iter , and remains consistenthereafter. Iteration E rr o r -9 -8 -7 -6 -5 -4 -3 -2 -1 Iteration G av g Iteration G m i n Iteration θ m i n (a) (b) (c)(d) (e) (f) Angle P e r ce n t a g e ( % )
50 100 150051015202530
Figure 6: Zalesak’s disk:(a) Initial particle distribution. The convergence histories of (b) E sys , (c) θ min ,(d) G avg , (e) G min . (f) Histogram of the angle distribution. igure 7: Zalesak’s disk: Particle distributions at (a) iteration 3000, (b) iteration 6000 and (c) iteration12000. (d)(e)(f) Delaunay triangulation of the final mesh and the zoomed-in views at iteration 20000.Table 2: Mesh quality of the Zalesak’s disk case G avg G min θ max θ min θ min θ < N tri runtime N iter Improved
Fu et al. ∼ ∼ In order to evaluate the proposed method in more complex scenarios,we consider the geometry Tyrannosaurus rex, referred as “T-Rex”. Thetarget-feature size function is defined on the geometry surface based on thesmoothed curvature field. We set h min = 0 . h max = 0 .
01, which resultsin a total number of 21,243 particles. In this case, only surface meshes aregenerated, i.e. only P s are considered in the simulation. The initial particlesampling is presented in Fig. 8 (a) and colored with the target feature size.32he convergence histories of E sys and θ min (Fig. 8 (b) and (c)) suggestthat both methods achieve the convergence successfully. The “ improved ”method achieves the target of Phase One after ca. 100s, while method [10]takes ca. 400s (see Fig 8 (b)), which results in a speedup factor of 4. Fromthe triangulated results in Fig. 9, highly similar meshes are generated byboth methods. A slightly better mesh quality of G min , θ max and θ min areobserved as shown in Table 3. Table 3: Mesh quality of the T-Rex case G avg G min θ max θ min θ min θ < N tri runtime N iter Fu et al.
Improved untime (s) E rr o r
200 400 600 800 100010 -8 -7 -6 -5 -4 -3 -2 -1 Fu et al.Improved
Runtime (s) θ m i n Fu et al.Improved
Angle P e r ce n t a g e ( % )
50 100 150051015202530354045
Fu et al.Improved (c) (d)(a) (b)
Figure 8: T-Rex:(a) Initial particle distribution. The convergence histories of (b) E sys , (c) θ min . (d)Histogram of the angle distribution. igure 9: T-Rex: (a) Delaunay triangulation of the final mesh and (b)(c)(d) the zoomed-in views atdifferent camera positions. Upper row: results calculated by the algorithm developed by Fu et al. [10].Bottom row: results calculated by the proposed algorithm. In this section, the proposed method is extended to generate 3D tetrahe-dral meshes. A sphere with two target feature-size distributions are consid-ered, which are referred as “Sphere01” and “Sphere02” respectively.For “Sphere01” case, we choose h t ( x, y ) = h max − h min R (cid:112) ( x − x ) + ( y − y ) + ( z − z ) + h min , (26)where R = 45 is the radius and ( x , y , z ) is a point on the sphere sur-face. The minimum feature-size h min = 0 .
375 and the maximum feature-size h max = 7 .
5. The total number of particles is 12,540. The initial particlesampling is plotted in Fig. 10 (a).From the time history of E sys (see Fig. 10 (b)), the convergence of Phase35ne is achieved after ca. 430s for the “ improved ” method and 1750s formethod [10]. The “ improved ” method reaches the optimized mesh qualityremarkably faster than method [10], which is revealed from the convergencehistories of θ min (Fig. 10 (c)) and θ < / θ < / θ < . (Fig. 10 (d)). Thehistograms of dihedral angle and radius ratio distribution (Fig. 11) showthat the underlying measured statistics remain almost constant after 4000iterations for the “ improved ” method. At 20000 iterations, both methodsachieve approximately the same distributions. From the snapshots of thesimulation at different instants (see Fig. 12), the same conclusions can bemade. Except for the improvement achieved in runtime, the “ improved ”method also exhibits slightly better mesh quality, which can be seen fromTable 4. Table 4: Mesh quality of the Sphere01 case θ min /θ max γ min /γ avg θ min θ < θ < θ < θ < N tet runtime N iter Fu et al.
Improved untime (s) E rr o r
500 1000 1500 200010 -9 -8 -7 -6 -5 -4 -3 -2 -1 Fu et al.Improved
Runtime (s) av g
500 1000 1500 200048505254565860
Fu et al.Improved
Runtime (s) N u m b e r o f T e t r a h e d r a l s
500 1000 1500 200010 <20<20<30<30<40<40 (a) (b)(c) (d) Figure 10: Sphere01:(a) Initial particle distribution. The convergence histories of (b) E sys , (c) θ min , and(d) θ < , θ < and θ < . ll Dihedral Angle P e r ce n t a g e ( % )
50 100 150051015
Fu et al.Improved
All Dihedral Angle P e r ce n t a g e ( % )
50 100 150051015
Fu et al.Improved
All Dihedral Angle P e r ce n t a g e ( % )
50 100 150051015
Fu et al.Improved
Radius Ratio P e r ce n t a g e ( % ) Fu et al.Improved
Radius Ratio P e r ce n t a g e ( % ) Fu et al.Improved
Radius Ratio P e r ce n t a g e ( % ) Fu et al.Improved (a) (b) (c)(d) (e) (f)
Figure 11: Sphere01: Histogram of the dihedral angle distribution at (a) iteration 4000, (b) iteration 8000,and (c) iteration 20000. Histogram of the radius ratio distribution at (d) iteration 4000, (e) iteration 8000,and (f) iteration 20000.Figure 12: Sphere01: Particle distributions at (a) iteration 1000, (b) iteration 4000, (c) iteration 8000and (d) iteration 20000. Upper row: results calculated by the algorithm developed by Fu et al. [10].Bottom row: results calculated by the proposed algorithm. Only particles belonging to positive cell ( C + )are plotted in (b)(c)(d) and colored by target mesh-size. igure 13: Sphere01: Tetrahedralization of the final mesh. (a) The full-size view. (b) The zoomed-inview. (c)(d) Clipped views. Upper row: results calculated by the algorithm developed by Fu et al. [10].Bottom row: results calculated by the proposed algorithm. For the case “Sphere02”, we follow [8]. The target feature-size functionis defined as h t ( x, y ) = 0 .
025 + 0 . | (cid:112) x + y + z − . | . (27)The minimum feature-size h min = 0 .
025 and the maximum feature-size h max =0 . GSM ”), the GSM coupled with a particle-based method [11] (denoted as “ particle+GSM ”) and the particle-basedmethod (denoted as “ particle ”).Fig. 14 (a) illustrates the time history of E sys . The target of Phase Oneachieves after ca. 240s. The dihedral angle and the radius ratio distribu-tion histograms from the “ improved ” method, the “ GSM ” method and the“ particle+GSM ” method are plotted together in Fig. 14 (b) and (c) respec-tively. It can be observed that the “ improved ” method generates the largest39roportion of regular tetrahedrons than the other two methods. For the num-ber of “slivers” generated, the “ improved ” method outperforms the “ parti-cle ” method and the “
GSM ” method (see Table 5). The “ particle+GSM ”presents the best mesh quality overall. However, in terms of computationalcost, the “ improved ” method only takes 236.42s to achieve the mesh quality,while 8,787.22s are required for the “ particle+GSM ” method (see Table 5).With slightly more computational effort, the proposed method achieves aspeedup factor of 37.
Table 5: Mesh quality of the Sphere02 case θ min /θ max γ min /γ avg θ < θ < θ < θ < N tet runtime N iter Improved
GSM particle+GSM particle The runtime of the “
GSM ”, “ particle+GSM ” and “ particle ” method are measured byan Intel R (cid:13) Xeon R (cid:13) CPU(E5645, 2.40GHz, 12 threads). igure 14: Sphere02: (a) The convergence histories of E sys . Comparison of (b) dihedral angle distributionhistogram and (c) radius ratio distribution histogram between the results in [8] and result obtained bythe proposed algorithm. (d) Tetrahedralization of the final mesh. (e) Clipped view of the final mesh. (f)Tetrahedrons with dihedral angle below 30 ◦ . Next, we consider a cube. In this case, multiple sharp edges and singu-larity points are presented. The target feature-size function is defined as h t ( x, y ) = h max − h min √ (cid:112) ( x − + ( y − + ( z − + h min , (28)where h min = 0 . h max = 4 .
88, and ( x, y, z ) ∈ [0 , improved ”method outperforms method [10] slightly in each category. Table 6: Mesh quality of the Cube case θ min /θ max γ min /γ avg θ min θ < θ < θ < θ < N tet runtime N iter Fu et al.
Improved E sys , (c) θ < , θ < and θ < , and (d) γ avg . e) (f) (g) (h) All Dihedral Angle P e r ce n t a g e ( % )
50 100 15005101520
Fu et al.Improved
All Dihedral Angle P e r ce n t a g e ( % )
50 100 15005101520
Fu et al.Improved
All Dihedral Angle P e r ce n t a g e ( % )
50 100 15005101520
Fu et al.Improved
All Dihedral Angle P e r ce n t a g e ( % )
50 100 15005101520
Fu et al.Improved
Radius Ratio P e r ce n t a g e ( % ) Fu et al.Improved
Radius Ratio P e r ce n t a g e ( % ) Fu et al.Improved
Radius Ratio P e r ce n t a g e ( % ) Fu et al.Improved
Radius Ratio P e r ce n t a g e ( % ) Fu et al.Improved (a) (b) (c) (d)
Figure 16: Cube: Histogram of the dihedral angle distribution at (a) iteration 7000, (b) iteration 10000,(c) iteration 20000, and (d) iteration 40000. Histogram of the radius ratio distribution at (e) iteration7000, (f) iteration 10000, (g) iteration 20000, and (h) iteration 40000.
The convergence curve of E sys in Fig. 15 (a) reveals that, benefitingfrom the correction term and the accumulation of particle momentum, the“ improved ” method converges significantly faster than method [10]. Fromthe snapshots of particle distributions at various instants of the simulation(see Fig. 17), we can see that at iteration No. 900, the equilibrium state ofall P s and P c has been achieved for the “ improved ” method. For method [10],only convergence of P c is achieved, since the particles are evolved following adimensional sequence. The same observation can be made for iteration No.2200, where P s converge to the optimized positions while P + are remainedin the original positions. Conversely, for the “ improved ” method, particlesof all features are evolved together. This can be observed from Fig. 17 (b)and (c). Consequently a faster convergence is achieved.43 igure 17: Cube: Particle distributions at (a) iteration 900, (b)(c) iteration 2200, (d) iteration 6400,(e) iteration 10000, and (f) iteration 40000. In each sub-figure, left: result calculated by the proposedalgorithm; right: result calculated by the algorithm developed by Fu et al. [10]. Only particles belongingto positive cell ( C + ) are plotted in (c)(d)(e)(f) and colored by target mesh-size. igure 18: Cube: Tetrahedralization of the final mesh. (a) The full-size view. (b) The zoomed-in view.(c) Clipped view. Upper row: results calculated by the algorithm developed by Fu et al. [10]. Bottomrow: results calculated by the proposed algorithm. Lastly a complex geometry of spur gear developed for gear lubricationtests [52] is considered. The size of computational domain is [86 . × . × . h min = 0 . h max = 3 .
2. The resulting total number of particles are 42,647 and 169,096respectively. The objective of this case is to test the robustness of the pro-posed feature-aware SPH in the mesh generation of complex geometries.The initial condition of “Gear01” is plotted in Fig. 19 (a). From the time45istories of θ < , G avg , and G min (see Fig. 19 (b-d)), both methods achievesan equilibrium state. Similar mesh qualities are observed for both methodsand 78% reduction in wall-clock time is accomplished by the proposed method(see table 7). From the particles distributions at different instants (Fig.20) and the final triangulated mesh (Fig. 21), the “ improved ” method canmaintain a stable simulation even when the accumulation of momentum isallowed and the particles of all features are evolved together. The geometryis fully recovered. 46 untime (s) G av g
100 200 300 400 500 6000.80.850.90.95
ImprovedFu et al.
Runtime (s) G m i n
100 200 300 400 500 6000.20.30.40.50.6
ImprovedFu et al.
Runtime (s) N u m b e r o f t r i a ng l es
100 200 300 400 500 60010 ImprovedFu et al. (a) (b)(c) (d)
Figure 19: Gear01:(a) Initial particle distribution. The convergence histories of (b) θ < , (c) G avg , and(d) G min . igure 20: Gear01: Particle distributions at (a) iteration 500, (b) iteration 1000, (c) iteration 2000 and(d) iteration 10000. Upper row: results calculated by the algorithm developed by Fu et al. [10]. Bottomrow: results calculated by the proposed algorithm.Figure 21: Gear01: Delaunay triangulation of the final mesh. (a)(b) The full-sized views. (c)(d) Thezoomed-in views at different camera positions. Upper row: results calculated by the algorithm developedby Fu et al. [10]. Bottom row: results calculated by the proposed algorithm.Table 7: Mesh quality of the Gear01 case G avg G min θ max θ min θ min θ < N tri runtime N iter Fu et al.
Improved
Figure 22: Gear02: (a) Initial particle distribution. (b) Tetrahedralization of the final mesh with full sizedview. (c) Tetrahedralization of the final mesh with clipped view. (d) Zoomed-in clipped view of the finalmesh. The convergence histories of (e) E sys and (f) γ avg . Histogram of (g) the dihedral angle distributionand (h) the radius ratio distribution.Table 8: Mesh quality of the Gear02 case θ min /θ max γ min /γ avg θ min θ < θ < θ < θ < N tet runtime N iter Fu et al.
Improved
5. Conclusions
In this paper, a feature-aware SPH formulation is proposed to achieveconcurrent and automated isotropic unstructured mesh generation. Variousnumerical experiments are investigated considering both triangular surface-mesh and tetrahedral volumetric-mesh generation of geometries consistingvarious sharp features (creases, sharp edges and singularity points). The49obustness of the particle relaxation procedure is improved and a remarkablespeedup for obtaining an optimized mesh vertex distribution is accomplished.The main contribution of the paper can be summarized as: Consistent cell-based feature definition system and target-information for-mulations are developed from previous work in [10][28]. The new fea-ture definition allows for more accurate characterization of complexgeometries and particle systems; A feature-aware correction term is introduced to the governing equations toresolve the issue of incomplete kernel support for particles close to theboundary vicinity. Intensive numerical validations demonstrate thatthe robustness of the simulation is improved. The constraint of nulli-fying particle momentum every timestep can be relaxed. Significantlyimproved convergence benefiting from the correction term is observed; A new measurement of the convergence error, i.e. E sys , is proposed basedon the concept of particle specific volume. E sys is evaluated based onparticle summation and no explicit geometric operations are needed.Numerical results show that the proposed criterion captures the con-vergence of the system successfully; A two-phase mesh generation model is proposed to merge the initial meshgeneration step and mesh-quality optimization step into the same setof governing equations characterized by different simulation parame-ters. Numerical experiments and results comparison reveal that theconvergence of feature-size adaptation target is significantly improved50nd the mesh-quality-optimization step achieves the same performanceas the original algorithm; The proposed algorithm is compared with other state-of-the-art variation-based mesh generators. Comparable mesh-quality is obtained by theproposed method with significantly smaller computational costs.In terms of future development, we will further improve the performanceof the proposed method by harnessing the computational power of GPU-based architectures. Moreover, extending the proposed algorithm to anisotropicunstructured-mesh generation is of high interest.
Acknowledgments
Zhe Ji is partially supported by China Scholarship Council (No. 201506290038).Xiangyu Hu acknowledges funding Deutsche Forschungsgemeinschaft (HU1527/10-1 and HU1527/12-1).
References [1] Q. Du, D. Wang, Recent progress in robust and quality Delaunay mesh generation, Journal ofComputational and Applied Mathematics 195 (1-2) (2006) 8–23.[2] J. Sch¨oberl, Netgen an advancing front 2D/3D-mesh generator based on abstract rules, Computingand visualization in science 1 (1) (1997) 41–52.[3] R. L¨ohner, Recent advances in parallel advancing front grid generation, Archives of ComputationalMethods in Engineering 21 (2) (2014) 127–140.[4] M. S. Shephard, M. K. Georges, Automatic three-dimensional mesh generation by the finite octreetechnique, International Journal for Numerical methods in engineering 32 (4) (1991) 709–749.[5] M. A. Yerry, M. S. Shephard, Automatic three-dimensional mesh generation by the modified-octreetechnique, International Journal for Numerical Methods in Engineering 20 (11) (1984) 1965–1990.[6] J. R. Shewchuk, Delaunay refinement algorithms for triangular mesh generation, Computationalgeometry 22 (1-3) (2002) 21–74.
7] L. P. Chew, Guaranteed-quality Delaunay meshing in 3D (short version), in: Proceedings of thethirteenth annual symposium on Computational geometry, ACM, 1997, pp. 391–393.[8] S. Ni, Z. Zhong, Y. Liu, W. Wang, Z. Chen, X. Guo, Sliver-suppressing tetrahedral mesh optimizationwith gradient-based shape matching energy, Computer Aided Geometric Design 52 (2017) 247–261.[9] Q. Du, V. Faber, M. Gunzburger, Centroidal voronoi tessellations: Applications and algorithms,SIAM review 41 (4) (1999) 637–676.[10] L. Fu, L. Han, X. Y. Hu, N. A. Adams, An isotropic unstructured mesh generation method basedon a fluid relaxation analogy, Computer Methods in Applied Mechanics and Engineering 350 (2019)396–431.[11] Z. Zhong, X. Guo, W. Wang, B. L´evy, F. Sun, Y. Liu, W. Mao, Particle-based anisotropic surfacemeshing, ACM Transactions on Graphics (TOG) 32 (4) (2013) 99.[12] J. R. Bronson, J. A. Levine, R. T. Whitaker, Particle systems for adaptive, isotropic meshing of CADmodels, in: Proceedings of the 19th International Meshing Roundtable, Springer, 2010, pp. 279–296.[13] N. Chrisochoides, Parallel mesh generation, in: Numerical solution of partial differential equationson parallel computers, Springer, 2006, pp. 237–264.[14] S. J. Owen, A survey of unstructured mesh generation technology., in: IMR, 1998, pp. 239–267.[15] P. J. Frey, P.-L. George, Mesh generation: application to finite elements, ISTE, 2007.[16] C. Jamin, P. Alliez, M. Yvinec, J.-D. Boissonnat, CGALmesh: a generic framework for Delaunaymesh generation, ACM Transactions on Mathematical Software (TOMS) 41 (4) (2015) 23.[17] L. Chen, J.-c. Xu, Optimal Delaunay triangulations, Journal of Computational Mathematics (2004)299–308.[18] J. Tournois, R. Srinivasan, P. Alliez, Perturbing slivers in 3D Delaunay meshes, in: Proceedings ofthe 18th international meshing roundtable, Springer, 2009, pp. 157–173.[19] S. Soner, C. Ozturan, Generating multibillion element unstructured meshes on distributed memoryparallel machines, Scientific Programming 2015 (2015) 8.[20] D. Feng, A. N. Chernikov, N. P. Chrisochoides, Two-level locality-aware parallel Delaunay image-to-mesh conversion, Parallel Computing 59 (2016) 60–70.[21] B. M. Klingner, J. R. Shewchuk, Aggressive tetrahedral mesh improvement, in: Proceedings of the16th international meshing roundtable, Springer, 2008, pp. 3–23.[22] J. Chen, D. Zhao, Y. Zheng, Y. Xu, C. Li, J. Zheng, Domain decomposition approach for parallelimprovement of tetrahedral meshes, Journal of Parallel and Distributed Computing 107 (2017) 101–113.[23] Q. Du, D. Wang, Tetrahedral mesh generation and optimization based on Centroidal Voronoi Tes-sellations, International journal for numerical methods in engineering 56 (9) (2003) 1355–1373.[24] M. Meyer, R. Whitaker, R. M. Kirby, C. Ledergerber, H. Pfister, Particle-based sampling and meshingof surfaces in multimaterial volumes, IEEE Transactions on Visualization and Computer Graphics [33] J. Nickolls, I. Buck, M. Garland, K. Skadron, Scalable Parallel Programming with CUDA, Queue6 (2) (2008) 40–53. doi:10.1145/1365490.1365500 .URL http://doi.acm.org/10.1145/1365490.1365500 [34] A. J. Crespo, J. M. Dom´ınguez, B. D. Rogers, M. G´omez-Gesteira, S. Longshaw, R. Canelas, R. Va-condio, A. Barreiro, O. Garc´ıa-Feal, DualSPHysics: Open-source parallel CFD solver based onSmoothed Particle Hydrodynamics (SPH), Computer Physics Communications 187 (2015) 204–216.[35] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of computationalphysics 117 (1) (1995) 1–19.[36] P. Incardona, A. Leo, Y. Zaluzhnyi, R. Ramaswamy, I. F. Sbalzarini, OpenFPM: A scalable openframework for particle and particle-mesh codes on parallel computers, Computer Physics Communi-cations.[37] L. Fu, X. Y. Hu, N. A. Adams, Adaptive anisotropic unstructured mesh generation methodbased on fluid relaxation analogy, Communications in Computational Physics, doi:10.4208/cicp.OA-2019-0049 .[38] J. J. Monaghan, Smoothed particle hydrodynamics, Annual review of astronomy and astrophysics30 (1) (1992) 543–574.[39] G. Contreras, M. Martonosi, Characterizing and improving the performance of intel threading build-ing blocks, in: Workload Characterization, 2008. IISWC 2008. IEEE International Symposium on,IEEE, 2008, pp. 57–66.
40] L. Fu, S. Litvinov, X. Y. Hu, N. A. Adams, A novel partitioning method for block-structured adaptivemeshes, Journal of Computational Physics 341 (2017) 447–473.[41] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based onHamilton-Jacobi formulations, Journal of computational physics 79 (1) (1988) 12–49.[42] L. Chiron, M. De Leffe, G. Oger, D. Le Touz´e, Fast and accurate sph modelling of 3d complex wallboundaries in viscous and non viscous flows, Computer Physics Communications 234 (2019) 93–111.[43] A. Amicarelli, G. Agate, R. Guandalini, A 3d fully lagrangian smoothed particle hydrodynamicsmodel with both volume and surface discrete elements, International Journal for Numerical Methodsin Engineering 95 (5) (2013) 419–450.[44] J. Feldman, J. Bonet, Dynamic refinement and boundary contact forces in SPH with applicationsin fluid flow problems, International Journal for Numerical Methods in Engineering 72 (3) (2007)295–324.[45] C. Hermange, G. Oger, Y. Le Chenadec, D. Le Touz´e, A 3d sph–fe coupling for fsi problems andits application to tire hydroplaning simulations on rough ground, Computer Methods in AppliedMechanics and Engineering 355 (2019) 558–590.[46] D. Shepard, A two-dimensional interpolation function for irregularly-spaced data, in: Proceedings ofthe 1968 23rd ACM national conference, 1968, pp. 517–524.[47] X. Y. Hu, N. A. Adams, A multi-phase sph method for macroscopic and mesoscopic flows, Journalof Computational Physics 213 (2) (2006) 844–861.[48] M. Matsumoto, T. Nishimura, Mersenne twister: a 623-dimensionally equidistributed uniformpseudo-random number generator, ACM Transactions on Modeling and Computer Simulation(TOMACS) 8 (1) (1998) 3–30.[49] L. Fu, Z. Ji, X. Y. Hu, N. A. Adams, Parallel fast-neighbor-searching and communication strategyfor particlebased methods, Engineering Computations 36 (3) (2019) 899–929.[50] Z. Ji, L. Fu, X. Y. Hu, N. A. Adams, A Lagrangian Inertial Centroidal Voronoi Particle methodfor dynamic load balancing in particle-based simulations, Computer Physics Communications 239(2019) 53–63.[51] H. Si, Tetgen, a Delaunay-based quality tetrahedral mesh generator, ACM Transactions on Mathe-matical Software (TOMS) 41 (2) (2015) 11.[52] B. R. Hoehn, P. Oster, T. Tobie, K. Michaelis, Test methods for gear lubricants, Goriva i maziva47 (2) (2008) 141–152.40] L. Fu, S. Litvinov, X. Y. Hu, N. A. Adams, A novel partitioning method for block-structured adaptivemeshes, Journal of Computational Physics 341 (2017) 447–473.[41] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based onHamilton-Jacobi formulations, Journal of computational physics 79 (1) (1988) 12–49.[42] L. Chiron, M. De Leffe, G. Oger, D. Le Touz´e, Fast and accurate sph modelling of 3d complex wallboundaries in viscous and non viscous flows, Computer Physics Communications 234 (2019) 93–111.[43] A. Amicarelli, G. Agate, R. Guandalini, A 3d fully lagrangian smoothed particle hydrodynamicsmodel with both volume and surface discrete elements, International Journal for Numerical Methodsin Engineering 95 (5) (2013) 419–450.[44] J. Feldman, J. Bonet, Dynamic refinement and boundary contact forces in SPH with applicationsin fluid flow problems, International Journal for Numerical Methods in Engineering 72 (3) (2007)295–324.[45] C. Hermange, G. Oger, Y. Le Chenadec, D. Le Touz´e, A 3d sph–fe coupling for fsi problems andits application to tire hydroplaning simulations on rough ground, Computer Methods in AppliedMechanics and Engineering 355 (2019) 558–590.[46] D. Shepard, A two-dimensional interpolation function for irregularly-spaced data, in: Proceedings ofthe 1968 23rd ACM national conference, 1968, pp. 517–524.[47] X. Y. Hu, N. A. Adams, A multi-phase sph method for macroscopic and mesoscopic flows, Journalof Computational Physics 213 (2) (2006) 844–861.[48] M. Matsumoto, T. Nishimura, Mersenne twister: a 623-dimensionally equidistributed uniformpseudo-random number generator, ACM Transactions on Modeling and Computer Simulation(TOMACS) 8 (1) (1998) 3–30.[49] L. Fu, Z. Ji, X. Y. Hu, N. A. Adams, Parallel fast-neighbor-searching and communication strategyfor particlebased methods, Engineering Computations 36 (3) (2019) 899–929.[50] Z. Ji, L. Fu, X. Y. Hu, N. A. Adams, A Lagrangian Inertial Centroidal Voronoi Particle methodfor dynamic load balancing in particle-based simulations, Computer Physics Communications 239(2019) 53–63.[51] H. Si, Tetgen, a Delaunay-based quality tetrahedral mesh generator, ACM Transactions on Mathe-matical Software (TOMS) 41 (2) (2015) 11.[52] B. R. Hoehn, P. Oster, T. Tobie, K. Michaelis, Test methods for gear lubricants, Goriva i maziva47 (2) (2008) 141–152.