Development of Octree-Based High-Quality Mesh Generation Method for Biomedical Simulation

Abstract

This paper proposes a robust high-quality finite element mesh generation method which is capable of modeling problems with complex geometries and multiple materials and suitable for the use in biomedical simulation. The previous octree-based method can generate a high-quality mesh with complex geometries and multiple materials robustly allowing geometric approximation. In this study, a robust mesh optimization method is developed combining smoothing and topology optimization in order to correct geometries guaranteeing element quality. Through performance measurement using sphere mesh and application to HTO tibia mesh, the validity of the developed mesh optimization method is checked.

1 Introduction

Numerical simulations of biomedical bodies are expected to bring new knowledge to medicine and bioengineering field through discovering pathology developing mechanisms and optimizing surgical procedures. In particular, biomedical simulations based on finite element (FE) analysis which can analyze arbitrarily shaped field have been an active research topic recently. Thanks to the emergence of medical image diagnostic apparatus such as CT scans and MRI scans, it has become relatively easy to acquire detailed geometric information of biomedical bodies including inner structures without dissection. For the realization of reliable biomedical FE simulations, the mesh generation of biomedical bodies, which have complex geometries and multiple materials generally came to be a large current bottleneck and development of robust mesh generation method which is capable of complex geometries and multiple materials is required.

Though some meshing methods tailored for specific problems such as modeling of soils [1], realized robust mesh generation, they are not capable of complex geometries and multiple materials. As a prototype of a robust generation method of mesh consisting of complex geometries and multiple materials, octree-based fast mesh generation method which can generate global mesh by tetrahedralizing cubes locally using look up table [2, 3] has been developed. This method has high affinity with a mesh generation of biomedical bodies thanks to the ability to use brightness values of two dimensional images as an input. However, mesh generated by this method contains small geometrical approximation, which leads to the bad convergence and local concentration of feature amount, such as stress and strain in analysis results. In order to use this method in biomedical simulations which requires high accuracy, modification of geometrical approximation remains to be settled.

Geometrical modification needs to be performed guaranteeing that invalid elements (i.e., inside out elements and greatly distorted elements) will not be generated which are unacceptable in numerical analysis. Therefore, simple geometrical modification which smooths only boundary nodes is not satisfactory due to the generation of invalid elements especially near boundary surfaces. Frequently boundary nodes get locked because of the restriction of the element quality guarantee and then geometrical modification does not progress any more. To escape from this situation topology optimization (i.e., changing node-element connectivity) is useful and it has been reported that combination of smoothing and topology optimization leads to generation of high quality mesh.

In this study, with the objective of the use in biomedical FE simulations, a robust mesh optimization method to modify geometric approximation guaranteeing mesh quality is developed. Consequently, combining with the previous method, a robust generation method of high-quality mesh whose geometries follow input geometries with multiple materials is established. Mesh optimization method is developed combining smoothing and topology optimization as previous studies. While modifying surface geometries of mesh by smoothing, situation that boundary nodes get locked because of generated bad elements is avoided by topology optimization.

The rest of this paper is organized as follows. Previous mesh generation method and mesh optimization method developed in this study are explained in Sect. 2. Performance measurement of developed method is performed in Sect. 3. Developed method is applied to tibia model after HTO surgery in Sect. 4. Section 5 summarizes this paper.

2 Methodology

Mesh generation method developed in this study follows the procedures shown in Fig. 1. In the first half, initial mesh is generated using octree-based mesh generation method. By allowing small geometric approximation, boundary surfaces of multiple materials are explicitly resolved with robustness, mesh conformity is ensured and mesh quality is guaranteed. However, geometric approximation, which leads to bad convergence and local concentration of feature amount in analysis results is undesirable. In the second half, mesh optimization reduces the magnitude of geometric approximation. This mesh optimization method repeats smoothing and topology optimization and thus the geometric approximation is reduced according to the number of iterations. In the following subsections, we explain the overview of octree-based mesh generation method and mesh optimization method.

Fig. 1.

Flowchart of developed mesh generation method

2.1 Overview of Octree-Based Mesh Generation Method

Using closed triangular patches discretizing surface information of each material region as an input, this method generates conforming mesh with linear tetrahedral elements robustly. The target domain is decomposed into multi-scale cubes using an orthogonal octree structure and then each cube is decomposed into tetrahedral elements. The former process is referred to as “cube generation” and the latter process as “cube decomposition” in this paper.

Cube Generation. Using closed triangular patches defining each material region as an input, the target domain (Fig. 2a) is recursively decomposed into multi-scale cubes by orthogonal octree structure (Fig. 2b). While fine cubes are allocated near material boundaries in order to resolve detailed geometries of boundary surfaces, coarse cubes are allocated inside material regions to generate mesh with minimum degree of freedoms. The size of boundary cubes, which have input material boundaries inside and are allocated near material boundaries is the smallest and referred to as “resolution” in this paper. Since large gap in size between neighboring cubes will directly lead to generation of bad quality tetrahedral elements and deterioration in analysis accuracy, size difference between neighboring cubes is restricted to be within one level. At the same time, sizes of cubes inside material regions are restricted to be smaller than the size prescribed for each material to generate fine enough mesh required for numerical analysis. Each node is allocated material information which identifies the material region the node belongs to.

Fig. 2.

Cube generation procedure in 2D

Cube Decomposition. Each cube is decomposed into tetrahedral elements referring to look up table using material information of nodes as an input (Fig. 3a–b). Mesh generation resolving multiple material boundaries is enabled by constructing boundary cube decomposition look up table using multiple material marching cubes method [4] which can define boundary surfaces of multiple materials inside structured grid. Moreover, by storing cube decomposition patterns in look up table which guarantee mesh conformity between neighboring cubes, tetrahedralization of each cube is realized.

Following procedures explained above, a conforming linear tetrahedral mesh with complex geometries and multiple materials is generated robustly. Uniqueness of all procedures involved in cube generation and decomposition assures the robustness in mesh generation. Though originally the input is material region definition in 3D, it is possible to allocate material information of each node from brightness values of image information in 2D. Therefore, it is not difficult to modify current source code to use 2D image as an input and this method has high affinity with mesh generation of biomedical bodies. However, mesh generated by this method contains geometric approximation whose size is half of the resolution at maximum. This is due to the location restriction of new boundary nodes to lattice points in 3D when decomposing boundary cubes with input material boundaries inside and defining new material boundaries (Fig. 4). Though guaranteeing mesh conformity and quality, this geometric approximation is undesirable since bad convergence and local concentration of feature amount will be brought together in FE analysis. In next subsection, mesh optimization method to reduce geometric approximation is developed.

Fig. 3.

Cube decomposition procedure

Fig. 4.

Geometric approximation

2.2 Mesh Optimization

Overview of Mesh Optimization. Mesh optimization can be classified into to two groups, smoothing and topology optimization according to the modification content of mesh information. The former smoothing changes positions of nodes while preserving mesh topology, or node-element connectivity. The latter topology optimization changes mesh topology while preserving node positions. Smoothing, which repositions node is available for geometric modification, but in many cases nodes get locked because of the generation of bad quality elements during the iteration of smoothing. Quality of elements with boundary nodes tends to be bad and in particular, it often becomes impossible to improve quality of elements with four boundary nodes by smoothing [5]. In such cases, topology optimization is effective to improve element quality and reportedly combination of smoothing and topology optimization leads to mesh with better quality [6, 7].

On the other hand, mesh optimization can be classified into global optimization and local optimization according to the size of domain to optimize in one process. While topology optimization is always performed locally since global topology optimization means global mesh generation, smoothing can be performed either locally or globally. Global smoothing, which repositions all nodes in whole mesh which are free to move at the same time is reduced to solve constrained non-linear programming problem in order to move nodes guaranteeing element quality firmly. However, in addition to the difficulty of the formulation of the objective function, obtaining global solution requires large computation and often fails since there is no guarantee of convergence in general cases. In local smoothing, which moves all nodes inside each small domain which are free to move at the same time, the optimization of the whole mesh is achieved by solving local optimization problems in each small domain. While having advantages such as the easy convergence of each local optimization which requires small computation and the suitability for parallelization, local optimization has disadvantages such as the increase of the overall computation because of the iteration of local optimization in each small domain [8] and convergence to local solution depending on the objective function [8].

In this study, mesh optimization method combining smoothing and topology optimization is developed following previous studies. Since topology on boundary surfaces is maintained while nodes are repositioned and element topology inside material regions is changed, mesh fineness on surface is preserved. Developed method is designed to modify geometry by iterating unit process consisting of smoothing and topology optimization. Smoothing improves surface geometries of mesh while topology optimization avoids such a situation that nodes get locked because of bad quality elements. Both smoothing and topology optimization in unit process are preformed locally and iteratively in many small domains. By adopting local processes which have small computation for smoothing and topology optimization, steady geometric improvement guaranteeing that no element violates the predetermined value of element quality metric. Moreover, room is left for acceleration envisaging the application to a large-scale mesh by adopting local processes which can be naturally parallelized.

In the following subsections, preparation required for mesh optimization is explained first, and then boundary node smoothing, interior node smoothing, and topology optimization of low quality elements which consists of unit process of mesh optimization are explained.

Preparation. First, as a preparation for mesh optimization, construction of boundary cube decomposition look up table with two materials which does not contain any element consisting of four boundary nodes and assignment of projection destinations to boundary nodes. In smoothing, tetrahedral elements consisting of four boundary nodes are undesirable. Since boundary nodes are moved toward projection destinations, quality of tetrahedral elements sharing boundary nodes tend to deteriorate and especially tetrahedral elements consisting of four boundary nodes should not exist. Therefore, in this study, decomposition look up table of boundary cube with two materials, which is frequently refereed to when generating initial mesh is improved to have no tetrahedral element consisting of four boundary nodes. Additionally it is required to choose projection destinations for each boundary node from input triangle patches defining material surfaces in advance of mesh optimization processes. This projection destination information will be used to move boundary nodes toward desirable position to improve mesh geometry. As for the assignment of projection destination of edge-centered boundary node, one triangle patch which has intersection with the edge on which the edge-centered node exists is chosen from input triangle patches. Projection destinations of face-centered boundary node is set to be the union of the projection destinations of the edge-centered boundary nodes which are on the same face. Therefore, face-centered boundary node can have four projection destinations at maximum. Based on the same idea, projection destinations of cube-centered boundary node is set to be the union of the projection destinations of the edge-centered boundary nodes which are on the same cube. Cube-centered boundary node can have twelve projection destinations at maximum.

Boundary Node Smoothing. Each boundary node has its own projection destination information. In this process, after Laplacian smoothing, boundary nodes are moved toward their projection destinations. Since the void domain of the target problem is also considered as one material, boundary nodes which are shared with air part and another material are also moved.

• Boundary node Laplacian smoothing

Position of boundary nodes is improved to mitigate the distortion of elements brought by boundary node repositioning. This process is essential to avoid such a situation that boundary nodes get locked because of the approach of boundary nodes in the iteration of mesh optimization. Laplacian smoothing is adopted as a determining method of new position on the ground of the easiness of implementation and small computation. However, simple movement to the arithmetic mean of all neighbouring nodes of the boundary node may result in further position from the projection destinations because of the neighbouring interior nodes. This phenomena will be remarkable with boundary nodes sharing air part (Fig. 5a–b). Therefore, in this study, boundary nodes are moved to the arithmetic mean of all neighbouring boundary nodes (Fig. 5a–c). Whether inside out element or low quality element is generated by this node repositioning is checked at every step to deal with the disadvantage of Laplacian smoothing. In case these invalid elements are generated, boundary nodes are moved by half of the movement vector. Thus, the boundary node Laplacian smoothing procedure is summarized as:

• 1. Calculate the arithmetic mean of neighbouring boundary nodes and movement vector from the current position.

• 2. Calculate new position by adding movement vector to the current position.

• 3. When moved to new position, if inside out element or low quality element whose element quality metric value violates predetermined value is generated, half the movement vector and go to 2.

• 4. Move to new position.

Note that this boundary node Laplacian smoothing is performed only in the first half of the mesh optimization iterations since it does not always move each boundary node toward its own projection destination.

• Movement of boundary node toward boundary surface

Boundary nodes are moved to boundary surfaces designated by projection destination information. New position is set to be the foot of a perpendicular line projected from the current position to the boundary surface (Fig. 5c–d). In case these invalid elements such as inside out elements or low quality elements are generated, boundary nodes are moved by half of the movement vector. Since there were some cases that a element get stuck in other element, the movement vector is also halved in these cases. Movement vector of boundary node which has two or more projection destinations is set to be the arithmetic mean of all movement vectors toward each boundary surfaces. Thus, movement procedure of boundary node toward boundary surface is summarized as:

• 1. Calculate movement vector from the current position to the foot of a perpendicular line on the boundary surface (Fig. 6). If boundary node has two or more projection destinations, calculate the arithmetic mean of all movement vectors toward each boundary surfaces.

• 2. Calculate new position by adding movement vector to the current position.

• 3. When moved to new position, if inside out element or low quality element whose element quality metric value violates predetermined value is generated or a element get stuck in other element, half the movement vector and go to 2.

• 4. Move to new position.

However, sometimes the above procedure cannot move boundary nodes sufficiently. One possibility of this phenomena is because of the activation of element quality restriction since there are other boundary nodes which share the boundary surface with the boundary node to move on the way to the new position. This problem stems from the fact that movement vector is simply set to be in the direction of the foot of a perpendicular line on the boundary surface. Therefore, when the above procedure failed to work, some perturbation is added to the movement vector and movement in the direction of boundary surface avoiding other boundary nodes is enabled. Perturbation vector is set by random number and its magnitude is set to be half of the magnitude of the movement vector. Since the magnitude ratio of movement vector and perturbation vector is constantly 1:2, new position is assured to be closer to the boundary surface than the current position. Thus, movement procedure of boundary node toward boundary surface with perturbation is summarized as:

• 1. Calculate movement vector from the current position to the foot of a perpendicular line on the boundary surface. If boundary node has two or more projection destinations, calculate the arithmetic mean of all movement vectors toward each boundary surfaces.

• 2. Set perturbation vector whose magnitude is half of the magnitude of movement vector using random number.

• 3. Calculate new position by adding movement vector and perturbation vector to the current position.

• 4. When moved to new position, if inside out element or low quality elements whose element quality metric value violates predetermined value is generated or a element get stuck in other element, half the movement vector and go to 3.

• 5. Move to new position.

Fig. 5.

Example of mesh near air part. Small dots represent boundary nodes. The boundary node represented by the red dot in (a) is going to be repositioned. While simple Laplacian smoothing results in (b) pulled by interior nodes, developed method results in (c) near the boundary surface. Next, projection to boundary surface results in (d). (Color figure online)

Fig. 6.

Movement vector toward boundary surface

Interior Node Smoothing. This process improves position of interior nodes. The objective of this process is to relieve the distortion of elements near boundary surfaces which is brought by moving boundary nodes to boundary surfaces in entire mesh. By repositioning interior nodes and mitigating the distortion near boundary surfaces, further movement of boundary nodes to boundary surfaces is intended. New position of interior nodes is determined by Laplacian smoothing and in case of the generation of invalid elements, interior nodes are moved by half of the movement vector. Thus, smoothing procedure of interior node is summarized as:

• 1. Calculate the arithmetic mean of neighbouring nodes and movement vector from the current position.

• 2. Calculate new position by adding movement vector to the current position.

• 3. When moved to new position, if inside out element or low quality element whose element quality metric value violates predetermined value is generated, half the movement vector and go to 2.

• 4. Move to new position.

Topology Optimization. Topology of elements whose quality became almost no better than the predetermined value in the iteration of mesh optimization is changed. Topology optimization involving more elements tends to bring about better results since the number of patterns to change topology increases. Therefore, in this study, the algorithm of topology optimization is designed to change topology of several elements at the same time. First, an interior node composing low quality element is chosen. The reason why boundary node is not chosen is to maintain the topology of material boundary surfaces. Next, the surface triangles of the polygon consisting of tetrahedral elements which share the interior node is extracted. This process deletes the interior node and all tetrahedral elements sharing this interior node. Finally, using constrained Delaunay tetrahedralization, the polygon is decomposed into tetrahedral elements preserving the surface triangles. New nodes generated in this process are restricted inside the polygon. However, if a tetrahedral element consisting of four boundary nodes or a low quality tetrahedral element is generated, topology optimization is not performed. Advantages of these processes are the topology preservation on boundary surfaces and the tendency to improve the element quality by changing the topology of several elements at the same time. Thus, topology optimization procedure is summarized as:

• 1. Choose an interior node composing a low quality element.

• 2. Extract all tetrahedral elements sharing the interior node and construct polygon

• 3. Extract surface triangles from the polygon.

• 4. Perform constrained Dealunay tetrahedralization preserving surface triangles.

• 5. If and only if no tetrahedra consisting of four boundary nodes or no low quality tetrahedra whose element quality metric value violates predetermined value is generated, change topology.

3 Performance Measurement

Performance of mesh optimization method developed in previous section is measured. Considering that geometric approximation can become large on curved surfaces, sphere with a diameter of 8 is used. Figure 8a shows initial sphere mesh with resolution 0.1 before mesh optimization. The mesh consists of 452,280 linear tetrahedral elements, 90,813 nodes, 272,439 degrees of freedom, and 30,070 boundary nodes. Geometric approximation is clearly visible. On the other hand, sphere mesh after mesh optimization is shown in Fig. 8b, which has smoother surface. In this mesh optimization, the total number of iterations consisting of smoothing and topology optimization set 50, while Laplacian smoothing of boundary nodes is performed only in the first 25 iterations. Allowable maximum aspect ratio is set 30. Aspect ratio, a, represents an indicator of element quality and an element with a smaller a value is considered to exhibit better quality (a becomes 1 for a regular tetrahedron). Figure 9a shows maximum distance and average distance between each node and boundary surface on which the node should be in terms of iterations. Before mesh optimization shown as 0th iteration, maximum distance is 4.96E−2 which almost match the half size of resolution 0.1 and average distance is 1.67E−2. While distance increases in some part of iterations until 25th iteration, after 26th iteration distance monotonically decreases since Laplacian smoothing of boundary nodes is not performed and after finishing 50th iteration both maximum and average distance are small enough. Moreover, after 38th iteration, maximum and average distance almost converged to 4.25E−15 and 1.9E−16 respectively. Figure 9b shows cumulative relative frequency distribution of element aspect ratio inside overall mesh. Max aspect ratio changed from 4.87 to 30 through mesh optimization. It can be understood that nodes moved to the limit of allowable maximum aspect ratio. Even after mesh optimization, 95% of all elements have aspect ratio smaller than 10 and overall mesh quality is still maintained (Fig. 9).

Fig. 7.

Topology optimization of low quality element. Change topology of elements sharing the interior node represented with the black dot.

Fig. 8.

Mesh optimization of sphere mesh (Res. 0.1)

Fig. 9.

Performance measurement results

4 Application Example

High tibial osteotomy (HTO) is a surgical procedure to disperse excessive load on the inner side of the thigh due to bow leg deformation toward the outer side of the thigh. For the realization of less invasive HTO surgery, quantitative evaluation of stress distribution applied to the bone, such as presence of stress concentration, is considered to be effective. This time, using two tibia meshes after HTO with and without mesh optimization respectively, linear elastic analysis under static load was performed and the validity of the developed mesh generation method was verified by comparing stress distributions between them.

The tibia of the right leg of an adult male was used [9, 10]. Setting the resolution to 0.05 cm, tibia mesh consisting of 5,459,851 linear tetrahedral elements, 1,073,949 nodes and 3,221,847 degrees of freedom was generated (Fig. 10). Initial mesh generation took 157 s using 64 cores and mesh optimization took 7,947 s using 1 core of SGI UV 300. SGI UV300 is a 512 core, 24.5-TB cache-coherent shared memory system consisting of 32 sockets of 16-core Intel Xeon E7-8867 v3 CPUs and 768 slots of 32 GB DDR4 DRAM connected with SGI NUMAlink ASIC technology. Smooth surface of the generated mesh can be seen.

Next, linear elastic structural analysis was performed using generated tibia mesh. As a boundary condition, static load was applied on the top face of tibia in downward direction and lower half part was fixed. Comparison of minimum principal stress distributions between two meshes is shown in Fig. 11. While contour line on the mesh before mesh optimisation is rough, contour line on the mesh after mesh optimisation is smooth. It is suggested that mesh optimization method developed in this study reduces local concentration of stress and enables quality assurance of numerical analysis by mesh refinement. Though the inner material of tibia used in this analysis was assumed to be uniform for simplicity, making use of the ability to resolve multiple material boundaries of developed mesh generation method, analysis considering heterogeneous material distribution inside bone would be easily realizable.

Fig. 10.

HTO tibia mesh after mesh optimization (cm)

Fig. 11.

Comparison of maximum principal stress distribution

5 Closing Remarks

In this study, by developing mesh optimization method which modifies geometric approximation of previous mesh generation method, robust generation method of high-quality mesh with complex geometries and multiple materials is established. Steady geometry modification guaranteeing element quality is realized by iterating local smoothing and topology optimization. Performance measurement using sphere mesh confirmed that geometry is finely improved. Comparing stress distributions of tibia meshes, it is suggested that mesh optimization method realizes numerical simulations whose quality can be assured by mesh refinement. One of the future works would be the acceleration of mesh optimization method. Current mesh optimization takes some days to complete optimization of mesh with $O(10^7)$ elements. Considering that initial mesh with $O(10^9)$ elements can be generated in several hours, by accelerating mesh optimization method whose room for parallelization possibility is left, mesh generation for large-scale biomedical simulations would be easily possible.