## 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 [4–9] and those based on templates [1,10–15]. 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 [11–15] 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 *F*_{O} (*F*_{O} 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:

Determination of quad sets to be inflated by integer optimization.

Overall mesh refinement by optimized sheet inflation.

## 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.

*f*

_{i}of mesh, a binary variable

*s*

_{i}is used to indicate whether the face is to be inflated or not. So the refinement function can also be defined as

*R*(

*f*

_{i}) =

*s*

_{i}. One of our objective for integer optimization is

*E*

_{s}to minimize the number of selected quads

*F*

_{inner}refers to the set containing all inner quads of mesh.

*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,

*s*

_{i}for each quad in

*P*is set to 1

*a*for every inner edge

*e*∈

*E*

_{inner}, the constraint can be expressed as follows:

*s*

_{i}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).

*p*

_{i}is used to denote every possible situation leading to singularities. Therefore, when singularities occur, the

*p*

_{i}is set to 1. The constraint can be expressed as follows:

*s*

_{0}and

*s*

_{1}correspond to two adjacent quads which belong to

*F*

_{inner}.

*c*

_{i}is specified by the user and the another objective of our integer optimization is set as follows:

## 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:

*Given a hexahedral mesh**M* = (*H*, *F*, *E*, *V*) *and an input quad set**Q* = (*F*_{q}, *E*_{q}, *V*_{q})*, for every node**v**in**V*_{q}, *Q**divides the hexahedron set**H*_{v} ∈ *H**into**n**parts* (*n* = 1, 2, 3, 4). *Each part represents a potential splitting direction of node**v*. *The local shrink set**S*_{v}*of node**v**could include several parts of**H*_{v}. *Figure 4 **shows two cases of local shrink set*.

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

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

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.

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.

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.

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

### 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 *V*_{1} and *V*_{2}, where *V*_{1} is at the original position as shown in Fig. 6(a) and *V*_{2} 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.

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:

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:

where$T1=\u2211e[(ne1\u22122)2+(ne2+nei3)2+(ne4\u22122)2]T2=\u2211e[(ne1\u22121)2+(ne2\u22121)2+(ne3\u22121)2+(ne4\u22121)2]$*e*refers to any mesh edge at the intersection line of quad sets,*T*_{1}refers to the quality score after using the template A as shown in Figs. 7(c)–7(d), and*T*_{2}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.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).

### 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.

## 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 *c*_{i} 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.

Example 1 | Example 2 | Example 3 | Example 4 | |
---|---|---|---|---|

Hex number | 338∖518 | 504∖699 | 3930∖4452 | 4578∖5503 |

(before∖after) | ||||

Scaled Jacobian of original mesh | 0.83∖0.12 | 0.93∖0.49 | 0.82∖0.21 | 0.95∖0.15 |

(mean∖min) | ||||

Scaled Jacobian of refined mesh | 0.82∖0.19 | 0.93∖0.47 | 0.81∖0.22 | 0.95∖0.21 |

(mean∖min) | ||||

Running time (s) | 8.1 | 5.4 | 30.2 | 627.2 |

Example 1 | Example 2 | Example 3 | Example 4 | |
---|---|---|---|---|

Hex number | 338∖518 | 504∖699 | 3930∖4452 | 4578∖5503 |

(before∖after) | ||||

Scaled Jacobian of original mesh | 0.83∖0.12 | 0.93∖0.49 | 0.82∖0.21 | 0.95∖0.15 |

(mean∖min) | ||||

Scaled Jacobian of refined mesh | 0.82∖0.19 | 0.93∖0.47 | 0.81∖0.22 | 0.95∖0.21 |

(mean∖min) | ||||

Running time (s) | 8.1 | 5.4 | 30.2 | 627.2 |

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 *c*_{i} with regard to the smoothness of quad set as shown in Figs. 15 and 16. The results show that larger parameter *c*_{i} values result in the introduction of fewer singularities and smoother quad sets. Therefore, if better quality refined meshes are required, the value of parameter *c*_{i} should be increased. In addition, the choice of *c*_{i} 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, *c*_{i} 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 *c*_{i} to a number greater than ten based on statistical data to ensure that no singularities are inserted as much as possible. When *c*_{i} is set to ten and the number of inserted singularities is more than expected, it is necessary to increase the value of *c*_{i}. When *c*_{i} is set to ten and the number of inserted elements is more than expected, it is necessary to reduce the value of *c*_{i}.

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:

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.

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

## 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

*arXiv preprint arXiv:2202.12670*.

*Gurobi Optimizer Reference Manual*