Abstract

One of the major challenges in mesh-based deformation simulation in computer graphics is to deal with mesh distortion. In this paper, we present a novel mesh-insensitive and softer method for simulating deformable solid bodies under the assumptions of linear elastic mechanics. A face-based 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 by-product 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 mesh-based and mesh-free 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 mesh-based 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 Courant-Friedrichs-Lewy (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 (S-FEM). 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 well-shaped or ill-shaped 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 domain-based 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 S-FEM proposed by Liu et al. [6] is briefly reviewed in Section 3. Our smoothed and corotational elasticity model (CS-FEM 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.

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 in-depth review [912]. 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 well-known 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 vertex-based stiffness warping to the element-based 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, 1618]. 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 higher-order terms of element rotation. Recently, Civit-Flores and Susín proposed a method to handle degenerate elements in isotropic elastic materials [22]. Our method extends the element-based stiffness warping method to the face-based stiffness warping method.

Because mesh quality is an important performance factor in mesh-based 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 mesh-free method (MFM) [27] was developed which took point-based representations for both the simulation volume and the boundary surface of elastically and plastically deforming solids. Later, many useful techniques have been developed in mesh-free methods (summarized in [10]). In computational science, Liu et al. combined the strain smoothing technique [28] used in mesh-free methods [29] into the FEM and formulated a cell/element-based smoothed finite element method (S-FEM) [6]. Their method inherits the mesh distortion insensitivity of mesh-free methods and the accuracy of the FEM. After the theoretical aspects of S-FEM 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 node-based S-FEM (NS-FEM) [31, 32], the alpha-FEM (-FEM) [33], the edge-based S-FEM (ES-FEM) [34, 35], and the face-based S-FEM (FS-FEM) [36]. As the earliest S-FEM model, the cell-based FEM converges in higher rate than the standard FEM in both displacement and energy norms. The NS-FEM 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 NS-FEM for dynamic problems. It proves to be stable and convergent, but its variational consistency depends on how it is formulated [37]. The ES-FEM-T3 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 NS-FEM and the cell-based S-FEM. The ES-FEM was further extended for 3D problems to form the FS-FEM. Similar to ES-FEM, the FS-FEM 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, S-FEM has been applied to a wide range of practical mechanics problems such as fracture mechanics and fatigue behavior [4043], nonlinear material behavior analysis [35, 38, 4448], plates and shells [4952], piezoelectric structures [43, 5355], heat transfer and thermomechanical problems [5659], vibration analysis and acoustics problems [39, 58, 6062], and fluid and structure interaction problems [6366]. We refer the reader to [67, 68] for recent in-depth reviews of S-FEM. As an alternative, Leonetti and Aristodemo proposed a composite mixed finite element model (CM-FEM) 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 FS-FEM because it is stable, accurate, and insensitive to mesh distortion in linear elasticity simulation using tetrahedral meshes. Essentially, our method uses both FS-FEM in linear elasticity modeling and stiffness warping in rotational distortion elimination. Like FS-FEM, it makes use of strain smoothing technique to avoid the problem related to element distortion instead of remeshing. Different from the original FS-FEM, stiffness warping is combined to produce 3D nonlinear deformable model using the corotational linear elasticity model. Moreover, the stiffness warping is converted from element-based stiffness warping to smoothing domain-based stiffness warping.

3. The Smoothed Finite Element Method

In this section, we will present the main idea of S-FEM, which is the foundation of our model. The common concepts and key equations are also listed here.

3.1. The Idea of S-FEM Models

The S-FEM, a numerical and computational method, was proposed by Liu et al. [6] in 2007 and based on FEM and some mesh-free techniques. In the standard FEM, once the displacement is properly assumed, its strain field is available using the strain-displacement 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 S-FEM 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 S-FEM and the FEM follow the same steps. S-FEM models building on different smoothing domains have different features and properties. Next, we will introduce the common concepts and steps in S-FEM models. In Section 4, we will take the FS-FEM as an example to illustrate the smoothing domain building, the strain filed construction, and discrete linear algebraic system establishment of S-FEM.

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 mesh-free 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 S-FEM models apply the smoothing operation on a strain field and get where is the smoothed strain-displacement 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 S-FEM.

For a linear elasticity model, the standard compatible strain-displacement 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 strain-displacement matrix:where and where is the number of nodes in a smoothing domain . Compared to the standard FEM models, the smoothed strain of the S-FEM 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 strain-displacement matrix is an average of the standard strain-displacement matrices over the smoothing domain .

3.4. Smoothed Stiffness Matrix

Then the S-FEM 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 S-FEM 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 well-posed 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 S-FEM 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 S-FEM models, the FS-FEM is both stable and computationally efficient for 3D problems, compared to the NS-FEM, the cell-based S-FEM, and the ES-FEM. So our method is based on the linear FS-FEM (LFS-FEM). Under large rotations, the linear elasticity model produces geometric distortion. Stiffness warping is good at solving this kind of problem. In the following, LFS-FEM is outlined first. Then stiffness warping is abstracted and extended from an element-based method to a smoothing domain-based 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. LFS-FEM

This section formulates the LFS-FEM presented by Nguyen-Thoi et al. [36] for 3D problems using tetrahedral elements. We first give the construction of face-based smoothing domain on a tetrahedral mesh and then show the equations of smoothed strain corresponding stiffness matrix definition.

Smoothing Domain Construction. In the FS-FEM, 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 strain-displacement matrix over a face-based smoothing domain using the compatible strain-displacement 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 strain-displacement 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 LFS-FEM 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 Domain-Based 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 element-based corotational formula for cloth simulation. Müller and Gross [7] and his colleagues extended it to the elastic solid simulation. Element-based corotational formula in practice gives stable simulations. For FS-FEM, the integration domains of stiffness and forces are based on smoothing domains; the element-based corotational formula cannot be used directly. So we extend the element-based corotational formula to smoothing domain-based corotational formula.

For completeness, the element-based 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 volume-weighted average operation is defined on the rotational matrices of elements associated with a smoothing domain. In LFS-FEM, 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 : element-based rotational matrix projection, matrix to quaternion transformation, quaternion interpolation, and quaternion to matrix transformation. Because the element-based 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 element-based 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 (L-FEM) 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).

(1) initialize ,
(2) forall element compute the standard strain-
   displacement matrix
(3) build smoothing domains
(4) forall smoothing domain
(5)  compute using Eq. (15)
(6)  compute using Eq. (17)
(7) endfor
(8)
(9) loop
(10) forall element compute using Eq. (20)
(11) forall smoothing domain compute
(12) assemble matrix
(13) assemble vector
(14) apply external force
(15) 
(16) 
(17) solve from
(18) 
(19) 
(20) endloop

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 CS-FEM is more than that of C-FEM. The execution time of both methods is log-log linear dependent on the number of elements. Compared to C-FEM, CS-FEM spends extra time on smoothing domain construction and smoothing domain-based data processing. The smoothing domain-based data include smoothing domain-based 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, CS-FEM only takes extra time computing .

We summarized the simulation time distribution in each step in Table 1. The table shows that CS-FEM spends less than 2% extra time () computing smoothing domain-based rotational matrix . The bottleneck of CS-FEM is linear algebraic equations assembly () and solving (). If the mesh is less than 6 K tetrahedron, CS-FEM reaches 60 fps and can be used in real-time applications.

Mesh Distortion Sensitivity under Pressure. This example was designed to measure the mesh distortion sensitivity of CS-FEM 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 10-node 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 CS-FEM using 4-node tetrahedral element is less accurate than that of CM-FEM using 10 displacement nodes and 4 stress regions () [71], but the displacement of CS-FEM is more accurate than that of . The results of CS-FEM are less accurate than the nonlinear FS-FEM (NLFS-FEM) but more accurate than those of the linear FEM (L-FEM). It is also shown that CS-FEM is less sensitive to mesh distortion than L-FEM in both displacement and strain energy.

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 FEM-H20 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 CS-FEM is rendered in Figure 7. The transient vertical displacements of node A were plotted in Figure 5.

The simulation frequencies of L-FEM and C-FEM begin to deviate from those of the regular mesh () if (lines in purple and green) and were as low as one-third of those of the regular mesh if (line in green), while the simulation frequency of LFS-FEM, NLFS-FEM, and CS-FEM 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 CS-FEM is less accurate than that of NLFS-FEM, while it is more accurate than that of LFS-FEM, L-FEM, and C-FEM, compared to that of the reference solution. Because the relative errors (refer to (%) in Table 2) of CS-FEM are always less than those of C-FEM and L-FEM, this proves that our method is less sensitive to mesh distortion than the methods without strain smoothing (C-FEM and L-FEM) under geometrically nonlinear deformations.

Geometric Distortion under Rotational Deformations. This example was designed to examine the volume gains of CS-FEM 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 CS-FEM and C-FEM are approaching zero. Without the help of corotational operation, the maximum total volume gains of L-FEM and LFS-FEM are both larger than 1.0. It is shown that, with the help of corotational operation, CS-FEM alleviates geometric distortion under rotational deformations.

Softening Effect. The goal of this test was to demonstrate the softening effect of our method compared to the C-FEM 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 CS-FEM (in red) is not less than those obtained from C-FEM using the same initial mesh (in gray). In Table 2 and Figure 8, the vertical displacement of CS-FEM is larger than that of C-FEM. The strain energy obtained from CS-FEM also is not less than that of the C-FEM as shown in Figure 11. The above results show the softening effect of CS-FEM.

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 element-based stiffness warping to smoothing domain-based 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 C-FEM without adding significant computational burden. Previous methods such as FEM and C-FEM 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 domain-based stiffness warping approach, the geometric distortion under large rotation is eliminated. The results have also shown that our method is softer than the C-FEM 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).