Abstract

In this paper, a novel hexahedral refinement method is introduced which finds the optimal tradeoff between the number of inserted elements and inserted singularities according to user-prescribed weighting. The input of our algorithm is the hexahedral mesh with quads tagged for refinement. The quad sets to be inserted by sheet inflation are determined by solving a integer optimization problem to minimize the number of inserted elements and inserted singularities while maintaining mesh consistency. Finally, we design the optimized sheet structures by selecting the optimal local shrink set at the cross section of intersecting quad sets to ensure the quality of mesh refinement. The refinement scheme can be applied iteratively until the mesh density meets the requirements. Experimental results for some mechanical parts verified the effectiveness of the proposed method.

1 Introduction

In finite element analysis, hexahedral meshes are more desirable than tetrahedral meshes, as they offer several advantages, such as providing shape functions with additional terms to increase the accuracy and directional sizing without loss of accuracy [1].

Many automatic all-hexahedral mesh generation methods have been developed over the past decade and some methods offer size control and adaptive meshes generation [2]. Mesh refinement is one of the commonly used techniques of mesh adaptation, which involves increasing the mesh density by inserting elements in a local region to improve the mesh quality and increase the accuracy of the analysis result [3]. However, it is still difficult to refine local regions of hexahedral meshes due to the non-local propagation of topological modification.

Hexahedral mesh refinement methods can be mainly divided into two categories: those based on topological operations [49] and those based on templates [1,1015]. Mitchell and Tautges [5] proposed a refinement operation to increase element quality by removing doublets, also known as pillowing. In order to solve the problem that the density difference between the source surface and the target surface is too large for the meshes generated on non-strict swept volume, Borden et al. [6] proposed the Cleave and Fill algorithm. By finding the shortest path on the surface and extending it to the interior for layer insertion operation, the local refinement is effectively carried out. However, the applicability of this method is limited as it is only suitable for swept volume meshes. Tchon et al. [3] proposed a quadrilateral and hexahedral mesh refinement method to meet the size requirements of an anisotropic field. In this work, pillowing operation is applied to refine the mesh, and sheets are continuously inserted in all directions according to a size specification map defined as a Riemannian metric. Aiming at the problem of continuous mesh deformation during forging, Park and Yang [8] developed a hexahedral mesh refinement technique characterized by an iterative process of inserting zero-thickness element layers and an optimal distribution of mesh density. Shen et al. [9] proposed a physics-based hexahedral mesh adaptation method based on posterior-error estimation using graph cuts algorithm for mesh refinement. Refinement template-based methods are widely used for octree-based hexahedral mesh generation to capture the details of geometrical models. Schneiders et al. [10] proposed four templates for the node, edge, face, and volume refinement. Although the efficiency of this method is high, it is unable to handle concavities. Zhang et al. [1] later expanded the number of templates to six, which not only led to a smoother refinement region, but also solved the problem of the concave areas. Other works including [1115] are basically on similar templates to refine the meshes, so as to improve the mesh quality. The problem of existing refinement methods is that they do not take into account the number of introduced singularities which will influence the topological structure and quality of hexahedral mesh.

In this paper, we present a hexahedral refinement method based on topological operations which allows the user to specify a tradeoff between the number of inserted singularities and number of inserted mesh elements. Compared with existing methods, the advantages of our algorithm are mainly reflected in the following two aspects. First, when determining the quad set, the number of singularities can be substantially controlled to obtain desired sheet. Second, the requirements for the input quads are flexible, such as fewer restrictions on the number and position of input quads. And the control over the number of inserted quads is added to the optimization objective.

2 Algorithm Overview

The input of the algorithm is a coarse hexahedral mesh H = (V, E, F, C) with vertices V, edges E, faces F, and hexahedral elements C, and the quads to be refined which called FO (FO is the subset of F). We define the refinement as a function over all faces R : F → {0, 1} where R(f) = 1 means that face f is going to be duplicated by sheet inflation operation and the edge perpendicular to f will be split. The quad set input for sheet inflation must initiate and terminate at a boundary.

In order to find the minimal-cost solution to refine the hexahedral mesh, the following two issues need to be addressed: (1) how to determine the quad sets to be inflated while balancing the number of inserted elements and inserted singularities?, and (2) how to design sheet generation schemes to guarantee high quality, especially for intersecting quad sets? To solve the first problem, we derive an integer program formulation which considers the tradeoff between the number of inserted elements and singularities according to user demands. To solve the second problem, we optimize the local shrink sets of the quad set based on a quality prediction, and then generate the optimized topological structure of the inserted sheet according to the designed schemes.

As illustrated in Fig. 1, the proposed method mainly consists of two following steps:

  1. Determination of quad sets to be inflated by integer optimization.

  2. Overall mesh refinement by optimized sheet inflation.

Fig. 1
Flowchart of local mesh refinement: (a) a geometric model, (b) initial coarse hexahedral mesh, (c) input quads selected for refinement, (d) quad sets determined using integer optimization, (e) sheets generated after sheet inflation, and (f) result mesh after sheet inflation
Fig. 1
Flowchart of local mesh refinement: (a) a geometric model, (b) initial coarse hexahedral mesh, (c) input quads selected for refinement, (d) quad sets determined using integer optimization, (e) sheets generated after sheet inflation, and (f) result mesh after sheet inflation
Close modal

3 Determination of Quad Sets to be Inflated by Integer Optimization

A conforming refinement can always be found by splitting every element to be refined, propagating along the splitting direction until reaching the boundary. In this way, a large part of mesh will be refined leading to a large increase in the number of hexahedral elements. To avoid this situation, we place the singularities at appropriate positions to make the refinement of local by determining the proper input quad sets of the sheet inflation.

To find a tradeoff between the number of inserted elements and the singularities when inserting sheets, we develop an integer optimization to obtain the optimal quad sets for sheet inflation.

When we use sheet inflation operation to achieve mesh refinement, a key problem is the selection of reasonable quad sets which satisfy both quality and size requirements of users. To insert a sheet at an appropriate position with the minimum number of singularities, we formulate this task as an integer optimization problem.

For every inner quad fi of mesh, a binary variable si is used to indicate whether the face is to be inflated or not. So the refinement function can also be defined as R(fi) = si. One of our objective for integer optimization is Es to minimize the number of selected quads
Es=fiFinnersi
(1)
where Finner refers to the set containing all inner quads of mesh.
Suppose the user defines a set P of quads as the input which must be inflated in order to refine the specified regions. The input set P may come from many application scenarios, such as mesh adaptation based on the analysis results to make the analysis more accurate, mesh refinement by selecting quads in areas with sparse density, mesh matching by selecting the quads which pass through the missing boundary chord, and topology improvement by selecting quads that pass through a specific type of topological structure. Then, si for each quad in P is set to 1
si=1fiP
(2)
In addition, to obtain the valid input quad set of sheet inflation, the quad set should initiate at the boundary and terminate at the boundary. That is to say, the number of selected quads adjacent to every inner mesh edge should be even. Using an auxiliary variable a for every inner edge eEinner, the constraint can be expressed as follows:
si=2aa=0,1,2eEinner
(3)
where the sum of si corresponds to the number of selected quads adjacent inner edge e. Here, we limit a to 0, 1, 2. When a is equal to 1, the total number of selected quads adjacent to the inner edge is 2, which means that there will be one sheet inserted as shown in Fig. 2(a). When a is equal to 2, the total number of selected quads adjacent to the inner edge is 4, which means that there will be a self-intersecting sheet inserted or two crossing sheets inserted as shown in Fig. 2(b).
Fig. 2
Different situations of inner edge in the quad sets: (a) bold inner edge of one simple quad set is adjacent to two quads and (b) every bold inner edge at the intersection of two quad sets is adjacent to four quads
Fig. 2
Different situations of inner edge in the quad sets: (a) bold inner edge of one simple quad set is adjacent to two quads and (b) every bold inner edge at the intersection of two quad sets is adjacent to four quads
Close modal
Since the topological quality of the quad set after sheet inflation is influenced by the inserted singularities, we want to detect and penalize the cost of introducing singularities. When two adjacent quads are selected to be inflated as shown in Fig. 3(a), then singularities will be introduced as shown in Fig. 3(b). To detect and penalize this situation, a binary variable pi is used to denote every possible situation leading to singularities. Therefore, when singularities occur, the pi is set to 1. The constraint can be expressed as follows:
2pis0+s12pi+1
(4)
where s0 and s1 correspond to two adjacent quads which belong to Finner.
Fig. 3
The case of introducing singularities: (a) two adjacent quads in the same element are chosen for inflation which is penalized in the formulation and (b) mesh after inflating two quads selected in (a) and introduction of singularities
Fig. 3
The case of introducing singularities: (a) two adjacent quads in the same element are chosen for inflation which is penalized in the formulation and (b) mesh after inflating two quads selected in (a) and introduction of singularities
Close modal
To balance two different objectives of minimizing the number of inserted mesh elements and that of inserted singularities, a parameter ci is specified by the user and the another objective of our integer optimization is set as follows:
Ep=pici
(5)
Based on the aforementioned approach, the constraints are combined together and the integer optimization program is formulated as
minimizeE=Es+Epsubjectto(2)(3)(4)
(6)
The solution provides the valid set of quads for sheet inflation. We solve this problem using Gurobi solver [16].

4 Overall Mesh Refinement by Optimized Sheet Inflation

The quads determined using the integer optimization program will be inflated into hexahedral elements by sheet inflation operation. To achieve effective mesh refinement, we divide the connected quads into groups according to their connectivity relationship. Then the shrink sets corresponding to quad sets are optimized determined by predicting the mesh quality after sheet inflation operation. Specifically, the intersecting quad sets are handled according to pre-designed schemes. The specific steps are as follows:

4.1 Determination of Connected Quad Sets and Shrink Sets.

During a sheet inflation operation, the nodes on the continuous quad surface are duplicated and reconnected to form a layer structure. Because the quads obtained in the above section cannot be inflated as a whole, we group the quads into different sets according to their mesh connectivity. Specifically, a quad set is formed by interactively collecting quads connected to each other.

Generally speaking, there are two directions in the two sides of a quad set that can be chosen for sheet inflation. The hexahedral set adjacent to one side of the quad set Q is called shrink set. The shrink set can guide the direction of node splitting, so it influences the topological structure of the sheet. To obtain high-quality sheets, we optimize the shrink sets by analyzing the geometry around the quad sets. In this paper, to determine the splitting direction of each node, we introduce the concept of local shrink set [17] as follows:

Definition 4.1

Given a hexahedral meshM = (H, F, E, V) and an input quad setQ = (Fq, Eq, Vq), for every nodevinVq, Qdivides the hexahedron setHvHintonparts (n = 1, 2, 3, 4). Each part represents a potential splitting direction of nodev. The local shrink setSvof nodevcould include several parts ofHv. Figure 4 shows two cases of local shrink set.

Fig. 4
Two cases of local shrink set: (a) a quad set, (b) quad set divides hexahedral elements into two parts and H1 or H2 could be the local shrink set of splitting node, and (c) quad set divides hexahedral elements into four parts and local shrink set could include H1, H2, and H3
Fig. 4
Two cases of local shrink set: (a) a quad set, (b) quad set divides hexahedral elements into two parts and H1 or H2 could be the local shrink set of splitting node, and (c) quad set divides hexahedral elements into four parts and local shrink set could include H1, H2, and H3
Close modal

The specific rules to determine the local shrink set are as follows:

  1. When the quad set Q has only one side, then the hexahedral set adjacent to Q is set as the local shrink set.

  2. When the mesh nodes on Q are on the geometric edge on one side of Q, then the hexahedral set on the other side of Q is selected as the local shrink set as shown in Fig. 5(a). Otherwise, the operation will lead to the poor-quality elements.

  3. When the mesh nodes on Q are on the geometric edge, then the local shrink sets of all the nodes both on Q and on the same geometric edge should be located on the same side of Q as shown in Fig. 5(b). Otherwise, the sheet inflation operation will generate an element across the geometric edge.

  4. When the mesh nodes on Q are on the concave geometric edge, then the hexahedral set facing the inside of the model is selected as the shrink set as shown in Fig. 5(c). Otherwise, the sheet inflation operation will fail.

  5. For other cases, the hexahedral set on one side of Q is randomly selected as shrink set.

Fig. 5
Three situations to determine the shrink set of the mesh nodes on the quad set (the vertical line represents the quad set, the horizontal line represents the geometric edge, and the dot represents the mesh nodes to be split): (a) mesh nodes on Q are on the geometric edge on one side of Q, (b) mesh nodes on Q are on the geometric edge, and (c) mesh nodes on Q are on the concave geometric edge
Fig. 5
Three situations to determine the shrink set of the mesh nodes on the quad set (the vertical line represents the quad set, the horizontal line represents the geometric edge, and the dot represents the mesh nodes to be split): (a) mesh nodes on Q are on the geometric edge on one side of Q, (b) mesh nodes on Q are on the geometric edge, and (c) mesh nodes on Q are on the concave geometric edge
Close modal

4.2 Generation of Sheet Structure.

For the general case of generating the structure of a sheet, we first split the nodes on the quad set to V1 and V2, where V1 is at the original position as shown in Fig. 6(a) and V2 is located in the direction of shrink set as shown in Fig. 6(b). Then, the hexahedral elements included in the shrink set will be deleted. Finally, the sheet is generated by connecting the splitting nodes and the nodes around the original mesh as shown in Fig. 6.

Fig. 6
General sheet inflation process: (a) input quad set of sheet inflation, (b) duplicated mesh nodes, and (c) generated sheet
Fig. 6
General sheet inflation process: (a) input quad set of sheet inflation, (b) duplicated mesh nodes, and (c) generated sheet
Close modal

To generate intersecting quad sets where the degree of inner edges is closer to four, we first pre-analyze the mesh topology around the quad set to determine a reasonable generation scheme and improve the quality of the sheet inflation operation. Then we generate the corresponding layer structure according to the template. The specific steps are as follows:

  1. Determination of local shrink set based on quality prediction. Let us assume that there is an intersection line on the self-intersecting quad set or at the intersection of two quad sets as shown in Figs. 7(b) and 8(b). The vertices on the intersection line split the local hexahedral elements into four parts. Based on the experimental experience, we devise two templates including A and B to generate the mesh topology at the intersection as shown in Fig. 7. The local shrink sets shown in Figs. 7(c) and 7(e) correspond to the topological structures as shown in Figs. 7(d) and 7(f), respectively. We compare the topological quality with different generation schemes according to following equations:

    T1=e[(ne12)2+(ne2+nei3)2+(ne42)2]T2=e[(ne11)2+(ne21)2+(ne31)2+(ne41)2]
    where e refers to any mesh edge at the intersection line of quad sets, T1 refers to the quality score after using the template A as shown in Figs. 7(c)7(d), and T2 refer to the quality score after using the template B as shown in Figs. 7(e)7(f). The surrounding hexahedral elements are divided into four parts whose element number is respectively. The smaller T is, the higher the quality after sheet inflation is. Therefore, we choose the better template option with smaller T. In addition, two different types of layer structures generated based on a self-intersecting quad set are illustrated in Fig. 8.
  2. Construction of intersecting sheets. Based on different shrink sets, we construct sheet structures of different templates as follows:

    • For the first scheme, the local shrink set at the intersection structure is shown as the red face set in Fig. 9(a), that is, the red vertex is split into three yellow vertices of Fig. 9(b). Then, the adjacent cells of the face set are deleted, as shown in Fig. 9(c). Finally, the split points and the surrounding points are connected according to the original grid connection relationship. Finally, the layer structure is generated, as shown in Fig. 9(d).

    • For the second scheme, the local shrink set at the intersection structure is shown as the red quad set in Fig. 10(a), where the red point is split into four yellow vertices of Fig. 10(b), Then, the adjacent cells of the face set are deleted, as shown in Fig. 10(c). Finally, the split points and the surrounding points are connected according to the original grid connection relationship. Finally, the layer structure is generated, as shown in Fig. 10(d).

Fig. 7
Two types of topological structures generated by inflating two intersecting quad sets: (a) original hexahedral mesh, (b) two intersecting input quad sets, (c) shrink set of the first sheet generation template A, (d) sheet structure after performing sheet inflation based on (c), (e) shrink set of the second sheet generation template B, and (f) sheet structure after performing sheet inflation based on (e)
Fig. 7
Two types of topological structures generated by inflating two intersecting quad sets: (a) original hexahedral mesh, (b) two intersecting input quad sets, (c) shrink set of the first sheet generation template A, (d) sheet structure after performing sheet inflation based on (c), (e) shrink set of the second sheet generation template B, and (f) sheet structure after performing sheet inflation based on (e)
Close modal
Fig. 8
Two types of topological structures by inflating one self-intersecting quad set: (a) original hex mesh, (b) one self-intersecting input quad set, and (c), (d) two sheet structures after performing sheet inflation based on different shrink sets
Fig. 8
Two types of topological structures by inflating one self-intersecting quad set: (a) original hex mesh, (b) one self-intersecting input quad set, and (c), (d) two sheet structures after performing sheet inflation based on different shrink sets
Close modal
Fig. 9
Two-dimensional schematic of the mesh generation of the intersecting quad sets for the first template: (a) quad sets and their local shrink sets, (b) mesh nodes and their split mesh nodes, (c) deletion of the elements related to the mesh nodes, and (d) new structure obtained by reconnecting the split mesh nodes and their surrounding nodes
Fig. 9
Two-dimensional schematic of the mesh generation of the intersecting quad sets for the first template: (a) quad sets and their local shrink sets, (b) mesh nodes and their split mesh nodes, (c) deletion of the elements related to the mesh nodes, and (d) new structure obtained by reconnecting the split mesh nodes and their surrounding nodes
Close modal
Fig. 10
Two-dimensional schematic of the mesh generation of the intersecting quad sets for the second template: (a) quad sets and their local shrink sets, (b) mesh nodes and their split mesh nodes, (c) deletion of the elements related to the mesh nodes, and (d) new structure obtained by reconnecting the split mesh nodes and their surrounding nodes
Fig. 10
Two-dimensional schematic of the mesh generation of the intersecting quad sets for the second template: (a) quad sets and their local shrink sets, (b) mesh nodes and their split mesh nodes, (c) deletion of the elements related to the mesh nodes, and (d) new structure obtained by reconnecting the split mesh nodes and their surrounding nodes
Close modal

4.3 Iterative Refinement.

In most cases, a single refinement round does not lead to satisfactory density results. Therefore, we can choose to either perform dicing algorithm [18], or to calculate new quad sets again for refinement until the mesh density meets our requirements. If the layer inserted in the previous round is located in an area that needs to be refined, the dicing algorithm is preferred from the perspective of efficiency. The dicing algorithm directly duplicates sheets within an existing mesh, and reconnects the surrounding hexahedral elements as shown in Figs. 11(a) and 11(b). The dicing algorithm can be used to refine the mesh on the original hexahedral mesh structure simply and effectively.

Fig. 11
Iterative mesh refinement using dicing algorithm: (a), (b) the selected sheet is divided into two, and (c) complete mesh refined after the second round
Fig. 11
Iterative mesh refinement using dicing algorithm: (a), (b) the selected sheet is divided into two, and (c) complete mesh refined after the second round
Close modal

If the input face set is re-selected, a new round of refinement is performed using the methods in Secs. 4.1 and 4.2.

5 Results

The proposed method was implemented using C++ programming language based on the ACIS geometric engine. Apart from the example displayed in Fig. 1 (Hereafter referred to as Example 1, the other preliminary experimental results are illustrated in Figs. 12, 13, and 14. The number of hexahedral elements, scaled Jacobian value and running time are presented in Table 1. The comparison of different penalties ci of Examples 2 and 3 is shown in Figs. 15 and 16. The comparison results of different refinement methods are shown in Fig. 18, and the specific statistical information is presented in Table 2. The comparison of different input quads of Example 2 is shown in Fig. 17. The run time was measured on a PC with an Intel Core i7 3.60 GHz CPU and 8GB RAM.

Fig. 12
Example 2: (a) initial coarse mesh, (b) input quads to be refined, (c) quad sets determined using integer optimization, and (d) mesh after sheet inflation
Fig. 12
Example 2: (a) initial coarse mesh, (b) input quads to be refined, (c) quad sets determined using integer optimization, and (d) mesh after sheet inflation
Close modal
Fig. 13
Example 3: (a) geometric model, (b) initial coarse mesh, (c) input quads to be refined, (d) quad sets determined using integer optimization, and (e) mesh after sheet inflation
Fig. 13
Example 3: (a) geometric model, (b) initial coarse mesh, (c) input quads to be refined, (d) quad sets determined using integer optimization, and (e) mesh after sheet inflation
Close modal
Fig. 14
Example 4: (a) geometric model, (b) initial coarse mesh, (c) input quads need to be refined, (d) quad sets determined using integer optimization, and (e) mesh after sheet inflation
Fig. 14
Example 4: (a) geometric model, (b) initial coarse mesh, (c) input quads need to be refined, (d) quad sets determined using integer optimization, and (e) mesh after sheet inflation
Close modal
Fig. 15
Example 2 for different singularity penalty (ci) values: (a) input quads to be refined, (b), (c) quad sets determined and the hex mesh after sheet inflation when ci = 0.01, (d) quad sets determined, and (e) the hex mesh after sheet inflation when ci = 10
Fig. 15
Example 2 for different singularity penalty (ci) values: (a) input quads to be refined, (b), (c) quad sets determined and the hex mesh after sheet inflation when ci = 0.01, (d) quad sets determined, and (e) the hex mesh after sheet inflation when ci = 10
Close modal
Fig. 16
Example 3 for different singularity penalty (ci) values: (a) input quads to be refined, (b), (c) quad sets determined and the hex mesh after sheet inflation when ci = 0.01, (d) quad sets determined, and (e) the hex mesh after sheet inflation when ci = 10
Fig. 16
Example 3 for different singularity penalty (ci) values: (a) input quads to be refined, (b), (c) quad sets determined and the hex mesh after sheet inflation when ci = 0.01, (d) quad sets determined, and (e) the hex mesh after sheet inflation when ci = 10
Close modal
Fig. 17
Example 2 for different input quads: (a) original mesh, (b), (f) two input quads with different degrees of topological smoothness, (c), (g) two quad sets determined for sheet inflation, (d), (h) two sheets after sheet inflation, and (e), (j) two result hexahedral mesh
Fig. 17
Example 2 for different input quads: (a) original mesh, (b), (f) two input quads with different degrees of topological smoothness, (c), (g) two quad sets determined for sheet inflation, (d), (h) two sheets after sheet inflation, and (e), (j) two result hexahedral mesh
Close modal
Fig. 18
Comparison of different refinement methods: (a) original mesh, (b) result mesh after using dicing method, (c) result mesh after using pillowing method, (d) result mesh after using template-based refinement method, and (e) result mesh after using graph-cut method
Fig. 18
Comparison of different refinement methods: (a) original mesh, (b) result mesh after using dicing method, (c) result mesh after using pillowing method, (d) result mesh after using template-based refinement method, and (e) result mesh after using graph-cut method
Close modal
Table 1

Mesh quality and running time of results

Example 1Example 2Example 3Example 4
Hex number338∖518504∖6993930∖44524578∖5503
(before∖after)
Scaled Jacobian of original mesh0.83∖0.120.93∖0.490.82∖0.210.95∖0.15
(mean∖min)
Scaled Jacobian of refined mesh0.82∖0.190.93∖0.470.81∖0.220.95∖0.21
(mean∖min)
Running time (s)8.15.430.2627.2
Example 1Example 2Example 3Example 4
Hex number338∖518504∖6993930∖44524578∖5503
(before∖after)
Scaled Jacobian of original mesh0.83∖0.120.93∖0.490.82∖0.210.95∖0.15
(mean∖min)
Scaled Jacobian of refined mesh0.82∖0.190.93∖0.470.81∖0.220.95∖0.21
(mean∖min)
Running time (s)8.15.430.2627.2
Table 2

Statistics for Example 2 mesh results using different refinement methods

MethodHex numberScaled Jacobian (mean∖min)Running time (s)
Dicing [18]10590.432∖0.9441.5
Pillowing [5]780−0.52∖0.571.4
Template-based [19]14500.33∖0.843.6
Graph-cut [9]6910.42∖0.8752.1
Ours6990.47∖0.935.4
MethodHex numberScaled Jacobian (mean∖min)Running time (s)
Dicing [18]10590.432∖0.9441.5
Pillowing [5]780−0.52∖0.571.4
Template-based [19]14500.33∖0.843.6
Graph-cut [9]6910.42∖0.8752.1
Ours6990.47∖0.935.4

As shown in Fig. 12, the input mesh of Example 2 was generated using the sweeping method. The input quads marked for inflation are along the outermost layer of the mesh boundary, as this is convenient for users to specify. In addition, based on the relatively little input information, topologically smooth quad sets are determined through the addition of the penalty term for introducing singularities into the integer program formulation.

The input meshes of Examples 3 and 4 were generated using the block decomposition method in Ref. [20]. There are fewer singularities in the original mesh, resulting in nonuniform density. Appropriate mesh refinement is performed in the coarse region to improve the scaled Jacobian value of the resulting mesh.

We compare different values of parameter ci with regard to the smoothness of quad set as shown in Figs. 15 and 16. The results show that larger parameter ci values result in the introduction of fewer singularities and smoother quad sets. Therefore, if better quality refined meshes are required, the value of parameter ci should be increased. In addition, the choice of ci mainly depends on the number of mesh elements and the number of input quads. When the number of input quads accounts for a large proportion of the whole model, ci should be increased as much as possible to reduce the introduction of singularities. In general, in order to get better quality refinement results, we tend to set ci to a number greater than ten based on statistical data to ensure that no singularities are inserted as much as possible. When ci is set to ten and the number of inserted singularities is more than expected, it is necessary to increase the value of ci. When ci is set to ten and the number of inserted elements is more than expected, it is necessary to reduce the value of ci.

Figure 17 shows two meshes after refinement for two types of input quads with different topologies. As shown in Figs. 17(e) and 17(j), the rougher input quads, the more singularities are introduced.

In addition, we compare our result with four other hexahedral refinement methods. These methods cannot control the increased number of elements and singularities as shown in Fig. 18. The dicing algorithm [18] can only duplicate the existing sheet which lacks flexibility, resulting in too fine refinement. The pillowing algorithm [5] is often used to optimize the hexahedral mesh on the boundary. As shown in Fig. 18(c), the poor elements are generated because the shape of shrink sets inside the model is too irregular. The template-based refinement method [19] cannot avoid generating singularities in many local regions because of some certain templates. Besides, determining quad set based on graph-cut method [9] cannot deal with input quads with intersections, and the resulting mesh in Fig. 18(e) does not contain the sheets generated by intersecting input quads. As shown in Table 2, our algorithm has advantages in the control of introduced elements number and mesh quality, but time efficiency is a limitation. The limitation comes from the fact that solving binary problem is in general NP-hard.

6 Conclusions

In this paper, a binary program formulation is presented for the purpose of local hexahedral mesh refinement. Compared with existing methods, the innovation of this paper lies in the following aspects:

  1. By introducing user-defined costs in integer program formulation, the completion of quad sets for sheet inflation with controllable quality is realized. Quality control is reflected in the balance between the number of inserted elements and the number of inserted singularities.

  2. By selecting the optimal direction of sheet inflation, the self-intersecting sheet is optimally inserted to ensure the quality of hexahedral mesh.

However, the proposed method also has some problems that need to be solved in the future work. First, the maximal edge valence or maximal self-intersections number of quads should be constrained. Second, the boundary quads could be allowed in the quad sets to refine the mesh boundary by sheet inflation. Third, we may limit the local area for refinement to alleviate the problem that the solution is too slow due to the large number of mesh elements. Finally, the mesh quality after mesh refinement should not be decreased and we will research the further optimization of quad set to ensure the quality of mesh refinement.

Acknowledgment

The authors are very grateful to the financial support from NSF of China (Grant Nos. 62202417, 61802211, and 61902099) and the Open Project Program of the State Key Lab of CAD&CG (Grant No. A2208), and Zhejiang University.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

References

1.
Zhang
,
H.
,
Zhao
,
G.
, and
Ma
,
X.
,
2007
, “
Adaptive Generation of Hexahedral Element Mesh Using An Improved Grid-Based Method
,”
Comput. Aided Des.
,
39
(
10
), pp.
914
928
.
2.
Pietroni
,
N.
,
Campen
,
M.
,
Sheffer
,
A.
,
Cherchi
,
G.
,
Bommes
,
D.
,
Gao
,
X.
,
Scateni
,
R.
,
Ledoux
,
F.
,
Remacle
,
J.-F.
, and
Livesu
,
M.
,
2022
, “
Hex-Mesh Generation and Processing: A Survey
,” arXiv preprint arXiv:2202.12670.
3.
Tchon
,
K.-F.
,
Dompierre
,
J.
, and
Camarero
,
R.
,
2002
, “
Conformal Refinement of All-Quadrilateral and All-Hexahedral Meshes According to An Anisotropic Metric
,”
11th International Meshing Roundtable
,
Ithaca, NY
,
Sept. 15–18
.
4.
Benzley
,
S. E.
,
Harris
,
N. J.
,
Scott
,
M.
,
Borden
,
M.
, and
Owen
,
S. J.
,
2005
, “
Conformal Refinement and Coarsening of Unstructured Hexahedral Meshes
,”
ASME J. Comput. Inf. Sci. Eng.
,
5
(
4
), pp.
330
337
.
5.
Mitchell
,
S. A.
, and
Tautges
,
T. J.
,
1995
, “
Pillowing Doublets: Refining a Mesh to Ensure That Faces Share at Most One Edge
,” Technical Report,
Sandia National Lab.(SNL-NM)
,
Albuquerque, NM
.
6.
Borden
,
M. J.
,
Benzley
,
S. E.
,
Mitchell
,
S. A.
,
White
,
D. R.
, and
Meyers
,
R. J.
,
2000
, “
The Cleave and Fill Tool: An All-Hexahedral Refinement Algorithm for Swept Meshes.
Proceedings of the 9th International Meshing Roundtable
,
New Orleans, LA
,
Oct. 2–5
, IMR, pp.
69
76
.
7.
Tchon
,
K.-F.
, and
Dompierre
,
J.
,
2004
, “
Automated Refinement of Conformal Quadrilateral and Hexahedral Meshes
,”
Int. J. Numer. Methods Eng.
,
59
(
12
), pp.
1539
1562
.
8.
Park
,
C.-H.
, and
Yang
,
D.-Y.
,
2006
, “
Adaptive Refinement of All-Hexahedral Elements for Three-Dimensional Metal Forming Analysis
,”
Finite Elements Anal. Des.
,
43
(
1
), pp.
22
35
.
9.
Shen
,
C.
,
Gao
,
S.
, and
Wang
,
R.
,
2022
, “
Hexahedral Mesh Adaptation Based on Posterior-Error Estimation
,”
Eng. Comput.
10.
Schneiders
,
R.
,
Schindler
,
R.
, and
Weiler
,
F.
,
1996
, “
Octree-Based Generation of Hexahedral Element Meshes
,”
Proceedings of the 5th International Meshing Roundtable
,
Pittsburgh, PA
,
Oct. 10–11
.
11.
Ito
,
Y.
,
Shih
,
A. M.
, and
Soni
,
B. K.
,
2009
, “
Octree-Based Reasonable-Quality Hexahedral Mesh Generation Using a New Set of Refinement Templates
,”
Int. J. Numer. Methods Eng.
,
77
(
13
), pp.
1809
1833
.
12.
Wada
,
Y.
, and
Okuda
,
H.
,
2002
, “
Effective Adaptation Technique for Hexahedral Mesh
,”
Concurr. Comput. Pract. Exp.
,
14
(
6–7
), pp.
451
463
.
13.
Kwak
,
D.-Y.
, and
Im
,
Y.-T.
,
2002
, “
Remeshing for Metal Forming Simulations-Part II: Three-Dimensional Hexahedral Mesh Generation
,”
Int. J. Numer. Methods Eng.
,
53
(
11
), pp.
2501
2528
.
14.
Parrish
,
M.
,
Borden
,
M.
,
Staten
,
M.
, and
Benzley
,
S.
,
2008
, “
A Selective Approach to Conformal Refinement of Unstructured Hexahedral Finite Element Meshes
,”
Proceedings of the 16th International Meshing Roundtable
,
Seattle, WA
,
Oct. 14–17
, Springer, pp.
251
268
.
15.
Zhang
,
Y.
,
Liang
,
X.
, and
Xu
,
G.
,
2013
, “
A Robust 2-Refinement Algorithm in Octree Or Rhombic Dodecahedral Tree Based All-Hexahedral Mesh Generation
,”
Comput. Methods. Appl. Mech. Eng.
,
256
, pp.
88
100
.
16.
Gurobi Optimization
,
2022
,
Gurobi Optimizer Reference Manual
, 9.5 ed.,
6
.
17.
Wang
,
R.
,
Ding
,
M.
,
Wu
,
H.
,
Zheng
,
Z.
, and
Gao
,
S.
,
2020
, “
Optimization Design Method of Complex Sheet Insertion Operation of Hexahedral Mesh
,”
J. Comput. Aided Des. Comput. Graph.
,
32
(
5
), pp.
846
856
.
18.
Melander
,
D. J.
,
Benzley
,
S. E.
, and
Tautges
,
T.
,
1997
, “
Generation of Multi-Million Element Meshes for Solid Model-Based Geometries: The Dicer Algorithm
,” Technical Report,
Sandia National Lab.(SNL-NM)
,
Albuquerque, NM
.
19.
Owen
,
S. J.
,
Shih
,
R. M.
, and
Ernst
,
C. D.
,
2017
, “
A Template-Based Approach for Parallel Hexahedral Two-Refinement
,”
Comput. Aided Des.
,
85
, pp.
34
52
.
20.
Wang
,
R.
,
Shen
,
C.
,
Chen
,
J.
,
Wu
,
H.
, and
Gao
,
S.
,
2017
, “
Sheet Operation Based Block Decomposition of Solid Models for Hex Meshing
,”
Comput. Aided Des.
,
85
, pp.
123
137
.