Abstract
One of the major challenges in meshbased deformation simulation in computer graphics is to deal with mesh distortion. In this paper, we present a novel meshinsensitive and softer method for simulating deformable solid bodies under the assumptions of linear elastic mechanics. A facebased strain smoothing method is adopted to alleviate mesh distortion instead of the traditional spatial adaptive smoothing method. Then, we propose a way to combine the strain smoothing method and the corotational method. With this approach, the amplitude and frequency of transient displacements are slightly affected by the distorted mesh. Realistic simulation results are generated under large rotation using a linear elasticity model without adding significant complexity or computational cost to the standard corotational FEM. Meanwhile, softening effect is a byproduct of our method.
1. Introduction
Physically based dynamic simulation of deformation has been an active research topic in computer graphics community for more than 30 years. Since Terzopoulos and his colleagues [1] introduced this field into graphics in 1980s, a large body of works have been published in pursuit of visually and physically realistic animation of deformable objects. According to the method for discretization of partial differential equations (PDEs) governing dynamic elasticity deformation, physically based methods are divided into meshbased and meshfree methods. Our method falls mainly within the first category by using the finite element method (FEM).
The FEM is one of the most widely used meshbased methods in solving PDEs. It discretizes the objects into a mesh of finite connected nodes and approximates the field within an element by interpolation of the values at the nodes of the element. As a result, the quality of the mesh plays an important role in the FEMs. Adopting either explicit integration or implicit integration methods, badly shaped elements directly lower down the accuracy and speed of numerical solutions of governing PDEs [2]. To make it worse, the mesh is not static and is changing with the object deformation. The traditional and most promising way to improve mesh quality is the adaptive remeshing. However, remeshing methods are daunting because they involve tedious mesh topology and geometry operations and projections of field variables from the previous mesh [3, 4]. Even excellent field variable projection schemes might lead to significant errors in displacements and velocities [5]. Remeshing methods include refinement and coarsening. Both of them modify the edge length of the mesh, which disturbs the CourantFriedrichsLewy (CFL) condition and introduces stability problems in turn. Our approach departs from this traditional viewpoint by smoothing the strain field on the mesh instead of smoothing the mesh. This idea is adopted by Liu and his colleagues [6] in their smoothed finite element methods (SFEM). This new point of view avoids the problems caused by remeshing such as instability. To make the simulation fast, we combine the strain smoothing technique with the stiffness warping approach [7]. Our results show that this method is capable of producing correct simulation results using either wellshaped or illshaped meshes under large rotations.
In summary, our main contributions are as follows:(i)A novel smoothed pseudolinear elasticity model is presented for deformable solid simulation. It adopts a linear elasticity model for nonlinear simulation and compensates for the error of linear calculation in the nonlinear simulation using stiffness warping. Different from the previous approaches, the pseudolinear elasticity model is built on the smoothed finite element method. It is the first time that the smoothed finite element method has been used for deformable solid simulation in computer graphics.(ii)A novel smoothing domainbased stiffness warping approach is proposed to accommodate the change of integration domains in our smoothed pseudolinear elasticity model.(iii)The strain smoothing technique adopted provides an alternative way to avoid problems related to mesh distortion met in deformable solid simulation in computer graphics.
The paper is structured as follows. First, closely related researches are presented in Section 2. Then the SFEM proposed by Liu et al. [6] is briefly reviewed in Section 3. Our smoothed and corotational elasticity model (CSFEM for short) is presented in Section 4. Next, the experiment results and analysis are delivered in Section 5. Finally, conclusions and future studies are sketched in Section 6.
2. Related Work
Several researchers have developed physically based models for deformable solid simulation since Terzopoulos and his colleagues introduced methods for simulating elastic [1] and inelastic [8] materials into the graphics community. We refer the reader to the survey papers that focus on deformation modeling in computer graphics for indepth review [9–12]. Here, we only focus on works on the FEM.
The FEM, one of the most popular methods, has been widely used in both elastic material (summarized in [10] and later reviews) and inelastic material simulations [13, 14]. Linear elasticity FEM models are applicable to interactive application because of their stability and efficiency. However, they are not suitable for large rotational deformations because of their wellknown geometric distortion. Capell et al. proposed dividing an object into small parts based on its skeleton manually [15]. Müller and his colleagues took a different approach: stiffness warping [16]. Then, they further improved the vertexbased stiffness warping to the elementbased stiffness warping and extended their approach to inelastic material simulation [7]. Their method produces fast and robust simulation results and has been widely used in interactive simulation [7, 14, 16–18]. Later, Chao et al. [19], McAdams et al. [20], and Barbic [21] gave an exact corotational FEM stiffness matrix for a linear tetrahedral element by adding the higherorder terms of element rotation. Recently, CivitFlores and Susín proposed a method to handle degenerate elements in isotropic elastic materials [22]. Our method extends the elementbased stiffness warping method to the facebased stiffness warping method.
Because mesh quality is an important performance factor in meshbased simulations, there is a huge body of literature on spatial or geometric adaptive methods. Remeshing is a widely used approach in maintaining mesh quality. Basis function refinement [23, 24], mesh embedding [25, 26], and other mixed models also work well in general. Manteaux and his colleagues gave a thorough review on them [4] recently in computer graphics. To avoid element distortion encountered in the FEM, a meshfree method (MFM) [27] was developed which took pointbased representations for both the simulation volume and the boundary surface of elastically and plastically deforming solids. Later, many useful techniques have been developed in meshfree methods (summarized in [10]). In computational science, Liu et al. combined the strain smoothing technique [28] used in meshfree methods [29] into the FEM and formulated a cell/elementbased smoothed finite element method (SFEM) [6]. Their method inherits the mesh distortion insensitivity of meshfree methods and the accuracy of the FEM. After the theoretical aspects of SFEM were clarified and its properties were confirmed by numerical experiments [30], the concept of smoothing was extended to formulate a series of smoothed FEM models such as the nodebased SFEM (NSFEM) [31, 32], the alphaFEM (FEM) [33], the edgebased SFEM (ESFEM) [34, 35], and the facebased SFEM (FSFEM) [36]. As the earliest SFEM model, the cellbased FEM converges in higher rate than the standard FEM in both displacement and energy norms. The NSFEM can alleviate the volumetric locking problem effectively. But it is usually less computationally efficient and temporally unstable. The FEM was proposed by Liu et al. to avoid the spurious nonzero energy modes in NSFEM for dynamic problems. It proves to be stable and convergent, but its variational consistency depends on how it is formulated [37]. The ESFEMT3 often offers superconvergent and better solution than the standard FEM. It also excels in handling the volumetric and/or bending locking problem using bubble functions [38]. It is immune from the “overly soft” problem met in NSFEM and the cellbased SFEM. The ESFEM was further extended for 3D problems to form the FSFEM. Similar to ESFEM, the FSFEM is more accurate than the standard FEM using the same T4 mesh for dynamic problems [39] and both linear and nonlinear problems [36]. Because of its excellent properties, SFEM has been applied to a wide range of practical mechanics problems such as fracture mechanics and fatigue behavior [40–43], nonlinear material behavior analysis [35, 38, 44–48], plates and shells [49–52], piezoelectric structures [43, 53–55], heat transfer and thermomechanical problems [56–59], vibration analysis and acoustics problems [39, 58, 60–62], and fluid and structure interaction problems [63–66]. We refer the reader to [67, 68] for recent indepth reviews of SFEM. As an alternative, Leonetti and Aristodemo proposed a composite mixed finite element model (CMFEM) to generate new smoothed operators [69]. A composite triangular mesh is assumed over the domain. An element is subdivided into three triangular regions and linear stress/quadratic displacement interpolations are used to approximate the exact stress/displacement fields. Their method is also insensitive to mesh distortion and applicable to the elastic and plastic analysis in plane problems [70]. Bilotta and his colleagues adopted the same idea and proposed a composite mixed finite element model for 3D structural problems with small strain fields [71, 72].
Our method is based on FSFEM because it is stable, accurate, and insensitive to mesh distortion in linear elasticity simulation using tetrahedral meshes. Essentially, our method uses both FSFEM in linear elasticity modeling and stiffness warping in rotational distortion elimination. Like FSFEM, it makes use of strain smoothing technique to avoid the problem related to element distortion instead of remeshing. Different from the original FSFEM, stiffness warping is combined to produce 3D nonlinear deformable model using the corotational linear elasticity model. Moreover, the stiffness warping is converted from elementbased stiffness warping to smoothing domainbased stiffness warping.
3. The Smoothed Finite Element Method
In this section, we will present the main idea of SFEM, which is the foundation of our model. The common concepts and key equations are also listed here.
3.1. The Idea of SFEM Models
The SFEM, a numerical and computational method, was proposed by Liu et al. [6] in 2007 and based on FEM and some meshfree techniques. In the standard FEM, once the displacement is properly assumed, its strain field is available using the straindisplacement relation, which is called fully compatible strain field. Then the standard FEM is formulated using the standard Galerkin formulation. The fully compatible FEM model leads to locking behavior for many problems. The assumed piecewise continuous displacement field induces a discontinuous strain field on all the element interfaces and inaccuracy in stress solutions in turn. In addition, the Jacobian matrix related to domain mapping becomes badly conditioned on distorted elements, leading to deterioration in solution accuracy. It also prefers quadrilateral elements and hexahedral elements and gives poor accuracy, especially for stresses, for easily obtained triangular and tetrahedral elements.
To avoid the upper problems of the FEM, the SFEM uses the smoothing technique to modify the fully compatible strain field or construct a strain field without computing the compatible strain field over all the smoothing domains. Then the smoothed Galerkin weak form, instead of the Galerkin weak form used in the standard FEM, is used to establish the discrete linear algebraic system of equations. Except the strain field construction and discrete linear algebraic system establishment, both the SFEM and the FEM follow the same steps. SFEM models building on different smoothing domains have different features and properties. Next, we will introduce the common concepts and steps in SFEM models. In Section 4, we will take the FSFEM as an example to illustrate the smoothing domain building, the strain filed construction, and discrete linear algebraic system establishment of SFEM.
3.2. Smoothing Operation
We first introduce the integral representation of the approximation of function :where is a prescribed smoothing function defined in the smoothing domain of point . Smoothing domain can be moved and overlapped for different . must follow the conditions of the partition of unity, positivity, and decay. Similarly, the integral form of the gradient of function can be represented asNote that (1) and (2) are the standard forms of smoothing operation and were widely used in the smoothed particle hydrodynamics and meshfree methods for field approximation. Applying Green’s divergence theorem on (2), we getwhere is the boundary of and is the outward normal on . Equation (3) holds if is continuous and at least piecewise differentiable.
For simplicity, a local constant smoothing function is used.. Substituting (4) into (3), we get the smoothed gradientThe integrals involving gradient of a function are recast into the boundary integrals involving only the function and the boundary normals by the gradient smoothing technique.
3.3. Strain Smoothing
The SFEM models apply the smoothing operation on a strain field and get where is the smoothed straindisplacement matrix of smoothing domain , is the nodal displacement matrix, and is the number of nodes in a smoothing domain. The smoothed strain is assumed to be constant in each smoothed domain in SFEM.
For a linear elasticity model, the standard compatible straindisplacement matrix iswhere is the symmetric gradient operator, is the finite element shape function matrix, and is the number of nodes in an element. Note that, applying smoothing operation on , we get the smoothed straindisplacement matrix:where and where is the number of nodes in a smoothing domain . Compared to the standard FEM models, the smoothed strain of the SFEM models only depends on the nodal displacements, the normal of the element boundary, and the shape function matrix. Usually one Gaussian point is used for line integration along each segment of boundary . From (8), we can see that the smoothed straindisplacement matrix is an average of the standard straindisplacement matrices over the smoothing domain .
3.4. Smoothed Stiffness Matrix
Then the SFEM models use the smoothed strain to compute the potential energy and define the smoothed Galerkin weak form:where is the material constant matrix, is a trial function, is a test function, is the external body force applied over the problem domain, and is the external force applied on the natural boundary.
The SFEM models formulated using (10) are variationally consistent and spatially stable if the number of linear independent smoothing domains is sufficient and converge to the exact solution of a physically wellposed linear elasticity problem. Substituting the assumed approximations and ,into (10), we have the standard discretized algebraic system of equations:where is the nodal displacement, is the nodal displacements for all the nodes in the SFEM model, and is the load vector:and is the global smoothed stiffness matrix defined by is a sparse matrix with nonzero entries if node and node share the same smoothing domain. For a static analysis problem, the nodal displacements can be obtained by solving (12).
4. Smoothed Pseudolinear Elasticity Model
In this section, our smoothed corotational elasticity model is formulated. In SFEM models, the FSFEM is both stable and computationally efficient for 3D problems, compared to the NSFEM, the cellbased SFEM, and the ESFEM. So our method is based on the linear FSFEM (LFSFEM). Under large rotations, the linear elasticity model produces geometric distortion. Stiffness warping is good at solving this kind of problem. In the following, LFSFEM is outlined first. Then stiffness warping is abstracted and extended from an elementbased method to a smoothing domainbased method. Next, the dynamics of the system are introduced to simulate the dynamic behavior of an object. In the end, the procedure of simulation is summed up.
4.1. LFSFEM
This section formulates the LFSFEM presented by NguyenThoi et al. [36] for 3D problems using tetrahedral elements. We first give the construction of facebased smoothing domain on a tetrahedral mesh and then show the equations of smoothed strain corresponding stiffness matrix definition.
Smoothing Domain Construction. In the FSFEM, both the strain smoothing and the stiffness matrix integration are based on local smoothing domains. These local smoothing domains are constructed based on faces of the neighboring elements such that and , , in which is the total number of smoothing domains (i.e., faces here) in . For tetrahedral elements, a smoothing domain is created by connecting three nodes of its associated face to the central point of two adjacent tetrahedral elements as shown in Figure 1.
The definition of the smoothed strain and stiffness matrix can be induced directly when the smoothing domains are constructed. From (8), we get the smoothed straindisplacement matrix over a facebased smoothing domain using the compatible straindisplacement matrix:where is the volume of smoothing domain , is the number of elements, is the volume of th element in smoothing domain , and , a matrix, is the standard compatible straindisplacement matrix of element in smoothing domain . is a matrix. for boundary faces and for inner faces.
The smoothed strain is written aswhere is the collection of the displacements of nodes in smoothing domain .
Then, the local smoothed stiffness matrix is given as follows:
The entries of are defined byThe local smoothed stiffness matrices are assembled to produce the entry of the global smoothed stiffness matrix :
We note that LFSFEM extends the standard FEM by applying gradient smoothing technique beyond elements. Linear elastic model only applies to small deformation and leads to visual artifacts under large rotational deformations. Next, we introduce stiffness warping approach to handle the elastic body suffering small deformation with large rotation.
4.2. Smoothing DomainBased Stiffness Warping
The stiffness warping method was proposed by Müller et al. [16] to remove the artifacts of growth in volume that linear elastic forces show while keeping the governing equation linear. The key to stiffness warping is how to extract the rotational components from deformations. The formula used to extract the rotational matrix is known as corotational formula. For stiff materials with little deformation but arbitrary rigid body motion, keeping track of rotations of a global rigid body frame would yield acceptable results as Terzopoulos and Witkin did [8], which still yields the typical artifacts of a linear model under large deformations other than rigid body modes. For nonstiff materials with large deformations, keeping track of individual rotations of every vertex alleviates the artifacts as Müller and his colleagues did [16], which brings possible ghost forces from the nonzero total elastic forces. Later, Etzmuß and his colleagues [17] proposed an elementbased corotational formula for cloth simulation. Müller and Gross [7] and his colleagues extended it to the elastic solid simulation. Elementbased corotational formula in practice gives stable simulations. For FSFEM, the integration domains of stiffness and forces are based on smoothing domains; the elementbased corotational formula cannot be used directly. So we extend the elementbased corotational formula to smoothing domainbased corotational formula.
For completeness, the elementbased corotational formula [7] is given below:where is an operator extracting rotational matrix from the deformation gradient matrix . , , is the deformed position of node in element , and is the undeformed position of node in element . There are three ways to define operator in computer graphics: QR factorization [73], singular value decomposition (SVD) [74], and polar decomposition [7]. Because polar decomposition is fast and stable [22], it is employed here.
To get the rotational matrix of the smoothing domain , a volumeweighted average operation is defined on the rotational matrices of elements associated with a smoothing domain. In LFSFEM, a smoothing domain is associated with elements. For smoothing domains associated with boundary faces, equals the rotational matrix of their unique associated element. For smoothing domains associated with inner faces, it takes four steps to get : elementbased rotational matrix projection, matrix to quaternion transformation, quaternion interpolation, and quaternion to matrix transformation. Because the elementbased rotational matrices are based on their own undeformed shape, a project operation should be applied to put them under the same coordinate system before interpolation. Assume that the local coordinate system of element in the smoothing domain is . is , is , and is . is the normalized . The projection matrix is defined as , where is the selected reference local coordinate system. Then the elementbased rotational matrix is projected by There is no meaningful interpolation formula which directly applies to rotation matrices. An alternative way is using quaternion interpolation. The quaternions are quite amenable to interpolation. The linear quaternion interpolation (i.e., lerp) yields a secant between the two quaternions, while spherical linear interpolation (i.e., slerp), performing the shortest great arc interpolation, gives the optimal interpolation curve between two rotations. So the projected matrix is first transformed to quaternion .
Then slerp is applied to get the interpolated quaternion :where the interpolation coefficient is .
In the fourth step, is transformed back to a matrix .
Applying the rotational matrix on each smoothing domain, we get the elastic forces acting on the nodes by a smoothing domain:where and are the deformed and undeformed nodal positions, respectively. is a matrix that contains copies of the matrix along its diagonal. The force offset vector is invariant and can be precomputed to accelerate the algorithm. By (23), the rotational part of deformation is cancelled by . The elastic forces are computed in the local coordinate of smoothing domain and then rotated back to the global coordinate by .
The global elastic forces acting on a node are obtained by summing the elastic forces (see (23)) from the node’s adjacent smoothing domains. Then the elastic forces of the entire mesh are reached bywhere the global stiffness matrix is the summation of the smoothing domain’s rotated stiffness matrix and the global force offset vector is the summation of the smoothing domain’s force offset vector .
4.3. Dynamics
The following governing equation describes the dynamics of the system:where is the current object state; and are the first and second derivatives of with respect to time. is the total number of nodes in the entire mesh. is the mass matrix and is the damping matrix. The lumped mass matrix is used for efficiency here. Rayleigh damping is used to compute the damping matrix , where and are known as the mass damping and stiffness damping coefficients, respectively. is the external loads vector. We omit the time parameter of using instead to lighten the notation.
Equation (25) is discretized in time domain and solved using implicit Euler. Substituting (24), , , and into (25), we get a group of linear algebraic equations:where is the velocity vector, is the time step, and stands for the th time step.
4.4. Dynamic Simulation Algorithm
The whole dynamic simulation process is summed up in the following pseudocode in Algorithm 1. It can be seen that the procedure of our algorithm is in great similarity to that of conventional linear FEM (LFEM) except the computation of smoothed strain matrix in line (5), the computation of smoothed stiffness matrix in line (6), and the computation of offset force in line (13). We use implicit Euler to solve the dynamic systems using lines (15–17), which makes our algorithm unconditionally stable. The traditional conjugate gradient solver is used for linear algebraic equation solving in line (17).

5. Results and Discussion
All of the experiments were executed on a machine with 2.3 GHz dual core CPU and 6 GB of memory.
Computing Time. To illustrate the computing time, we simulated a beam fixed in one end and deformed under gravity with Poisson’s ratio and elasticity modulus Pa. It took 50 iterations to solve the linear equations (see (26)). The CPU time of the overall simulation was measured and plotted as a function of elements number in the simulation mesh in Figure 2. It is shown that the execution time of CSFEM is more than that of CFEM. The execution time of both methods is loglog linear dependent on the number of elements. Compared to CFEM, CSFEM spends extra time on smoothing domain construction and smoothing domainbased data processing. The smoothing domainbased data include smoothing domainbased stiffness matrix and rotational matrix . Because the local stiffness matrix is constant, it is computed at the preprocessing step and reused in the following steps. During each iteration, CSFEM only takes extra time computing .
We summarized the simulation time distribution in each step in Table 1. The table shows that CSFEM spends less than 2% extra time () computing smoothing domainbased rotational matrix . The bottleneck of CSFEM is linear algebraic equations assembly () and solving (). If the mesh is less than 6 K tetrahedron, CSFEM reaches 60 fps and can be used in realtime applications.
Mesh Distortion Sensitivity under Pressure. This example was designed to measure the mesh distortion sensitivity of CSFEM under small strains.
The results were measured by energy error norm and displacement error norm, respectively [36]. The energy error norm was defined as :where is the numerical solution of strain energy and is the exact solution of strain energy. The displacement error norm of node was defined bywhere is the numerical solution of displacement and is the exact solution of displacement.
In this experiment, a 3D cubic cantilever subjected to a uniform pressure on its upper face was considered. Its size is () and it is discretized using a mesh including 216 nodes and 625 tetrahedral elements (refer to Figure 3). The related parameters were taken as Poisson’s ratio , elasticity modulus , and pressure . The exact solution of the problem is unknown; we used the reference solution provided by [36]. The reference solution therein was found by using standard FEM with a mesh including 30,204 nodes and 20,675 10node tetrahedron elements. The reference solution of the strain energy is and the deflection at point B is 3.3912.
The experiment was taken using five distorted meshes. Those meshes were generated using an irregular factor between and , random numbers , , and between and , and the edge lengths , , and of a grid in a cubic mesh in the following fashion: The meshes were generated using distortion coefficient from 0.0 to 0.4.
Figures 4(a) and 4(b) show the energy error norm and displacement error norm at point B. It is shown that the strain energy of CSFEM using 4node tetrahedral element is less accurate than that of CMFEM using 10 displacement nodes and 4 stress regions () [71], but the displacement of CSFEM is more accurate than that of . The results of CSFEM are less accurate than the nonlinear FSFEM (NLFSFEM) but more accurate than those of the linear FEM (LFEM). It is also shown that CSFEM is less sensitive to mesh distortion than LFEM in both displacement and strain energy.
(a)
(b)
Mesh Distortion Sensitivity under Large Rotation. The aim of this experiment was to examine the mesh sensitivity of our method under large rotational deformations. A 3D cantilever beam subjected to gravity was considered. The size of the beam was (m × 0.3 m × 0.3 m) and it is discretized using a mesh including 160 nodes and 405 tetrahedral elements. The related parameters were taken as Pa and . The exact solution of the problem is unknown. The reference solution was found by using the nonlinear FEMH20 with a fine mesh including 45,376 nodes and 10,125 elements. The reference solution of vertical displacement at node A at time 0.25 s is 16.8977 cm. The experiment was taken using five distorted meshes with 0.0 to 0.4. The mesh plotted in Figure 6 is the distorted mesh generated with .
The configuration at = 0.25 s generated by CSFEM is rendered in Figure 7. The transient vertical displacements of node A were plotted in Figure 5.
(a) LFEM
(b) CFEM
(c) LFSFEM
(d) NLFSFEM
(e) CSFEM
The simulation frequencies of LFEM and CFEM begin to deviate from those of the regular mesh () if (lines in purple and green) and were as low as onethird of those of the regular mesh if (line in green), while the simulation frequency of LFSFEM, NLFSFEM, and CSFEM does not shift from that of the regular mesh until . The vertical displacement of A at time 0.25 s was listed in Table 2. It shows that the solution of CSFEM is less accurate than that of NLFSFEM, while it is more accurate than that of LFSFEM, LFEM, and CFEM, compared to that of the reference solution. Because the relative errors (refer to (%) in Table 2) of CSFEM are always less than those of CFEM and LFEM, this proves that our method is less sensitive to mesh distortion than the methods without strain smoothing (CFEM and LFEM) under geometrically nonlinear deformations.
Geometric Distortion under Rotational Deformations. This example was designed to examine the volume gains of CSFEM under large rotational deformations. A 3D cantilever beam subjected to gravity was considered. The size of the beam was (2.0 m × 0.5 m × 0.5 m) and it is discretized using a mesh including 475 nodes and 1440 tetrahedral elements with and Pa. The simulation results were rendered in Figure 8. Figure 9 shows that the total volume gains of CSFEM and CFEM are approaching zero. Without the help of corotational operation, the maximum total volume gains of LFEM and LFSFEM are both larger than 1.0. It is shown that, with the help of corotational operation, CSFEM alleviates geometric distortion under rotational deformations.
(a)
(b)
(c)
(a) CFEM
(b) CSFEM
(c) LFEM
(d) LFSFEM
Softening Effect. The goal of this test was to demonstrate the softening effect of our method compared to the CFEM model. A bridge was fixed at the bottom and subjected to gravity and a static load ( N). The bridge was discretized using a mesh including 3,923 nodes and 8,680 tetrahedral elements with , Pa, and mass damping coefficient = 1.0. The simulation results were rendered in Figure 10. It is shown that the deformation obtained from CSFEM (in red) is not less than those obtained from CFEM using the same initial mesh (in gray). In Table 2 and Figure 8, the vertical displacement of CSFEM is larger than that of CFEM. The strain energy obtained from CSFEM also is not less than that of the CFEM as shown in Figure 11. The above results show the softening effect of CSFEM.
6. Conclusion
Our paper has provided a novel method for simulating elastic solid bodies. The key to our technique is a smoothed pseudolinear elasticity model utilizing the stiffness warping approach, which has been extended from elementbased stiffness warping to smoothing domainbased stiffness warping to accommodate the integration domain remodeling. To our knowledge, it is the first time to apply the smoothed finite element method in deformable solid body simulation in computer graphics and also the first time to combine the smoothed finite element method with the stiffness warping method. Our method achieves the same results as the CFEM without adding significant computational burden. Previous methods such as FEM and CFEM are sensitive to mesh distortion and produce shifted displacements during deformation. We have experimentally verified that our method minimizes the impact of mesh distortion on deformation simulation. Using the smoothing domainbased stiffness warping approach, the geometric distortion under large rotation is eliminated. The results have also shown that our method is softer than the CFEM in terms of displacements and strain energy. In the future, we plan to combine the exact rotational matrix and degenerate element handling technique to produce realistic simulation even under extreme stretch and inverted shapes.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Special thanks should go to Doctor Antonio Bilotta, Assistant Professor of DIMES, University of Calabria, for providing the results of his study. The authors also would like to thank the authors of the original studies included in this analysis. This work was partially supported by the National Natural Science Foundation of China (Grant no. 61170203) and Beijing Natural Science Foundation (4174094).