Prismatic mesh generation based on anisotropic volume harmonic field

Abstract

In this paper, we present an effective prismatic mesh generation method for viscous flow simulations. To address the prismatic mesh collisions in recessed cavities or slit areas, we exploit 3D tensors controlled anisotropic volume harmonic field to generate prismatic meshes. Specially, a well-fitting tetrahedral mesh is first constructed to serve as the discrete computation domain of volume harmonic fields. Then, 3D tensors are exploited to control the volume harmonic field that better fits the shape geometry. From the topological perspective, the generation of the prismatic mesh can be treated as a topology-preserved morphing of the viscous wall. Therefore, iso-surfaces in the volume harmonic field should be homeomorphic to the viscous wall while fitting its shapes. Finally, a full prismatic mesh can be induced by estimating the forward directions and visible regions in the volume harmonic field. Moreover, to be compatible with different simulation situations, the thickness of prismatic meshes should be variable. Our approach provides local adjustable thickness of prismatic meshes, which can be achieved by controlling local 3D tensors. The proposed anisotropic volume harmonic field based prismatic meshes are efficient and robust, and a full prismatic mesh can be guaranteed without low quality collisions. Various experiments have shown that our proposed prismatic meshes have obvious advantages in terms of efficiency and effectiveness.

1 Introduction

Mesh generation is fundamental to physical simulation and plays a critical role in modern industry, especially in the aerospace field. The shape design of aircraft’s major components, such as wings, emages, engine inlet, etc., requires numerical simulation tools to determine the aerodynamic force and torque for them [1]. As a discretized computational domain in Computational Fluid Dynamics (CFD), the quality of mesh directly affects the calculation accuracy and convergence rate of the numerical simulation. For the complex geometries, it is still a significant challenge to generate high-quality Reynolds-Averaged Navier–Stokes (RANS) mesh due to the existence of recessed cavities or slit areas [2].

The prismatic mesh generation can be treated as surface offset or morphing of original meshes, how to control the directions and movement distances of front nodes are the key parts. During the aerodynamic simulation of high Reynolds number flow, it is essential to adopt layered anisotropic prisms perpendicular to the object to capture the boundary layer near the viscous wall. Most researches are mainly devoted to advancing front methods [3–13] and a few approaches focus on solving the partial differential equation (PDE) [12–23]. Recently, the prismatic hybrid mesh is a prevailing grid type for constructing the external flow field because of its good performance in the balance of precision and efficiency [3, 4, 14, 24–27]. The central idea of front methods is to march outward along the weighted normal of each front node until it terminates when the current front node collides with another front node. It should be noted that a large amount of self-intersection detection and optimization are required to preserve mesh quality, the advancing front methods are difficult to meet the application requirements to some extent, especially for complex models. The PDE based methods generally solve the PDEs, such as Laplace equation [22] or Eikonal equation [17, 28], as the governing equation, and assign the gradient of solutions to the marching directions. Generally, the PDE-based methods are time-consuming [16–18, 22, 28], but the directions of front nodes are smoother and the possible global collision positions could be directly figured out [14, 16]. In this approach, our method is inspired by the above PDE-based methods, and we are working to promote efficiency and mesh quality.

From the perspective of topology, prismatic mesh generation is equivalent to a special topology-preserved morphing of the viscous wall. A full prismatic mesh can be obtained if there is a family of offset surfaces that are homeomorphic to the viscous wall. Similarly, if the iso-surfaces of volume harmonic filed satisfy the above conditions, we can generate a full prismatic mesh. Therefore, how to find a rational volume harmonic field [29, 30] in the background mesh is the key. The layer thicknesses of prismatic meshes are related to the adjacent iso-surfaces of the corresponding volume harmonic field. To make the prismatic mesh thickness controllable, inspired by the research of [31], we introduce 3D tensors as constraints to locally adjust the volume harmonic field anisotropically without bringing in the saddle points. The anisotropic harmonic field presented in this approach is a general method that can be extended to different situations. For example, for two objects close to each other, anisotropic tensor constraints can be applied to delay the occurrence of self-intersection. Moreover, different from the traditional PDE-based methods, we do not directly adopt the solution gradient as the marching direction, but to be the initial direction. The iso-surfaces of volume harmonic field and visible regions [10, 11] are considered as the constraints of marching distance and direction for front nodes severally. These strategies ensure the effectiveness of prismatic mesh generation and establish the connection between the mesh quality and the volume harmonic field.

In addition, the background mesh used to calculate volume harmonic field also has an effect on prismatic mesh quality. If the topology of the envelope surface and the viscous wall are inconsistent, the volume harmonic field inevitably has saddle points, where the topological transformation will occur between different iso-surfaces. To avoid the problem of saddle points, we exploit a Minkowski Sum boundary surface [32] as the envelope surface. There is always a Minkowski Sum boundary surface that has the same topology as the viscous wall, and the boundary surface and the viscous wall will be infinitely close in extreme cases. Compared to the ball or box, the background mesh constructed by Minkowski Sum not only preserves the topology, but also eliminates vast redundant cells, which greatly improves the efficiency of calculation and reduces memory consumption. In practice, the construction of unstructured tetrahedral meshes is accomplished by TETGEN software [33], which provides a strong foundation for the practicability of our algorithm with its comprehensive functions and robustness.

The remaining part of the paper proceeds as follows. Section 2 presents an overview of the proposed pipeline. Section 3 briefly introduces the theory-background of the anisotropic volume harmonic field and the corresponding discretization formula. Section 4 expands the implementation details, such as the construction of the Minkowski Sum boundary surface, the computation of the tensor constraints, and the generation of anisotropic volume harmonic filed and prismatic mesh. The prismatic mesh generation process of our samples with varying degrees of complexity is demonstrated in Section 5. Section 6 is dedicated to the outcomes of the study.

2 Overview of proposed framework

2.1 Algorithm inputs

Our algorithm supports discrete tessellation surfaces as data inputs, including triangular meshes and quadrilateral meshes, and the final outputs correspond to tri-prism meshes and hexahedral meshes, respectively. In our framework, the input surface is taken as the viscous wall.

The parameter input contains :

• 1 Layers number =n, the required number of prismatic mesh layers.

• 2 First layer thickness =h₁, the thickness of first layer mesh closest to the viscous wall.

• 3 Stretch =α, the stretch factor between the thickness of the k layer and the k+1 layer, k=1,2,...,n−1.

Given the above parameters, the thickness of each layer can be calculated. We select appropriate vertices on the viscous wall and trace them to the envelope surface along the gradient of the volume harmonic field. Then, according to the thickness h_i, we calculate the sampling energy e_i along the trajectories, which will be used to compute the iso-surfaces. This strategy establishes a relation between the thicknesses of mesh layers and iso-surfaces of volume harmonic field. Due to the characteristics of the volume harmonic field, the generated prismatic mesh remains smooth, avoiding the inevitable collision when the front nodes advance equidistantly at the concave. Moreover, the local thickness can be adjusted directly by locally controlling the shape of the iso-surfaces.

2.2 Algorithm framework

The framework of the entire algorithm is presented in Fig. 1.

Fig. 1

Prismatic mesh generation framework

The main steps of prismatic mesh generation based on the anisotropic volume harmonic field are outlined as follows:

• 1 Calculating the expected thickness h_i,i=1,⋯,n according to the input parameters.

• 2 Computing the spacing between the envelope mesh and the viscous wall, formulated as $h_{e}=t \ast (\sum _{i=0}^{n-1} \alpha ^{i} \ast h_{1})$, which is t times the overall height of the prismatic mesh. The scaling factor t is a number not less than 1 (t=1.1 by default).

• 3 Setting h_e to be the radius of a ball in the calculation of Minkowski Sum boundary surface, which is the envelope surface. The envelope surface may have inferior quality cells, which need to be improved by Laplace smooth.

• 4 Generating a background tetrahedral mesh between the envelope surface and the viscous wall by TETGEN software with a tessellation coefficient d (d=1.1 by default for a balance between accuracy and efficiency).

• 5 Constructing tensor constraints according to the input parameters. There is a 3D tensor for each point, and it can be used to control the harmonic field locally.

• 6 Calculating an anisotropic volume harmonic field by optimizing the Laplace equation using iteration method. Here, the energy of vertices on envelope surface is set to 0 and the energy of viscous wall vertices is set to 1.

• 7 Calculating the sampling energy e_i according to the expected thickness h_i.

• 8 Computing the marching directions and distances of front nodes, and constructing the final prismatic mesh.

3 Anisotropic volume harmonic field

Given a viscous wall, there is an appropriate boundary constraint so that a volume harmonic field without saddle points can be constructed. Then, the iso-surfaces can induce a generation of full prismatic mesh. Generally, it is easy to construct a volume harmonic field, and induce a prismatic mesh with smooth layer thickness. However, the local control of layer thickness is necessary to better fit the shape and applications, therefore, a local controllable anisotropic volume harmonic field is needed. In this section, the theory involved in the anisotropic volume harmonic field is explained, and the implementation is detailed in Section 4.1.

3.1 Volume harmonic field

In this subsection we will briefly introduce the traditional volume harmonic field, which can be treated as a steady-state heat conduction in a three-dimensional manifold. In mathematics, a volume harmonic field H:V↦R is a scalar field defined on the volumetric mesh M, where V represents the vertex set of M. H is a solution of Laplace equation

$$\begin{array}{@{}rcl@{}} &\Delta H = div(\nabla H) = 0, \end{array} $$

(1)

subject to Dirichlet boundary conditions

$$\begin{array}{@{}rcl@{}} &H(v_{i}) = c_{i}, &v_{i} \in V_{B} \subseteq V, \end{array} $$

(2)

where c_i is a constant and V_B is denoted as the set of sources and sinks. The harmonic field construction can be regarded as a heat distribution when the heat diffuses from sources to sinks until it reaches a steady-state. To solve the Laplace equation in a discretized computational domain, the piece-wise linear Laplace operator is expressed as

$$\begin{array}{@{}rcl@{}} &\Delta H(v_{i}) = \sum\limits_{v_{j} \in N(v_{i})} W(e_{{ij}}) (H(v_{i}) - H(v_{j})), \end{array} $$

(3)

where N(v_i) is the neighborhood of v_i,W(e_ij) is a scalar weight assigned to the edge e_ij. Then, the discrete formulation of Laplace Eq. 1 can be written into a sparse linear system

$$\begin{array}{@{}rcl@{}} &\mathbf{L}\mathbf{H} = 0, \end{array} $$

(4)

where L is a matrix that expresses the discrete Laplace-Beltrami operator

$$\begin{array}{@{}rcl@{}} &L_{{ij}}=\left\{ \begin{array}{rcl} \sum\limits_{v_{k} \in N(v_{i})}W(e_{{ik}}) & & {i = j,} \\ -W(e_{{ij}}) & & {v_{j} \in N(v_{i}),} \\ 0 & & {otherwise.} \end{array} \right. \end{array} $$

(5)

The weighting coefficient W(e_ij) commonly uses the standard cotangent scheme [29]. It is noteworthy that different weight schemes and Dirichlet conditions will lead to different harmonic fields. In the following section, we focus on this perspective and expand the description of the anisotropic volume harmonic field construction.

3.2 Anisotropic volume harmonic field with tensor constraints

The harmonic field can be treated as a heat diffusion, which can be obtained by solving an elliptic PDE with the constraint of the heat sources and boundaries. From the viewpoint of geometry, the harmonic field has some good properties, such as no-curl, no-divergence and smoothness, and they can be exploited to tackle lots of graphic issues. Most of the PDE-based methods only take into account the isotropy. For clarity, the isotropy mentioned here is used for the construction of traditional PDE equation without anisotropic constraints. However, the local layer thickness of prismatic mesh induced by the isotropic harmonic field lacks control. Normally, it is possible to change the boundary constraints to control the layer thickness globally, but still difficult to local layer thickness. To cope with the control of local layer thickness and more actual demands, we propose a 3D tensor guided anisotropic volume harmonic field.

A 3D tensor is defined on each vertex in a three-dimensional manifold, and the tensor can be represented by a normalized orthogonal frame [x₁,x₂,x₃] and the corresponding scalar factors [γ₁,γ₂,γ₃]

$$\begin{array}{@{}rcl@{}} &\mathbf{T}(v_{i}) = \gamma_{1} \mathbf{x_{1}} \mathbf{x_{1}}^{T} + \gamma_{2} \mathbf{x_{2}} \mathbf{x_{2}}^{T} + \gamma_{3} \mathbf{x_{3}} \mathbf{x_{3}}^{T}. \end{array} $$

(6)

T(v_i) is a 3×3 symmetric matrix, and the Laplace Eq. 4 with tensor constraints can be reformulated as

$$\begin{array}{@{}rcl@{}} &\mathbf{\widetilde{L}}\mathbf{\widetilde{H}} = 0, \end{array} $$

(7)

where

$$\begin{array}{@{}rcl@{}} &\widetilde{L}_{{ij}}=\left\{ \begin{array}{rcl} \sum\limits_{v_{k} \in N(v_{i})} \widetilde{W}(e_{{ik}}) & & {i = j,} \\ - \widetilde{W}(e_{{ij}}) & & {v_{j} \in N(v_{i}),} \\ 0 & & {otherwise,} \end{array} \right. \end{array} $$

(8)

and

$$\begin{array}{@{}rcl@{}} &\widetilde{W}(e_{{ij}}) = exp \left(\ \left(\frac{(v_{i} - v_{j})^{T}}{\| v_{i} - v_{j} \|} (\mathbf{T}(v_{i}) + \mathbf{T}(v_{j})) \frac{(v_{i} - v_{j})}{\| v_{i} - v_{j} \|} \right) \ / \ \delta \ \right). \end{array} $$

Where δ is a penalty factor to characterize the influence of the tensor field on the harmonic field. The anisotropic harmonic field $\mathbf {\widetilde {H}}$ is the solution of the Laplace Eq. 7 under the Dirichlet boundary conditions ($\mathbf {\widetilde {H}}(v_{i}) = c_{i}, v_{i} \in V_{B}$).

4 Prismatic mesh generation framework

4.1 Construction of anisotropic volume harmonic field

This section presents the implementation details of constructing an anisotropic volume harmonic field, including the generation of the envelope surface, the tetrahedral background mesh and the tensor constraints.

4.1.1 Minkowski Sum boundary surfaces

PDE-based methods need to construct background grids as the calculation domain. Generally, the viscous wall corresponds to the inner boundary of the background mesh, and the envelope surface corresponds to the outer boundary. Ellipsoids and cubes are commonly used to be the envelope surfaces in traditional researches. However, there may be large distortions of iso-surfaces if the inner boundary and the outer boundary have big shape differences, especially for the topology changes. The distorted iso-surfaces will induce low quality of layer meshes, resulting in non-full prismatic meshes. Therefore, a proper outer boundary that better fits the shape and topology of the viscous wall is necessary. In this approach, we exploit Minkowski Sum to obtain the outer boundary surface, which can be guaranteed to maintain the topological consistency of the inner boundary surface, enabling us to avoid the negative effects of saddle points in the harmonic field. Moreover, the background mesh generated by the proposed Minkowski Sum boundary surface has fewer redundant tetrahedral cells, which also improves calculation efficiency and reduces memory consumption.

Given two polyhedral grids P and Q, the Minkowski Sum of them is defined as

$$\begin{array}{@{}rcl@{}} & P \oplus Q = \{ p + q | p \in P, q \in Q \}, \end{array} $$

(9)

the Minkowski Sum boundary is represented as ∂(P⊕Q). Since the calculation of a volume harmonic field on the background mesh is not sensitive to the accuracy of the outer boundary surface, we adopt an efficient and robust convolution approach proposed in [32] to extract the Minkowski Sum boundary surface. Also, a Laplace smoothing is utilized to improve the quality of the obtained surface.

Given a viscous wall and a sphere mesh, the Minkowski Sum boundary surface of them is related to the radius of the sphere. In our framework, the sphere radius is formulated as

$$\begin{array}{@{}rcl@{}} & r_{e} = \sum_{i=0}^{n-1} \alpha^{i} \ast h_{1}, \end{array} $$

(10)

where n, h₁ and α are the number of layers, the first layer thickness and the stretch factor, respectively. Note that a preliminary judgment on the topology change between the Minkowski Sum boundary surface and the viscous wall is implemented. If there are topology changes, we reduce the sphere radius locally to preserve the topological consistency. We demonstrate the Minkowski Sum boundary surface of a plane model and a sphere mesh, as shown in Fig. 2.

Fig. 2

Illustration of the Minkowski Sum boundary surface (envelope surface). (a) The original surface mesh (viscous wall). (b) The Minkowski Sum boundary surface. (c) The combination of the Minkowski Sum boundary surface and the original surface

After obtaining the boundary surfaces, we will construct the background meshes for the calculation of volume harmonic fields. Here, TetGen [33] is exploited to generate the tetrahedral background mesh, because it is a highly mature tetrahedral mesh generation software that can produce high quality geometry-aware tetrahedral mesh with local tessellation controllable. To balance the accuracy and efficiency, the density coefficient is set to 1.1 in all our experiments.

4.1.2 Anisotropic volume harmonic field

A proper anisotropic volume harmonic field can be used to generate a prismatic mesh that better fits the geometry and topology of viscous wall in the background mesh. The iso-surfaces of the harmonic field determine the layer thickness and quality of prismatic meshes. Therefore, 3D tensors are exploited to control volume harmonic field locally, and they make the prismatic mesh generation more flexible and controllable to tackle complex requirements. This subsection elaborates on the construction of an anisotropic harmonic field.

Construction of 3D tensors

In our method, the 3D tensors are constructed to control the variation rate of the harmonic field gradient. Essentially, the gradient variation rate in a standard volume harmonic field could characterize the geometric and topological changes between the viscous wall and envelope surface. The 3D tensors are derived from the standard volume harmonic field, and then used to guide an anisotropic volume harmonic field that can avoid collisions and reduce mesh distortion.

A 3D tensor can be constructed for each vertex in the background mesh

$$\begin{array}{@{}rcl@{}} &\mathbf{T}(v_{i}) = \sum_{j=1}^{N(v_{i})} e^{\lambda_{j}} \mathbf{e_{{ij}}}\mathbf{e_{{ij}}}^{T}/{||\mathbf{e_{{ij}}}\mathbf{e_{{ij}}}^{T}||}, \end{array} $$

(11)

where N(v_i) is the neighborhood vertex set of v_i,λ_j=(e_ij^T∇H(v_i)e_ij)^β,e_ij is the edge vector from v_i to v_j,∇H(v_i) is the gradient at vertex v_i in the isotropic volume harmonic field H obtained by directly computing the Laplace equation with common cotangent weight, β is a regulating factor (−1 by default).

Compared with the standard volume harmonic field, the one constructed with the proposed tensors can better control the anisotropic iso-surfaces that result in high-quality and easy-controllable prismatic meshes. Three models with different types of features are tested to illustrate the validity, as shown in Fig. 3.

Fig. 3

Illustration of the anisotropic volume harmonic field. (a) The original models with Minkowski Sum boundary surfaces; (b) The section of iso-surfaces induced by standard volume harmonic field; (c) The section of iso-surfaces induced by the proposed anisotropic volume harmonic field

Discrete calculation of anisotropic volume harmonic field

In our method, the tetrahedral background mesh serves as the discrete computational domain. The projection of 3D tensors on edges is treated as the edge weight to calculate the anisotropic volume harmonic field. Also, to improve the calculation efficiency as well as reduce memory consumption, we exploit the iteration method to obtain an approximate solution of the Laplace Eq. 7 instead of directly solving a huge sparse linear system.

Given boundary conditions, $\widetilde {H}(v_{p}) = 1, v_{p} \in V_{{inner}}$ and $\widetilde {H}(v_{q}) = 0, v_{q} \in V_{{outer}}$, where V_inner,V_outer are vertices on the viscous wall and envelope surface, respectively. Based on previous work [29], we define harmonic energy with tensor constraints as

$$\begin{array}{@{}rcl@{}} &K(\widetilde{H}) = &\sum\limits_{e_{{ij}} \in E} \widetilde{W}(e_{{ij}}) \left(\widetilde{H}(v_{i}) - \widetilde{H}(v_{j})\right)^{2}, \end{array} $$

(12)

where

$$\begin{array}{@{}rcl@{}} &\widetilde{W}(e_{{ij}}) = exp(\ T(e_{{ij}}) \ / \ \delta ), \end{array} $$

(13)

where T(e_ij) is the projection of the corresponding tensor T(v_i) on edge e_ij used to control diffusion velocity along e_ij.

The process of minimizing the energy presented in Formula 12 is to make $\widetilde {H}$ harmonic, which is equivalent to updating $\widetilde {H}(v_{i})$ iteratively to achieve $\Delta \widetilde {H}(v_{i}) = 0$. In this sense, each iteration of updating $\widetilde {H}(v_{i})$ can be treated as the energy at vertex v_i spread to vertices in N(v_i) with the diffusion velocity T. Hence, the formula for $\widetilde {H}(v_{i})$ can be written as

$$\begin{array}{@{}rcl@{}} &\widetilde{H}(v_{i}) = \frac{\sum\limits_{v_{j} \in N(v_{i})} \widetilde{W}(e_{{ij}}) \widetilde{H}(v_{j})}{\sum\limits_{v_{j} \in N(v_{i})} \widetilde{W}(e_{{ij}})}. \end{array} $$

(14)

The detailed implementation of anisotropic harmonic field is described in Algorithm 1.

4.2 Generation of prismatic meshes

Prismatic meshes can be generated by researching the marching direction and marching distance of front nodes in the volume harmonic field.

Marching directions

The marching directions can be obtained by calculating the gradients in the anisotropic volume harmonic field. Most PDE-based methods directly adopt gradients as the marching directions so that the generated meshes have pretty orthogonality. However, due to the discrete gradient computation, the gradients at the steep concave parts are usually inaccurate and even have a large deviation. If the gradients are directly taken as the marching directions, then low-quality cells will be introduced. As a piecewise linear representation, each tetrahedral cell of the background mesh is a linear space. If the scalar values of four cell vertices are not equal to each other in the harmonic field, there is a unique iso-surface whose normal is the gradient of harmonic field. Therefore, we first calculate iso-surfaces whose normal is the gradient of harmonic field. Then, a smooth operation is executed on the iso-surfaces to improve the gradient collisions, especially at the concave parts. Finally, the modified gradients used to be the marching directions can be calculated from the iso-surfaces. Note that we need to restrict the marching direction to the visible cone to prevent from bringing negative volume cells straightly.

Marching distances

In our method, the marching distance depends on the distance between adjacent iso-surfaces. Due to the complex geometries, the marching distances in the same layer might not be consistent, especially for steep concave parts and slit features. Therefore, we adjust the iso-surfaces to better fit the expected marching distances. The expected marching distance of each layer can be formulated as

$$\begin{array}{@{}rcl@{}} &h_{i} = \alpha \ast h_{i-1} \quad \quad i = 2, \cdots, n. \end{array} $$

(15)

We establish a connection between the marching distances and the iso-surfaces. Given a vertex v_k∈V_inner, there is a gradient line from the current vertex to the envelope surface. According to {h_i} on the gradient line, we can extract a set of corresponding sampling energy {e_j},j=1,⋯,n. Then, each sampling energy e_i will guide a corresponding iso-surface. Finally, the combination of gradient lines and iso-surfaces yields a high-quality prismatic mesh. Due to the high local-controllability of anisotropic harmonic field, the prismatic meshes could well fit the geometry at complex parts, such as at steep concave parts and slit features. As shown in Fig. 4, the prismatic meshes induced by anisotropic harmonic field have higher quality and better flexibility than the meshes induced by standard harmonic field.

Fig. 4

Illustration of the prismatic meshes. (a) The prismatic meshes induced by standard volume harmonic field; (b) The prismatic meshes induced by the proposed anisotropic volume harmonic field

5 Numerical experiment and discussion

In this section, we briefly illustrate the experimental results. All the experiments were conducted on a PC with 1.60GHz Intel(R) core(TM) i5-8520U CPU, 16 GB RAM and 64-bit Windows 10 operating system. To describe the effectiveness, we applied our algorithm to various models with complex geometric information. Also, the comparisons of prismatic meshes obtained between the standard harmonic field based method and the proposed method are demonstrated.

5.1 General models with complex geometry

The missile model is a closed genus 0 surface with sharp angles and concave edges. The original surface with its envelope surface is demonstrated in Fig. 5a. A high-quality prismatic mesh requires the initial mesh layer should be close to the viscous wall as homogeneous as possible, especially at steep concave parts. The marching distances of front nodes induced by the iso-surfaces in traditional standard harmonic field are much larger at the concave parts, as shown in Fig. 5b. In our method, we automatically add local 3D tensors to adjust the diffusion velocity at the concave regions, and obtain anisotropic iso-surfaces that could better fit the viscous wall, as shown in Fig. 5c Left. Therefore, we could obtain higher quality prismatic mesh than the one based on standard harmonic field, as shown in Fig. 5c Right. Moreover, the prismatic mesh obtained using the proposed method is a full prismatic mesh that contains 650,010 prisms grids.

Fig. 5

Comparison of prismatic meshes. (a) The original surface and envelope surface; (b) and (c) The iso-surfaces and prismatic meshess generated using the standard harmonic field and the proposed anisotropic harmonic field, respectively

We also test our method on more models, as shown in Figs. 6 and 7. The submarine is a closed genus 0 model with concave parts. The unmanned aerial vehicle is a closed genus 0 model with concave parts and narrow slits. We can see that our method generates full prismatic meshes with higher quality, especially for the complex parts, such as concave parts and narrow slits.

Fig. 6

Prismatic mesh of submarine model. (a) The original surface and envelope surface; (b) and (c) The iso-surfaces and prismatic meshes generated using the standard harmonic field and the proposed anisotropic harmonic field, respectively

Fig. 7

Prismatic mesh of unmanned aerial vehicle model. (a) The original surface and envelope surface; (b) and (c) The iso-surfaces and prismatic meshes generated using the standard harmonic field and the proposed anisotropic harmonic field, respectively

The proposed method can be treated as a glorified normal smoothing approach, but it also works well for models with the concave corner less than 90 degrees. One of the leading factors is that the controllable anisotropic volume harmonic field can provide a larger space for marching directions. Figure 8 demonstrates the applicability of our proposed approach to models with small-angle concave features.

Fig. 8

Prismatic mesh of a model with the concave corner less than 90 degrees. It is a full prismatic mesh without collisions

5.2 Complex models with genus

The proposed framework can also be applied to complex models with complex geometries and genus. As shown in Fig. 9, the plane is a closed genus 2 model with concave parts and circular ring structures. Due to the Minkowski Sum based envelope surface, the topology of the obtained prismatic mesh coincides with that of the original models. Therefore, our proposed framework avoids complex segmentation and merging operations for high-genus models. Moreover, due to the local controllability of anisotropic volume harmonic field, our method can also deal much better with the parts with complex geometries. In all our experiments, our proposed framework generates high-quality full prismatic meshes automatically.

Fig. 9

Prismatic mesh of plane model with complex geometries and genus. (a) The original surface and envelope surface; (b) and (c) The iso-surfaces and prismatic meshes generated using the standard harmonic field and the proposed anisotropic harmonic field, respectively

6 Conclusion

Our work, by automatically adding local tensors to adjust the harmonic field, provides more freedom for the mesh generation based on the traditional PDE methods. The thickness of mesh induced by standard harmonic field can only be regulated globally. But the control and flexibility of local mesh are low, which makes the prisms at concave edges and narrow slits inevitably distort too much. Anisotropic harmonic field, to a certain extent, makes up for the defects of the traditional PDE-method and matches well with the prismatic grid generation task. The standard harmonic field can provide mass information, and we utilize it to automate the addition of local tensor constraints. Although the time-consuming is high under this strategy, it is easy to find parallel solutions to improve efficiency.