#### Abstract

In some simulations like virtual surgery, an accurate surface deformation method is needed. Many deformation methods focus on the whole model swing and twist. Few methods focus on surface deformation. For the surface deformation method, two necessary characteristics are needed: the accuracy and real-time performance. Some traditional methods, such as position-based dynamics (PBD) and mass-spring method (MSM), focus more on the real-time performance. Others like the finite element method (FEM) focus more on the accuracy. To balance these two characteristics, we propose a hybrid mesh deformation method for accurate surface deformation based on FEM and PBD. Firstly, we construct a hybrid mesh, which is composed of a coarse volume mesh and a fine surface mesh. Secondly, we implement FEM on coarse volume mesh and PBD on fine surface mesh, and the deformation of fine surface mesh is constrained by the displacement of the coarse volume mesh. Thirdly, we introduced a small incision process for our proposed method. Finally, we implemented our method on a simple deformation simulation and a small incision simulation. The result shows an accurate surface deformation performance by implementing our method. The incision effect shows the compatibility of our proposed method. In conclusion, our proposed method acquires a better trade-off between accuracy and real-time performance.

#### 1. Introduction

The mesh deformation algorithm plays a very important role in computer simulation technique, such as computer animation and soft tissue deformation simulation in virtual reality. In some simulations like virtual surgery, the surface deformation plays an important part in realistic perform effect. Generally speaking, two major problems need to be solved for surface deformation algorithm: low time consumption and high accuracy. Although many methods focus on one of these two problems, there is a huge demand for high accuracy deformation under real-time performance. Therefore, a suitable deformation algorithm needs a trade-off between the accuracy requirement and the real-time performance.

The finite element method (FEM) is usually used for the numerical discrete model of constitutive equation in elastic mechanics [1]. The standard method is used to create physical elements according to each geometric element and then to calculate the displacement of relevant element nodes according to the applied load [2]. The boundary conditions need to be well treated. The FEM can build different kinds of physical elements according to corresponding geometric units, such as triangle element in 2D and tetrahedral element in 3D [3]. The advantage of FEM is that it is close to the real world since it utilizes the actual physical parameters to get output with high accuracy of deformation. It is also able to simulate the complex mesh model [4]. However, FEM costs more computational time. The real-time performance needs to be improved [5]. At present, FEM is widely used in mesh deformation simulations. Paulus et al. proposed a new re-mesh method and rapid FEM for soft tissue simulation [6]. Courtecuisse et al. demonstrated a real-time tissue simulation based on the implicit numerical integration method in nonlinear FEM [7]. Haouchine et al. proposed a method to enhance the internal structure of the liver mesh model in real time [8]. Tang and Wan proposed a constraint FEM for the interactive simulation for virtual surgery [9]. Kugelstadt et al. proposed a fast co-rotated FEM using operator splitting for whole mesh model swinging and twisting [10]. It acquires a good trade-off between accuracy and low time consumption but is not compatible with mesh dissection (model incision).

Position-based dynamics (PBD) has been widely used in recent years. Classical model methods such as FEM and MSM calculate the deformation according to Newton’s second law. Compared with these methods, PBD directly calculates the position change based on the solution of the quasi-static problem. Therefore, PBD has better real-time performance and stability [11]. However, PBD is generally not as accurate as classical model methods. In recent years, PBD has been widely used in deformation simulations. Kubiak et al. proposed a surgical line simulation method based on the PBD method of Müller et al. [12, 13]. In order to assist the robot-assisted rehearsal and planning of partial nephrectomy, Camara et al. proposed a real-time simulation platform that allows surgeons to quickly construct a biomechanical model [14]. Pan et al. proposed an interactive method of the hybrid soft tissue model based on extended PBD [15]. Liu et al. proposed a PBD method for tetrahedral element mesh model to simulate the deformation on soft tissue [16].

For the mass-spring method (MSM), it has many advantages such as easy operation, simple modeling, and low time consumption. However, MSM needs suitable constraints. Unsuitable constraints will lead to fault and divergence of deformation results. In MSM, it is important to determine the spring parameters which include stiffness and damping coefficient. These parameters determine the accuracy of MSM. Due to the random selection of these parameters, the accuracy of deformation results is not high. Two ways of parameter selection are data measurement and formula analysis. Data measurement estimates the parameters by measuring the deformation data of real biological tissues. The formula analytical method is based on the constitutive equation of the material to deduce the parameters [17]. The MSM also has a wide range of applications. Skornitzke et al. proposed a mitral valve MSM model [18]. Gao et al. simulated the nonlinear anisotropy of biological materials based on MSM [19]. Zhang et al. proposed a three-stage method based on MSM to simulate soft tissue [20].

In this paper, we propose a hybrid mesh deformation method for accurate surface deformation based on FEM and PBD. Our proposed method combined the advantages of FEM and PBD to satisfy the requirement of accuracy and real time. First, we construct the hybrid mesh which composes a fine surface mesh and a coarse volume mesh. Then, the FEM and PBD is utilized to simulate deformation on the volume mesh and surface mesh, respectively. Finally, the deformation of surface mesh is constrained on the volume mesh. We also introduced a small incision process which is compatible with our proposed method.

The reminder of this paper is organized as follows: Section 2 briefly explains purpose of our proposed method. Sections 3 and 4 describe the mathematical schemes including FEM and PBD. Section 5 introduces the construction of the hybrid mesh. Section 6 presents the algorithm flow of our proposed method. Section 7 introduces the small incision processing for our proposed method. Section 8 introduces the experiment environments, shows the rendering effects and results of deformation, and demonstrates the compatibility of small incision processing. In Section 9, we conclude our work.

#### 2. Overview

In previous work, FEM is used for high accuracy in fine mesh deformation but usually has high time consumption. PBD has the characteristics of low time consumption and low accuracy. In our method, FEM is used to calculate the deformation of the coarse tetrahedral mesh, and the PBD is used to perform the deformation of the fine surface triangle mesh. The nodes of the volume mesh which connect to the surface mesh are control points (black points in Figure 1). These points constrain the surface mesh and make its deformation within a more accurate range. Meanwhile, since the FEM is applied to the coarse tetrahedral mesh and the surface mesh uses the PBD, the time consumption is not too high. Thus, the real-time requirements can be met.

#### 3. Finite Element Method

Our method uses a set of tetrahedral elements to separate the 3D domain. Based on this spatial discretion, the node displacement interpolation method is used to approximate the continuous displacement field. Specifically, according to equation (1), using the node of elements, the shape function is interpolated to form the displacement field function of the mesh element.

is the shape function matrix of elements, and is the displacement of the mesh element node. We choose linear interpolation as the shape function of the tetrahedron discretion. Through the shape function, the relationship between the strain and the displacement of the element node is established:

is the strain matrix of elements. According to the principle of virtual work, the relationship between nodal displacement and external nodal force of elements **f**^{e} is inwhere **K**^{e} is the stiffness matrix of elements (equation (4)) and **D** is the elastic matrix:

We combine each element stiffness matrix and element external force according to the global index of the mesh nodes to obtain the global stiffness matrix **K** and global external force **f** for the whole mesh. The relationship between global node displacement and global external force vector **f** can be obtained in

When it is necessary, solve equation (5) to calculate the deformation displacement **u**. **K** is determined in preprocessing which costs no time in deformation calculation.

#### 4. Position-Based Dynamics

In PBD, it is necessary to predict the position and velocity of each node and update the nodal displacement to proper position by constraint function. The establishment of the constraint function affects the stability and efficiency of the model. Thus, solving the constraint function becomes the most important part of PBD. Let the constraint function be *C*(**P**) = 0. For the nodal predicted position **P**, the nodal position projected by the constraint functions is **P** + ∆**P**. So, constraint function *C*(**P** + ∆**P**) iswhere ∆**P** is the correction of **P** after the constraint projection and ∇_{P}*C*(**P**) is the gradient of the constraint function. In order to maintain the conservation of momentum and angular momentum, ∆**P** and ∇_{P}*C*(**P**) have the same direction. Thus, a Lagrangian multiplier *λ* can be used to express ∆**P** and ∇_{P}*C*(**P**) relationship:

In our method, we utilize two constraints to perform the deformation of fine surface mesh: stretch constraint and bending constraint. Two constraints are as follows.

##### 4.1. Stretch Constraint

As Figure 2 shows, the stretch constraint simulates the elastic force between two nodes.

The stretch constraint represents equation (8), where *l* is the initial length of **P**_{1}**P**_{2}.

According to equation (7), we obtained ∆**P**_{1} and ∆**P**_{2} in equation (9), where = 1/*m*_{i} (*i* = 1, 2). *m*_{1} and *m*_{2} are the mass of **P**_{1} and **P**_{2}:

##### 4.2. Bending Constraint

As Figure 3 shows, the constraint simulates the bending between two triangle elements.

The bending constraint function iswhere **n**_{1} and **n**_{2} are the normal of ∆**P**_{1}**P**_{3}**P**_{2} and ∆**P**_{1}**P**_{2}**P**_{4} and *φ* is the initial angle between ∆**P**_{1}**P**_{3}**P**_{2} and ∆**P**_{1}**P**_{2}**P**_{4}. According to equation (7), we obtained ∆**P**_{i} (*i* = 1, 2, 3, 4) inwhere = 1/*m*_{i}, *m*_{i} is the mass of **P**_{i}, and **q**_{i} (*i* = 1, 2, 3, 4) is

#### 5. Construction of the Hybrid Mesh

The hybrid mesh includes a fine triangle surface mesh and a matching internal coarse tetrahedral volume mesh. Now, given a fine triangle surface mesh, the steps to generate the hybrid mesh are as follows:(i)Step 1: input the fine triangle surface mesh *M*_{S}(ii)Step 2: simplify *M*_{S} to a coarse triangle surface mesh (iii)Step 3: according to Delaunay algorithm offered in open-source software Tetgen, use to generate a coarse tetrahedral volume mesh *M*_{V}(iv)Step 4: combine *M*_{S} and *M*_{V} to obtain the hybrid mesh

The Illustration of these steps is shown in Figure 4. Since Step 3 uses the open-source algorithm, we do not introduce it here [21]. Now, Step 2 is introduced as follows. Given *M*_{S}, it is necessary to select some nodes from *M*_{S} as the nodes of the coarse surface triangle mesh . The selection of these nodes is based on the principle of minimum vertices cover; i.e., some nodes in *M*_{S} are selected to a node set *V*, so that each edge in *M*_{S} contains at least one node in *V* and the number of nodes in *V* is minimized. We use the method in [22] to obtain the *V*.

After *V* is obtained, the way to simplify *M*_{S} to is introduced in the following. As shown in Figure 5, for one node *p*_{i} in *M*_{S} and not in *V*, a node *p*_{j} which is in *V* and adjacent to *p*_{i} is chosen. After that, *p*_{i} is merged with *p*_{j} and the edge *p*_{i}*p*_{j} is removed. For all adjacent vertices *p*_{k} of *p*_{i}, if *p*_{k} is also adjacent to *p*_{j}, the edges *p*_{i}*p*_{k} are also removed. Otherwise, a new edge *p*_{k}*p*_{j} will be created after *p*_{i}*p*_{k} removed.

Then, the selection of *p*_{j} will become a problem. In order to find a proper *p*_{j}, we first calculate the average distance *d*_{avg} between *p*_{i} and *p*_{k}:where *d*_{ik} is the distance between *p*_{i} and *p*_{k}. And then, we select *p*_{j} that satisfies

The node *p*_{j} is selected in this condition in order to guarantee a better mesh quality after simplification. When all the nodes *p*_{i} are simplified in this way, a coarse surface triangle mesh is obtained. Step 2 is now complete.

#### 6. Deformation Algorithm Based on the Hybrid Mesh

The algorithm includes processing of external force, computing the deformation of the coarse volume mesh (*M*_{V}) using FEM, deformation of *M*_{V}, position prediction of each node on the fine surface mesh (*M*_{S}), and constraint projection on *M*_{S}. Figure 6 shows the illustration of the deformation algorithm on the 2D sectional view of the hybrid mesh. The black points indicate the shared nodes of *M*_{V} and *M*_{S}. These nodes are obtained from the *M*_{S} simplification in Section 5. We call these nodes *control nodes/points*. Control nodes share the same position information with nodes on the boundary of *M*_{V}. The white points represent normal nodes in *M*_{S}. *M*_{V} is represented by dotted lines, which means that it is not rendered in graphics. The solid line at the bottom of *M*_{V} indicates that it is subject to displacement constraints.

##### 6.1. Processing of External Force

Figure 6(a) shows that a vertical upward force is applied to the hybrid mesh. Since *M*_{S} has many nodes, the default load point is on the node. If not, snap the load point to the closest node. After snapping, as Figure 6(b) shows, the position of load point may not be on the control node. In this condition, as Figure 6(c) shows, the equivalent force on the surrounding nodes of the external force is calculated for *M*_{V}. The relationship between equivalent force and external force satisfies linear interpolation as the shape function of the triangular element.

##### 6.2. Deformation of *M*_{V}

After processing of external force, deformation of *M*_{V} is solved according to equation (5). Control nodes share the same position with the nodes on the boundary of *M*_{V}, so the position of control nodes is also obtained. Solving equation (5) will not take too long to run since *M*_{V} is coarse.

##### 6.3. Position Prediction of *M*_{S}

For position prediction of *M*_{S}, like PBD, the node position changes under the external load, without considering the influence of internal forces between the nodes in *M*_{S}. This is node position prediction. The final node position of *M*_{S} will be determined after constraint projection. Take the *M*_{S} in Figure 6(d) as an example, where only one node is affected by the external load. Then, its position changes. The other vertices are not affected by external loads, so their positions have not changed. In this time, the nodes on *M*_{S} are ready for constraint projection.

##### 6.4. Constraint Projection of *M*_{S}

As Figure 6(e) shows, since control nodes share the same position on *M*_{S} and the boundary of *M*_{V}, the displacement of control nodes in *M*_{S} is determined when deformation of *M*_{V} completes. In addition, constraint projection is also conducted for position correction, and the deformation position of normal nodes in *M*_{S} is determined. Since only *M*_{S} participates in graphic rendering and *M*_{V} does not, this is the final state of hybrid mesh at the current moment.

##### 6.5. Algorithm Flow

*Step 1*. Initialization: compute the global stiffness matrix **K** of *M*_{V}, and initialize the position **p**_{i} and velocity of each node of *M*_{S} and the external force **f**. *Step 2*. In the current ∆*t*, get **f**, compute the deformation of *M*_{V} according to equation (11), obtain the position of nodes on the boundary of *M*_{V}, and update the velocity of the nodes on *M*_{S} according to where = 1/*m*_{i} (*m*_{i} is the node mass of *M*_{S}). *Step 3*. Predict the position of each node on *M*_{S} according to *Step 4*. Project constraints on *M*_{S}: this step has two cases. For control nodes, = . For other nodes on *M*_{S}, solve position correction ∆**p**_{i} on *M*_{S} according to equations (9) and (11) and obtain the according to *Step 5.* According to equation (18), update the **p**_{i} and of each node on *M*_{S}: *Step 6*. Start the next ∆*t* and return to Step 2.

#### 7. Small Incision Processing

In some deformation simulations, operations such as incision on mesh are necessary. Generally, the incision on the surface of a model is small, so the incision displacement of *M*_{V} can be ignored. Moreover, *M*_{V} is not rendered in graphics and thus we focus on the incision of *M*_{S}. To this end, we take the incision treatment as a linear segment which is intersected by *M*_{S} and a sweep surface generated by the virtual cutting tool. When the operation of incision happens, calculate the intersection points on *M*_{S}. There are two types of the intersection points: edge points and facet points. When the operation ends, subdivide the triangular elements and update topology change. Besides, the update of constraint on *M*_{S} is also needed.

As Figure 7 shows, there are three element types for mesh topology change. We use the method of [23] to subdivide the triangular mesh.

**(a)**

**(b)**

**(c)**

To update the constraint functions, the original constraints of the elements in the incision area are deleted. And, the new constraints of the newly subdivided triangular elements are generated. The new constraints generation is introduced in Figure 8.

**(a)**

**(b)**

For the stretch constraint, as Figure 8(a) shows, the intersection P is an edge point. Calculate the barycentric coordinate of *P* (*λ*_{1}, *λ*_{2}, *λ*_{3}) (in Figure 8(a), *λ*_{3} = 0). The initial length of AP *l*_{AP} iswhere *l*_{AB} is initial length of AB in the original constraint. For stretch constraint, as Figure 8(b) shows, the intersection *P* is the facet point. The barycentric coordinate of *P* is (*λ*_{1}, *λ*_{2}, *λ*_{3}). Take the calculation of initial length AP *l*_{AP} as an example. According to law of cosine, angle *α* iswhere *l*_{AB}, *l*_{BC}, and *l*_{AC} are the initial lengths in original constraints AB, BC, and AC, respectively. The initial length of DC *l*_{DC} is

According to the law of cosine, the initial length of AD *l*_{AD} is

The initial length of AP *l*_{AP} is

For the bending constraint, the subdivision of triangle elements does not change the initial angle of two original triangles. Thus, the initial angle of triangular subelements remains unchanged. Since a control point is in the intersection segment, the control point *P* is duplicated due to the element subdivision. These two newly duplicated points are no longer control points.

#### 8. Experiments and Results

##### 8.1. Experimental Environment

We conduct our experiments on a computer using Intel Core i3 (2.7 GHz), 4 GB RAM, HD530 integrated graphics, and Windows 10. The software platform consists of VC++2019, open-source SOFA framework, and OpenGL graphics library.

##### 8.2. Deformation Simulation

In order to present the mesh deformation effect, we implemented our proposed method on a simple deformation simulation.

In simulation, we use a pin-head to apply the external force **f** on the mesh. As Figure 9(a) shows, at the beginning, the mouse cursor controls the pin-head. When the cursor (pin-head) does not touch the mesh, the position of the pin-head is consistent with the position of the cursor. Once the cursor gets close to the node of mesh on the surface, as Figure 9(b) shows, the pin-head snaps to the closest node if the left mouse button is pressed down. In this situation, the position of pin-head is consistent with this node, not the cursor, and it becomes the external force-applied node. Then, as Figure 9(c) shows, move the cursor to generate external force **f**.

**(a)**

**(b)**

**(c)**

**x** is the current position of the cursor. **x**_{0} is the original cursor position when the pin-head snaps to the mesh node. The coefficient *k* can be set by user’s preference. The **f** is provided to the deformation algorithm. If the left mouse button is released, the external force is withdrawn and the position of pin-head is consistent with the cursor again. The cursor is not rendered in graphics. The flow chart is shown in Figure 10.

Figure 11 shows the deformation effect on two shape models (cube and sphere) and two real mesh models (liver and steak) with tension and pressure. The two real mesh models are from turbosquid.com. We set the material properties as follows: Young’s modulus: 5.5 × 10^{3} Pa; Poisson’s ratio: 0.3. By observing the four models, the surface deformation effects are realistic and the graphic rendering is good.

##### 8.3. Surface Deformation Performance

###### 8.3.1. Accuracy Verification

In order to show the accuracy of our proposed method, we compare the deformation result of liver and steak mesh simulated by our method with FEM [10], PBD [14], and MSM [17]. The mesh used in FEM is a fine tetrahedral volume mesh, which is directly generated by *M*_{S} through the Delaunay algorithm in Tetgen. The reason we use fine volume mesh in FEM is that the deformation experimental results will be used as the reference in experiment for its high accuracy. The mesh used in PBD and MSM is the same surface fine mesh model as the *M*_{S} in our hybrid mesh. The information of cube, sphere, liver, and steak mesh models used in different methods is listed in Table 1.

The specific process is as follows: First, we directly applied vertical upward tension and vertical downward pressure on the liver and steak mesh. The magnitudes are 0.5 N, 1.0 N, and 2.0 N for small deformation and 5.0 N, 10.0 N, and 20.0 N for large deformation. Second, we choose 10 points near the force-applied points of the mesh as sample points. Finally, the positions of the sample points after deformation of two shape models (cube and sphere) and two real mesh models (liver and steak) are recorded, and the deformation displacement *d* are calculated by usingwhere represents the position of the sample points *i* after the deformation and represents the initial position of the sample points *i*. The deformation displacement radar diagrams of cube, sphere, liver, and steak models are shown in Figures 12–15 respectively.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

Meanwhile, we calculate the RMSE (root mean square error) of the deformation displacement of our method, PBD, and MSM. Note that the RMSEs are obtained based on the FEM as reference:where *r* represents the RMSE value, *n* represents the number of the deformation sample points, represents the position of the deformation sample points *i* after the deformation of our method, PBD, and MSM, respectively, and represents the position of the deformation sample points *i* after the deformation of FEM as reference. For RMSE calculation, we choose 50 points nearest to the load-applied point (include the load-applied point) as sample points. Table 2 shows RMSE of deformation displacement of our method, PBD, and MSM on liver and steak mesh models. Table 3 shows RMSE of deformation displacement of our method, PBD, and MSM on cube and sphere shapes.

From the diagrams shown in Figures 12 and 13, for both small deformation (0.5 N, 1.0 N, 2.0 N) and large deformation (5.0 N, 10.0 N, 20.0 N), the deformation displacement curve of our method is closer to that of FEM than PBD and MSM. According to the definition of RMSE, the closer the RMSE value is to zero, the closer the deformation displacements is to the value of FEM. The RMSE value of our method is closer to zero than the value of PBD and MSM on both small deformation and large deformation. The above results verify the accuracy of our method.

###### 8.3.2. Time Verification

For time verification, we compared the consumed time of our method with that of PBD, MSM, and FEM in calculating the deformation results. As shown in Tables 4 and 5, the time consumed by our method, PBD, MSM, and FEM under different forces is recorded on two shape models (cube and sphere) and two real mesh models.

From the result shown in Tables 3 and 4, on the one hand, compared with PBD and MSM, our method costs more computation time due to the addition of the time consumed for *M*_{V} deformation, but this time increment is not too much. On the other hand, compared with FEM, our method reduces large computation cost, since we use coarse volume mesh instead of fine volume mesh in FEM for deformation calculation. Although the study in [18] has reduced some computational cost of FEM by optimizing the numerical schemes of dealing with the stretch part of deformation energy, it needs more elements to present the accurate surface deformation effect than ours. Our proposed method needs fewer parameters of volume element and surface element to present the effect. The real-time performance of our proposed method is guaranteed.

###### 8.3.3. Discussion of the Trade-Off between Accuracy and Time Consumption

In accuracy verification, by comparing the deformation displacement and RMSE of our method with the ones of PBD, MSM, and FEM, the results indicates that the deformation of our proposed method is similar to the one of FEM, which we take as reference to evaluate the accuracy. The deformation result of our method is more accurate, due to the constraint on *M*_{S} from the displacement of *M*_{V}. For the computation time, by comparing the consumed time of our methods with that of PBD, MSM, and FEM, the result indicates that our method has lower time consumption than FEM. Although our method has more time consumption than PBD and MSM, the time increment is not too much since computing the deformation for *M*_{V} does not take too much time. Comprehensively considering the deformation result and time consumption, our proposed method achieves the trade-off between the high accuracy and low time consumption.

##### 8.4. Small Incision Effect

We also developed a small incision simulation using our hybrid mesh model. We use a pin-head as a virtual cutting tool. If the pin-head is close to the mesh surface and the left mouse button is pressed down, the incision operation begins. As the pin-head moves, the intersection happens on the mesh surface. When the left mouse button is released, the incision operation is finished. Figure 16 shows the small incision simulation effects on the liver and steak mesh model. For the incision to be visible, we applied a tiny force on a node near the incision on the surface mesh. From Figure 16, we can see that the incisions of these two shape models and two real mesh models are smooth. The incision processing is compatible with our proposed method.

**(a)**

**(b)**

**(c)**

**(d)**

#### 9. Conclusions

In this paper, we propose a deformation method based on the hybrid mesh model which achieves the trade-off between the high accuracy and low time consumption. First, we proposed the construction of the hybrid mesh from fine surface triangular mesh. Its most important step is the simplification process of the fine surface mesh. Second, we combined the FEM for coarse volume mesh deformation and PBD deformation method for fine surface mesh deformation and proposed the algorithm flow based on the hybrid mesh model. Third, we implemented our method on a simple deformation simulation. The experimental result of the deformation algorithm is presented. Compared with the PBD, MSM, and FEM, our method verifies the high accuracy and low time consumption. Finally, we developed an incision simulation which is compatible with our proposed method.

In the future, our proposed method can be improved further in discontinuity processing. The cutting on the hybrid mesh model is a major problem. It involves the subdivision of *M*_{V}, the way to calculate the deformation of the subdivided *M*_{V}, and how the surface node is constrained on the subdivided *M*_{V}. When a large cut happens, the graphic rendering of incision on *M*_{V} is a problem. It is a challenging task to maintain the trade-off between high accuracy and low time consumption in deformation simulation with subdivision of *M*_{S} and *M*_{V}.

#### Data Availability

No data were used to support this study.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the NSFC Joint Fund under Grants GG2090090072, U1332130, and U1713206; 111 Projects under Grant B07033; 973 Project under Grant 2014CB931804; and Key Research and Development Plan of Anhui Province under Grant 1704a0902051.