# A Smoothed Finite Element-Based Elasticity Model for Soft Bodies

Authors: Juan Zhang; Mingquan Zhou; Youliang Huang; Pu Ren; Zhongke Wu; Xuesong Wang; Shi Feng Zhao
Journal: Mathematical Problems in Engineering
 (2017)

---

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


---

## Body

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


## 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 in-depth 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 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, 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 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 [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 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 Section4, 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 functionw∈Ω:(1)whx=∫ΩxwξΦx-ξdΩ,where Φ(x) is a prescribed smoothing function defined in the smoothing domain Ωx∈Ω of point x. Smoothing domain Ωx can be moved and overlapped for different x. Φx must follow the conditions of the partition of unity, positivity, and decay. Similarly, the integral form of the gradient of function wh can be represented as(2)∇whx=∫Ωx∇wξΦx-ξdΩ.Note 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 get(3)∇whx=∫ΓxwξnξΦx-ξdΓ-∫Ωxwξ∇Φx-ξdΩ,where Γx is the boundary of Ωx and n(ξ) is the outward normal on Γx. Equation (3) holds if w is continuous and at least piecewise differentiable.For simplicity, a local constant smoothing function is used.(4)Φξ=1Vxξ∈Ωx0ξ∉Ωx.Vx=∫Ωx‍dξ. Substituting (4) into (3), we get the smoothed gradient(5)∇¯whx=1Vx∫ΓxwξnξdΓ.The integrals involving gradient of a function w 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(6)ε¯kh=B¯kdk,where B¯k=[B¯k,1⋯B¯k,N(n,k)] is the smoothed strain-displacement matrix of smoothing domain k, dk=[dk,1⋯dk,Nn,k]T is the nodal displacement matrix, and N(n,k) is the number of nodes in a smoothing domain. The smoothed strain ε¯kh is assumed to be constant in each smoothed domain in S-FEM.For a linear elasticity model, the standard compatible strain-displacement matrix is(7)Bx=∇sNx,where ∇s is the symmetric gradient operator, N=[N1⋯NN(n,e)] is the finite element shape function matrix, and N(n,e) is the number of nodes in an element. Note that, applying smoothing operation on B, we get the smoothed strain-displacement matrix:(8)B¯kx=∫ΩkBξΦx-ξdΩ=1Vk∫ΓkNξnξdΓ,where Vk=∫Ωk‍dΩ and (9)B¯k,ix=1Vk∫ΓkNiξnξdΓ,i∈1⋯Nn,k,where N(n,k) is the number of nodes in a smoothing domain k. 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 Ωk.


### 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:(10)∫Ωδε¯TDε¯dΩ-∫ΩδuTbdΩ-∫ΓδuTtdΓ=0,where D is the material constant matrix, u is a trial function, δu is a test function, b is the external body force applied over the problem domain, and t 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 uh and δuh,(11)uh=∑i=1Nn,eNixdi,δuh=∑i=1Nn,eNixδdi,into (10), we have the standard discretized algebraic system of equations:(12)K¯d=f,where di is the nodal displacement, d is the nodal displacements for all the nodes in the S-FEM model, and f is the load vector:(13)f=∫ΩNiTxbdΩ+∫ΓNiTxtdΓand K¯ is the global smoothed stiffness matrix defined by(14)K¯ij=∫ΩB¯iTDB¯jdΩ.K¯ is a sparse matrix with nonzero entries K¯ij if node i and node j share the same smoothing domain. For a static analysis problem, the nodal displacements can be obtained by solving (12).


## 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 Section4, 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 functionw∈Ω:(1)whx=∫ΩxwξΦx-ξdΩ,where Φ(x) is a prescribed smoothing function defined in the smoothing domain Ωx∈Ω of point x. Smoothing domain Ωx can be moved and overlapped for different x. Φx must follow the conditions of the partition of unity, positivity, and decay. Similarly, the integral form of the gradient of function wh can be represented as(2)∇whx=∫Ωx∇wξΦx-ξdΩ.Note 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 get(3)∇whx=∫ΓxwξnξΦx-ξdΓ-∫Ωxwξ∇Φx-ξdΩ,where Γx is the boundary of Ωx and n(ξ) is the outward normal on Γx. Equation (3) holds if w is continuous and at least piecewise differentiable.For simplicity, a local constant smoothing function is used.(4)Φξ=1Vxξ∈Ωx0ξ∉Ωx.Vx=∫Ωx‍dξ. Substituting (4) into (3), we get the smoothed gradient(5)∇¯whx=1Vx∫ΓxwξnξdΓ.The integrals involving gradient of a function w 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(6)ε¯kh=B¯kdk,where B¯k=[B¯k,1⋯B¯k,N(n,k)] is the smoothed strain-displacement matrix of smoothing domain k, dk=[dk,1⋯dk,Nn,k]T is the nodal displacement matrix, and N(n,k) is the number of nodes in a smoothing domain. The smoothed strain ε¯kh is assumed to be constant in each smoothed domain in S-FEM.For a linear elasticity model, the standard compatible strain-displacement matrix is(7)Bx=∇sNx,where ∇s is the symmetric gradient operator, N=[N1⋯NN(n,e)] is the finite element shape function matrix, and N(n,e) is the number of nodes in an element. Note that, applying smoothing operation on B, we get the smoothed strain-displacement matrix:(8)B¯kx=∫ΩkBξΦx-ξdΩ=1Vk∫ΓkNξnξdΓ,where Vk=∫Ωk‍dΩ and (9)B¯k,ix=1Vk∫ΓkNiξnξdΓ,i∈1⋯Nn,k,where N(n,k) is the number of nodes in a smoothing domain k. 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 Ωk.


## 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:(10)∫Ωδε¯TDε¯dΩ-∫ΩδuTbdΩ-∫ΓδuTtdΓ=0,where D is the material constant matrix, u is a trial function, δu is a test function, b is the external body force applied over the problem domain, and t 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 uh and δuh,(11)uh=∑i=1Nn,eNixdi,δuh=∑i=1Nn,eNixδdi,into (10), we have the standard discretized algebraic system of equations:(12)K¯d=f,where di is the nodal displacement, d is the nodal displacements for all the nodes in the S-FEM model, and f is the load vector:(13)f=∫ΩNiTxbdΩ+∫ΓNiTxtdΓand K¯ is the global smoothed stiffness matrix defined by(14)K¯ij=∫ΩB¯iTDB¯jdΩ.K¯ is a sparse matrix with nonzero entries K¯ij if node i and node j 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 Ω=∑k=1Nf‍Ωk and Ωp∩Ωq≠∅, p≠q, in which Nf is the total number of smoothing domains (i.e., faces here) in Ω. For tetrahedral elements, a smoothing domain Ωks is created by connecting three nodes of its associated face to the central point of two adjacent tetrahedral elements as shown in Figure 1.Figure 1
Domain discretization and a smoothing domain (light pink) associated with a face (dark brown) in the FS-FEM.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:(15)B¯k=1Vk∑j=1Ne,k14Vk,jBk,j,j∈1,Ne,k,where Vk is the volume of smoothing domain k, N(e,k) is the number of elements, V(k,j) is the volume of jth element in smoothing domain k, and B(k,j), a 6×12 matrix, is the standard compatible strain-displacement matrix of element j in smoothing domain k. B¯k=[B¯(k,1)⋯B¯N(n,k)] is a 6×3N(n,k) matrix. N(n,k)=4 for boundary faces and N(n,k)=5 for inner faces.The smoothed strain is written as(16)ε¯k=B¯kdk,where dk=[u1,v1,w1,…,uN(n,k),vN(n,k),wN(n,k)]T is the collection of the displacements of nodes in smoothing domain k.Then, the local smoothed stiffness matrixK¯k is given as follows:(17)K¯k=VkB¯kTDB¯k.The entries ofK¯k are defined by(18)K¯k,ij=VkB¯k,iTDB¯k,j.The local smoothed stiffness matrices are assembled to produce the entry of the global smoothed stiffness matrix K¯ij:(19)K¯ij=∑k=1NfK¯k,ij.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:(20)R=GF,where GF is an operator extracting rotational matrix from the deformation gradient matrix F. F=x01x02x03x010x020x030-1, xij=xj-xi, xi is the deformed position of node i in element e, and xi0 is the undeformed position of node i in element e. There are three ways to define operator G 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 domainR¯ks, 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 N(e,k) elements. For smoothing domains associated with boundary faces, R¯ks equals the rotational matrix of their unique associated element. For smoothing domains associated with inner faces, it takes four steps to get R¯ks: 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 j in the smoothing domain k is S(k,j)=e1e2e3. e1 is x^100, e2 is x^200, and e3 is e2×e1. x^ is the normalized x. The projection matrix is defined as Prj=S(k,j)(S(k,r))T, where S(k,r) is the selected reference local coordinate system. Then the element-based rotational matrix is projected by (21)R~k,j=Rk,jPrj-1,r,j∈1,Ne,k.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 R~(k,j) is first transformed to quaternion q(k,j).Thenslerp is applied to get the interpolated quaternion qk:(22)qk=slerpqk,r,qk,j;λ=qk,rqk,r-1qk,jλ,where the interpolation coefficient is λ=V(e,j)/Vk.In the fourth step,q is transformed back to a matrix R¯ks.Applying the rotational matrix on each smoothing domain, we get the elastic forcesfk acting on the nodes by a smoothing domain:(23)fk=RkK¯kRk-1xk-xk0=RkK¯kRk-1xk-RkK¯kxk0=RkK¯kRk-1xk+Rkfk0,where xk=[x1,y1,z1,…,xN(n,k),yN(n,k),zN(n,k)]T and xk0=[x10,y10,z10,…,xN(n,k)0,yN(n,k)0,zN(n,k)0]T are the deformed and undeformed nodal positions, respectively. Rk is a 3N(n,k)×3N(n,k) matrix that contains N(n,k) copies of the 3×3 matrix R¯ks along its diagonal. The force offset vector fk0=-K¯kxk0 is invariant and can be precomputed to accelerate the algorithm. By (23), the rotational part of deformation is cancelled by (Rk)-1. The elastic forces are computed in the local coordinate of smoothing domain k and then rotated back to the global coordinate by Rk.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 by(24)f′=K¯′x+f0′,where the global stiffness matrix K¯′ is the summation of the smoothing domain’s rotated stiffness matrix K¯k′=RkK¯k(Rk)-1 and the global force offset vector f0′ is the summation of the smoothing domain’s force offset vector fk0′=Rkfk0.


### 4.3. Dynamics

The following governing equation describes the dynamics of the system:(25)Mx¨t+Cx˙t+f′=fext,where xt=[x1,y1,z1,…,xN(n),yN(n),zN(n)]T is the current object state; x¨(t) and x˙(t) are the first and second derivatives of x(t) with respect to time. N(n) is the total number of nodes in the entire mesh. M is the 3N(n)×3N(n) mass matrix and C is the 3N(n)×3N(n) damping matrix. The lumped mass matrix is used for efficiency here. Rayleigh damping is used to compute the damping matrix C=αM+βK, where α and β are known as the mass damping and stiffness damping coefficients, respectively. fext is the 3N(n)×1 external loads vector. We omit the time parameter of x(t) using x instead to lighten the notation.Equation (25) is discretized in time domain and solved using implicit Euler. Substituting (24), vi+1=x˙i+1, v˙i+1=x¨i+1, and xi+1=xi+Δtvi+1 into (25), we get a group of linear algebraic equations:(26)M+ΔtC+Δt2K¯′vi+1=Mvi-ΔtK¯′xi+f0′-fext,where v is the velocity vector, Δt is the time step, and i=t/Δt stands for the ith time step.


### 4.4. Dynamic Simulation Algorithm

The whole dynamic simulation process is summed up in the following pseudocode in Algorithm1. 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 B¯k in line (5), the computation of smoothed stiffness matrix K¯k in line (6), and the computation of offset force f0′ 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).Algorithm 1:The simulation procedure.
(1) initializex0, v0
(2)forall element e compute the standard strain-
displacement matrixB
(3) build smoothing domains
(4)forall smoothing domain k
(5)  computeB-k using Eq. (15)
(6)  computeK-k using Eq. (17)
(7)endfor
(8)i←0
(9)loop
(10)forall element e compute R using Eq. (20)
(11)forall smoothing domain k compute Rk
(12) assemble matrixK-′←∑k‍RkK-k(Rk)-1
(13) assemble vectorf0′←-∑k‍RkK-kxk0
(14) apply external forcefext
(15)A←M+ΔtC+Δt2K-′
(16)b←Mvi-Δt(K-′xi+f0′-fext)
(17) solvevi+1 from Avi+1=b
(18)xi+1←xi+Δtvi+1
(19)i←i+1
(20)endloop


## 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 Ω=∑k=1Nf‍Ωk and Ωp∩Ωq≠∅, p≠q, in which Nf is the total number of smoothing domains (i.e., faces here) in Ω. For tetrahedral elements, a smoothing domain Ωks is created by connecting three nodes of its associated face to the central point of two adjacent tetrahedral elements as shown in Figure 1.Figure 1
Domain discretization and a smoothing domain (light pink) associated with a face (dark brown) in the FS-FEM.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:(15)B¯k=1Vk∑j=1Ne,k14Vk,jBk,j,j∈1,Ne,k,where Vk is the volume of smoothing domain k, N(e,k) is the number of elements, V(k,j) is the volume of jth element in smoothing domain k, and B(k,j), a 6×12 matrix, is the standard compatible strain-displacement matrix of element j in smoothing domain k. B¯k=[B¯(k,1)⋯B¯N(n,k)] is a 6×3N(n,k) matrix. N(n,k)=4 for boundary faces and N(n,k)=5 for inner faces.The smoothed strain is written as(16)ε¯k=B¯kdk,where dk=[u1,v1,w1,…,uN(n,k),vN(n,k),wN(n,k)]T is the collection of the displacements of nodes in smoothing domain k.Then, the local smoothed stiffness matrixK¯k is given as follows:(17)K¯k=VkB¯kTDB¯k.The entries ofK¯k are defined by(18)K¯k,ij=VkB¯k,iTDB¯k,j.The local smoothed stiffness matrices are assembled to produce the entry of the global smoothed stiffness matrix K¯ij:(19)K¯ij=∑k=1NfK¯k,ij.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:(20)R=GF,where GF is an operator extracting rotational matrix from the deformation gradient matrix F. F=x01x02x03x010x020x030-1, xij=xj-xi, xi is the deformed position of node i in element e, and xi0 is the undeformed position of node i in element e. There are three ways to define operator G 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 domainR¯ks, 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 N(e,k) elements. For smoothing domains associated with boundary faces, R¯ks equals the rotational matrix of their unique associated element. For smoothing domains associated with inner faces, it takes four steps to get R¯ks: 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 j in the smoothing domain k is S(k,j)=e1e2e3. e1 is x^100, e2 is x^200, and e3 is e2×e1. x^ is the normalized x. The projection matrix is defined as Prj=S(k,j)(S(k,r))T, where S(k,r) is the selected reference local coordinate system. Then the element-based rotational matrix is projected by (21)R~k,j=Rk,jPrj-1,r,j∈1,Ne,k.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 R~(k,j) is first transformed to quaternion q(k,j).Thenslerp is applied to get the interpolated quaternion qk:(22)qk=slerpqk,r,qk,j;λ=qk,rqk,r-1qk,jλ,where the interpolation coefficient is λ=V(e,j)/Vk.In the fourth step,q is transformed back to a matrix R¯ks.Applying the rotational matrix on each smoothing domain, we get the elastic forcesfk acting on the nodes by a smoothing domain:(23)fk=RkK¯kRk-1xk-xk0=RkK¯kRk-1xk-RkK¯kxk0=RkK¯kRk-1xk+Rkfk0,where xk=[x1,y1,z1,…,xN(n,k),yN(n,k),zN(n,k)]T and xk0=[x10,y10,z10,…,xN(n,k)0,yN(n,k)0,zN(n,k)0]T are the deformed and undeformed nodal positions, respectively. Rk is a 3N(n,k)×3N(n,k) matrix that contains N(n,k) copies of the 3×3 matrix R¯ks along its diagonal. The force offset vector fk0=-K¯kxk0 is invariant and can be precomputed to accelerate the algorithm. By (23), the rotational part of deformation is cancelled by (Rk)-1. The elastic forces are computed in the local coordinate of smoothing domain k and then rotated back to the global coordinate by Rk.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 by(24)f′=K¯′x+f0′,where the global stiffness matrix K¯′ is the summation of the smoothing domain’s rotated stiffness matrix K¯k′=RkK¯k(Rk)-1 and the global force offset vector f0′ is the summation of the smoothing domain’s force offset vector fk0′=Rkfk0.


## 4.3. Dynamics

The following governing equation describes the dynamics of the system:(25)Mx¨t+Cx˙t+f′=fext,where xt=[x1,y1,z1,…,xN(n),yN(n),zN(n)]T is the current object state; x¨(t) and x˙(t) are the first and second derivatives of x(t) with respect to time. N(n) is the total number of nodes in the entire mesh. M is the 3N(n)×3N(n) mass matrix and C is the 3N(n)×3N(n) damping matrix. The lumped mass matrix is used for efficiency here. Rayleigh damping is used to compute the damping matrix C=αM+βK, where α and β are known as the mass damping and stiffness damping coefficients, respectively. fext is the 3N(n)×1 external loads vector. We omit the time parameter of x(t) using x instead to lighten the notation.Equation (25) is discretized in time domain and solved using implicit Euler. Substituting (24), vi+1=x˙i+1, v˙i+1=x¨i+1, and xi+1=xi+Δtvi+1 into (25), we get a group of linear algebraic equations:(26)M+ΔtC+Δt2K¯′vi+1=Mvi-ΔtK¯′xi+f0′-fext,where v is the velocity vector, Δt is the time step, and i=t/Δt stands for the ith time step.


## 4.4. Dynamic Simulation Algorithm

The whole dynamic simulation process is summed up in the following pseudocode in Algorithm1. 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 B¯k in line (5), the computation of smoothed stiffness matrix K¯k in line (6), and the computation of offset force f0′ 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).Algorithm 1:The simulation procedure.
(1) initializex0, v0
(2)forall element e compute the standard strain-
displacement matrixB
(3) build smoothing domains
(4)forall smoothing domain k
(5)  computeB-k using Eq. (15)
(6)  computeK-k using Eq. (17)
(7)endfor
(8)i←0
(9)loop
(10)forall element e compute R using Eq. (20)
(11)forall smoothing domain k compute Rk
(12) assemble matrixK-′←∑k‍RkK-k(Rk)-1
(13) assemble vectorf0′←-∑k‍RkK-kxk0
(14) apply external forcefext
(15)A←M+ΔtC+Δt2K-′
(16)b←Mvi-Δt(K-′xi+f0′-fext)
(17) solvevi+1 from Avi+1=b
(18)xi+1←xi+Δtvi+1
(19)i←i+1
(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 v=0.33 and elasticity modulus E=1×106 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 K¯k and rotational matrix Rk. Because the local stiffness matrix K¯k 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 Rk.Figure 2
Let execution time be a function of elements.We summarized the simulation time distribution in each step in Table1. The table shows that CS-FEM spends less than 2% extra time (tIR) computing smoothing domain-based rotational matrix Rk. The bottleneck of CS-FEM is linear algebraic equations assembly (tA) and solving (teq). If the mesh is less than 6 K tetrahedron, CS-FEM reaches 60 fps and can be used in real-time applications.Table 1
The simulation time statistics: vertex numberV, tetrahedron number T, face or smoothing domain number F, preprocessing time for generating smoothing domain and computing smoothed stiffness matrix tpre, smoothing domain-based rotational matrix computation time tIR in each iteration, stiffness matrix assembly time tA and linear algebraic equation solving time teq in each iteration, the average computation time tstep in an iteration, and the frame per second fps.
Number
V [K]
T [K]
F [K]
t
p
r
e [ms]
t
R [ms]
t
I
R [ms]
t
A [ms]
t
e
q [ms]
t
s
t
e
p [ms]
fps[1]
(1)
0.16
0.41
0.94
0.52
0.02 (1.71%)
0.01 (0.43%)
0.65
0.72
1.43
701.08
(2)
0.93
3.24
6.984
1.07
0.29 (3.82%)
0.13 (1.78%)
4.05
2.88
7.50
133.42
(3)
2.80
10.94
23.00
4.32
0.83 (3.29%)
0.24 (0.96%)
14.31
9.23
25.12
39.81
(4)
6.25
25.92
53.86
11.13
1.52 (2.71%)
0.67 (1.19%)
33.46
18.79
56.08
17.83
(5)
11.77
50.62
104.40
20.16
2.93 (2.74%)
1.46 (1.37%)
65.95
33.70
107.13
9.33Mesh 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 ee:(27)ee=E~-E1/2,E~=uTK¯u2,where E~ is the numerical solution of strain energy and E is the exact solution of strain energy. The displacement error norm of node i was defined by(28)edi=u~i-uiui×100%,where u~i is the numerical solution of displacement and ui 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 (1.0×1.0×1.0) 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 v=0.25, elasticity modulus E=1.0, and pressure p=1.0. 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 E=0.9486 and the deflection at point B (1.0,1.0,0.5) is 3.3912.Figure 3
A 3D cantilever subjected to a uniform pressure using distorted meshes withαir = 0.The experiment was taken using five distorted meshes. Those meshes were generated using an irregular factorαir between 0.0 and 0.49, random numbers rx, ry, and rz between -1.0 and 1.0, and the edge lengths Δx, Δy, and Δz of a grid in a cubic mesh in the following fashion: (29)x=x+Δx∗rx∗αir,y=y+Δy∗ry∗αir,z=z+Δz∗rz∗αir.The meshes were generated using distortion coefficient αir from 0.0 to 0.4.Figures4(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 (MT10/4) [71], but the displacement of CS-FEM is more accurate than that of MT10/4. 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.Figure 4
Displacement and energy error norm of a 3D cubic cantilever subjected to a uniform pressure versus the distortion coefficients.
(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 (0.9m × 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 E=0.4×106 Pa and v=0.33. 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 (0.9,0.0,0.3) at time 0.25 s is 16.8977 cm. The experiment was taken using five distorted meshes with αir= 0.0 to 0.4. The mesh plotted in Figure 6 is the distorted mesh generated with αir=0.4.The configuration att = 0.25 s generated by CS-FEM is rendered in Figure 7. The transient vertical displacements of node A were plotted in Figure 5.Figure 5
Transient vertical displacements of a 3D cantilever beam subjected to gravity.
(a)
L-FEM
(b)
C-FEM
(c)
LFS-FEM
(d)
NLFS-FEM
(e)
CS-FEMFigure 6
A mesh of a 3D cantilever beam generated withαir = 0.4.Figure 7
The initial and final configuration (t = 0.25 s) of the 3D cantilever beam subjected to gravity using CS-FEM.The simulation frequencies of L-FEM and C-FEM begin to deviate from those of the regular mesh (αir=0) if αir≥0.3 (lines in purple and green) and were as low as one-third of those of the regular mesh if αir=0.4 (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 αir=0.4. 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.Table 2
Vertical displacement of A in irregular meshes at time 0.25 s with varying distortion coefficients. The number in bracket stands for the relative error (%) between the numerical results atαir>0.0 and those at αir=0.0. The reference solution of vertical displacement of A is 16.8977 cm.
αir=0.0
αir=0.1
α
i
r
=
0.2
αir=0.3
αir=0.4
CS-FEM
17.32 cm
17.23 cm (0.51%)
17.06 cm (1.50%)
16.70 cm (3.58%)
16.3 cm (5.89%)
LFS-FEM
17.75 cm
17.65 cm (0.56%)
17.49 cm (1.46%)
17.11 cm (3.60%)
16.54 cm (6.82%)
NLFS-FEM
16.74 cm
16.71 cm (0.17%)
16.49 cm (1.49%)
16.34 cm (2.38%)
15.95 cm (4.72%)
L-FEM
16.33 cm
16.23 cm (0.61%)
15.94 cm (2.39%)
15.3 cm (6.31%)
5.77 cm (64.67%)
C-FEM
15.99 cm
15.89 cm (0.63%)
15.63 cm (2.25%)
15.14 cm (5.31%)
5.81 cm (63.66%)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 v=0.33 and E=1×105 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.Figure 8
Simulation results of beams subjected to gravity in timet=0 s (a), 0.5 s (b), and 1.0 s (c) using LFS-FEM (in blue), C-FEM (in green), L-FEM (in purple), and CS-FEM (in red).
(a)
(b)
(c)Figure 9
Volume gains under rotational deformation.
(a)
C-FEM
(b)
CS-FEM
(c)
L-FEM
(d)
LFS-FEMSoftening 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 (0,0,-50 N). The bridge was discretized using a mesh including 3,923 nodes and 8,680 tetrahedral elements with v=0.45, E=5×106 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.Figure 10
Simulation results of a bridge subjected to gravity and a static load.Figure 11
Comparison of strain energy obtained from CS-FEM and C-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.



---

*Source: 1467356-2017-03-29.xml*