Development of an adaptive explicit algorithm for static simulation using the vector form intrinsic finite element method

Thevectorformintrinsicfiniteelement(VFIFE)methodisasolutiontechniquefornonlinearstructuralproblems,whichdescribesacontinuous body using a set of particles instead of a mathematical function. Thus, a dynamic particle equation can be established by Newton’s law of motion, and a viscous or kinetic damping can be introduced to obtain the steady state of the structure. This paper focuses mainly on the development of a stability condition regarding the explicit central difference method used in VFIFE to guarantee the system’s convergence. The process is established and evaluated in combination with a dynamic relaxation method with kinetic damping and discrete control theory. Four numerical examples of structure nonlinear problems are used to verify the accuracy, stability and efficiency of the method.


INTRODUCTION
Nowadays, solving nonlinear structural problems frequently involves using a numerical method, such as the finite element method (FEM) or the finite difference method. In the FEM, the object is discretized into elements such that the element displacement vector and the force vector are related by an element stiffness matrix. However, numerically manipulating structural stiffness often causes stability and convergence issues in the procedure. For geometrically nonlinear problems, the solution is affected by rigid body motion; therefore, a geometric stiffness must be established to eliminate in-equilibrium force. Thus, alternative methods, such as the updated Lagrangian formulation, total Lagrangian formulation [1] and co-rotational formulation [2], have been proposed and investigated. Regardless, certain situations, such as fracturing, flexible multibody, or highly nonlinear problems, can lead to challenges. To overcome these issues, Day [3] proposed the dynamic relaxation method (DRM), which converts a static system to an artificially dynamic problem and has better computational efficiency. The DRM has been applied to a wide range of problems, including highly nonlinear problems [4], form-finding analysis [5][6][7][8], the post-buckling problem [9,10] and fracture mechanics [11,12]. The explicit time integration method in the DRM is conditionally stable, and its time step size is narrowly restricted. When the structure triggers high-frequency oscillations, the time step size must decrease to an extremely small value, leading to a low convergence rate. Thus, many studies have aimed to improve the convergence rate of DRM [13][14][15][16]. Extensive research has also been conducted to enhance the effectiveness of the DRM's time integration scheme. For example, Ramirez [17] derived the exact coefficient of the discrete transfer function for the Newmark family integration algorithm. Chang [18] proposed an unconditionally stable explicit pseudodynamic algorithm where the presence of high-frequency modes can be performed without the problem from the use of a very small time step. Chen and Ricles [19] developed a direct integration algorithm of unconditional stability based on the discrete transfer function approach. Alamatian [20] used the Gershgorin circle theorem to improve the convergence rate while guaranteeing the numerical stability of DRM.
Ting et al. [21] proposed the vector form intrinsic finite element (VFIFE) method to solve nonlinear structural problems. In VFIFE, the object is represented by a group of particles with mass/inertia and the motion of the particle is described within a path unit. The element that connects particles carries the material and geometric characteristics. Thus, the resistant forces inside the element can be evaluated by the deformation of the element within a path unit. Consequently, the motion of each particle can be formulated by vector mechanics using Newton's second law of motion without invoking a continuum governing equation. The VFIFE approach has been documented [22] and applied to planar elements [23], a two-dimensional (2D) frame structure [24], a three-dimensional (3D) membrane [25], 3D solids [26] and 3D frame [27]. In the analysis of mechanics of structure, both the stability and efficiency are crucial indicators for the algorithm. In light of this necessity, the goal of this work is to investigate the stability and convergence of the integration scheme used in VFIFE, namely the central difference algorithm, where the stability and efficiency were rarely discussed. We conduct stability analysis in combination with the DRM and discrete control theory. We derive conditions of stability via the discrete transfer function. Finally, we provide several numerical examples using frame element to verify our work and show the effectiveness of the approach in solving nonlinear problems.

MATHE MATICAL BACKGROUND
The VFIFE algorithm describes the motion of a continuous body by using a finite number of points (or particles) associated with the theory of vector form analysis. The description of point values includes the allocation of mass points to the structure and definition for paths of point motions. Next, constitutive conditions are selected and assigned to connect the mass points in deformable structures that yield the formulation. Finally, along each path element, the motion for each free particle is established by Newton's law. In this section, we provide a brief mathematical background for VFIFE; for more details, see [21,22,28].

Point value description
As shown in Fig. 1, a set of points representing the configuration of a frame structure is displaced from one position to the other position, where m denotes the mass of an arbitrary point in the structure and x denotes the position of the point of interest. Assume that the interaction among points is modeled by a beam element. The properties of each element can be characterized by nodal displacements and corresponding nodal forces. Further, the element displacement functions must fulfill continuity conditions. In other words, the element nodal displacement is determined by the particle motion, and the resistant forces induced by element deformations are represented by element internal nodal forces that are applied to those particles connected with the element. The element lumps the mass at node points and does not possess mass. On the other hand, external forces may come from either (1) concentrated forces directly applied to the particle of interest or (2) forces on the element that must be transformed equivalently to nodal external forces.

Discrete path and governing equation
Referring to Fig. 1, assume the particle in the structure undergoes a motion from time t a to t whose path is called the path element and within which the deformation of the structure is infinitesimal. Then, it is feasible to apply infinitesimal strain and engineering stress to evaluate continuum stress and compute virtual work. The formulation hereafter comprises the process of displacement and deformation, the calculation of equivalent internal forces and external forces, and the solving of the equation of motion for each particle.
First, within each path element, the motion of a particle can be described by Newton's law of motion: where d is the displacement vector including translational and rotational displacements of the particle, F ext denotes the external forces, including the concentrated external force/moment vectors applied at the particle and equivalent nodal force/moment vectors from the load on the element, f i int denotes the equivalent nodal forces, moments due to axial deformation and rotations at the two nodes of the element, n e is the number of elements incident with the particle of interest and M is the mass matrix whose elements contain the mass and mass moment of inertia of the particle.

Evaluation of deformation and internal forces
VFIFE includes a technique to evaluate the deformation of an element. Figure 2a shows a two-node frame element (1, 2) at two positions within a time segment t a -t. The relative change for nodal displacements (u 1 , u 2 ) and nodal angles (β 1 , β 2 ) can be calculated as where (x i , θ i ) are the position and angle of node i at time t and (x ia , θ ia ) are the position and angle of node i at time t a . The nodal positions, nodal internal forces and properties of the element must be prescribed at time t a . It is assumed in VFIFE that the deformation of the element is so small that the element's geometric shape remains unchanged during time segment t a -t. Thus, the internal forces in the beam element can be calculated via beam theory.
To further realize the motion of the element within the time segment, we first deal with the rigid body motion of the element. As shown in Fig. 2b, a local coordinate system (e a 1 , e a 2 , e a 3 ), which is parallel to the principal coordinate system, is established on (1 a , 2 a ) and a local coordinate system (e t 1 , e t 2 , e t 3 ) is established on (1 t , 2 t ). Note that the two coordinate systems can be related by a rotation matrix, as will be discussed later. The rigid body motion of the element may involve translation and rotation. The translation can be calculated by Eq. (2a), while the rotation may be viewed as the summation of the rotation of the element about the axis e a 1 and the rotation about the axis perpendicular to the plane formed by e a 1 and e t 1 . Thus, the rigid body rotation, γ t , for the element can be written as where β 1 x = β 1 · e a 1 is the value of the nodal angle β 1 projected on the axis e a 1 and As the element (1, 2) undergoes a rigid body rotation motion about the rotation axis, γ t , during the time segment, a rotation matrix, R t a , of the coordinate system (e a 1 , e a 2 , e a 3 ) with respect to the coordinate system (e t 1 , e t 2 , e t 3 ) can be established in terms of the parameters of γ t via Rodrigues' formula (see Appendix A). Similarly, a transformation matrix, a , that relates the global coordinate system and the local coordinate system (e a 1 , e a 2 , e a 3 ) can be established such that To determine the deformation and the relationship between stress and strain, we assume a reverse motion of the element from time t back to time t a , including translation and rotation. Referring to Fig. 3, the element (1, 2) first undergoes a reverse translation (−u 1 ) such that node 1 t coincides with 1 a . The relative displacement between the two nodes is Second, the amount of reverse rotation from state t to state t a is exactly the value described by Eq. (3), but in the opposite direction.
After the reverse motion, the nodal deformation can be obtained by comparing the position and orientation in the reverse state and those in the t a state. Thus, the change of translational displacement can be obtained as where η 2d is the axial deformation vector along the beam axis and η 2r is the displacement due to reverse rigid body motion. Further, the magnitude of η 2d can be obtained as whereˆ e represents the axial stretch of the element, and l t and l a are the element lengths at times t and t a , respectively. Additionally, the deformation of each nodal angleβ id (i = 1, 2) during the time segment can be obtained aŝ whereβ i is total nodal angle i described by Eq. (2). After denotingγ t = [β 1 x ,θ at,y ,θ at,z ] T ,β id can be expressed aŝ Thus, once the nodal angleβ i and angle for rigid body motionγ t are known, we can obtain the deformation of the nodal angles. Note that the above relations are calculated in the local coordinate system (e a 1 , e a 2 , e a 3 ).

Internal force
The determination of element internal forces proposed by VFIFE utilizes the deformations rather than the total displacement, which is used in the traditional FEM. A coordinate system, called the deformation coordinates, is defined for the shape functions. Then, the increment of internal forces and moments due to deformation can be derived via virtual work. A detailed derivation can be found in [21]. The results are shown in the following equation: where E a is Young's modulus, G a is the shear modulus, A a is the area of the cross-section,Ĵ ya andĴ za are the area moment of inertia andĴ xa is the polar moment of inertia. The remaining increment of force/moment components can be obtained via static equilibrium condition [21]. After the increments of internal forces and moments are obtained, the internal forcesf i and momentm i are calculated with respect to the local coordinate system (e a 1 , e a 2 , e a 3 ) at a fictitious state, i.e.
Again, the results must be transformed to the global coordinate system at time t by

Time integration
To solve the equation of motion of particles within a path element in Eq. (1), it is necessary to utilize the techniques of time integration.
In this work, we apply the central difference explicit scheme. Further, if a static solution is required, the mass damping scheme associated with a mass damping force f d = −ζ c mḋ(t ) can be included in the right-hand side of Eq. (1). Derivation of the displacement equation from Eq. (1) can be referred to [21,28]. The resulting scalar equation for the displacement can be written as where ζ c is the mass damping coefficient and the subscript "n" represents the nth time step. Alternatively, if we wish to apply DRM with kinetic damping, the kinetic energy of the particle during the time segment must be derived. In such an approach, the static equilibrium corresponds to the condition where the potential energy reaches a minimal value. Let the displacement equation without damping (ζ c = 0) be written as The nodal velocities for the next step can be represented by the half-interval method aṡ Then, the minimum potential energy can be detected by checking whether the inequality below holds at each instance: If it holds, then setḋ 2 n−1/2 = 0 to obtain d n = d n−1 . This yields a new displacement: By alternately iterating Eqs (17) and (21) in accordance with the condition of kinetic energy, we can achieve the static solution of a dynamic system.

STABLE E XPLICIT ALGORITHM VIA CONTROL THEORY 3.1 Properties of transfer functions
The stability of a system can be observed from the pole position of the system's transfer function, as adopted from control theory. Many studies have applied this concept to develop new time integration algorithms with good stability [14][15][16]. Here, we develop an algorithm based on the central difference scheme associated with the discrete control theory, which guarantees the convergence of the integration. Consider a second-order, single-degree-of-freedom (SDOF) system under external force; the equation can be expressed as The SDOF system will have the following transfer function G(s) in the Laplace domain: where ω n is the natural frequency, ζ is the damping ratio, s is the variable in the Laplace domain, and X(s) and F(s) are the Laplace transforms of d(t) and f(t), respectively. Then, the poles of the transfer function can be written as The system is stable only if the real part of the poles is negative-valued or located in the left half plane of the complex plane. The stability for the difference equation can also be evaluated via the above concept, but it needs to be represented in a discrete system via a z-transform. For a discrete function f(k), the z-transform is defined as [29] where f(k) is a sample of f(t) and k is the discrete sample time. By shifting theorem, we have where m is an integer. Without loss of generality, the equation of motion can be expressed in terms of displacement equation in difference form as Applying Eq. (26) to Eq. (27a) and letting the z-transform of d n be X(z), Eq. (27a) becomes Thus, the general form of discrete transfer function G(z), which is the ratio of X(z) to F(z), can be written as where F(z) is the z-transform f n , and (B 2 , B 1 , B 0 ) and (A 2 , A 1 , A 0 ) are coefficients of the numerator and denominator of G(z), respectively. Detailed derivations can also be referred to [17].

Consistency of the continuous system and discrete system
In control theory, to ensure the conversion of the discrete system from the continuous system, two systems must follow a discretization method called "matched pole and zero". The equation that can relate the characteristics in the z-domain to those in the s-domain can be written as z = e s t .
Rather than using the exponential function, Tustin's method [30] is used to approximate the discretization for the continuous system by relating the z and s domains, namely This can also be written as After substituting Eq. (29c) into Eq. (23), the following discrete transfer function corresponding to the continuous transfer function can be obtained: .
The poles of the transfer function of Eq. (30) can be obtained by solving the characteristic equation of the transfer function, i.e. (4 + ω 2 n t 2 + 4ζ ω n t )z 2 + (−8 + 2ω 2 n t 2 )z + (4 − 4ζ ω n t + ω 2 n t 2 ) = 0. (31a) This yields In control theory, to ensure the stability of the discrete system, the poles of the corresponding transfer function need to lie on or inside the unit circle of the discrete domain.

Development of stability for central difference algorithms
In this study, the integration algorithm to evaluate the central time difference is designed by assuming the following variation in velocity and acceleration. The procedure starts with the following assumption in the nth time step: where α 1 and α 2 are two parameters to be determined. Also note that the velocity and acceleration can be, respectively, expressed aṡ

Numerator Denominator
The kinematic quantities (displacement, velocity and acceleration) and external force are then associated with different time steps, i.e. n and (n + 1). Substituting Eqs (33) and (34) into Eq. (32) to eliminate velocity and acceleration and transforming those quantities into the z-domain using Eq. (25) yields where the coefficients of the numerator and denominator for central difference are listed in Table 1.
Equating the poles of Eq. (35) to Eq. (31b) and solving for α 1 and α 2 yields Once α 1 and α 2 are specified according to t and ω n , the stability of the integration algorithm can be preserved. Conversely, if α 1 and α 2 are prespecified, the mass for the kinetic DRM can be adjusted from Eq. (37) as Since the central difference algorithm is an explicit time integration method, it is conditionally stable. Equation (38) may also imply that to achieve the stability in a dynamic system, the fictitious mass/inertia can be adjusted adaptively with the material stiffness k, time step t and the coefficient α 2 .

Stability verification by an SDOF system
The stability of the time integration can be usually verified by using an SDOF system. Let the system equation be expressed in statespace form as where A is the amplification matrix, u n is the input vector, B is the input matrix, and X n+1 and X n are state vectors at the (n + 1)th and (n)th time steps, respectively. The stability criterion can be determined by the eigenvalues of the matrix A, as the roots of the characteristic polynomial derived from the determinant of (A − λI): where I is the identity matrix and λ is the eigenvalue. The eigenvalues of matrix A are the same as the poles of the discrete transfer function. Hence, the properties of the poles can also be applied to those of eigenvalues of matrix A. Therefore, with the central difference algorithm, the displacement equation, Eq. (32), having no external load can be written as Combining Eq. (41) and identity equation d n = d n into the state-space form yields The amplification matrix A for the system can thus be obtained as The two eigenvalues for A are Consider the system response in the period of T = 1 s. The parameter α 2 is defined by various time steps when the control scheme is used. Figure 4a shows the comparison of magnitude of the eigenvalue with and without control scheme. It can be seen that the system diverges when the time step is larger than 0.318 if the calculation uses the traditional central difference where α 2 is not adaptively adjusted. Figure 4b shows the positions of poles of the system with various time steps. It can be seen that the poles lie on or inside the  unit circle only if the time step is smaller than 0.318. However, if the time step is larger than 0.318, the poles can be retained on the unit circle by using the control scheme as the green diamond shows.

Fictitious mass and inertia in VFIFE
In FEM, the deformation requires the overall structure's system stiffness matrix due to the interaction of internal forces between the members; thus, the analysis can be thought for a multiple-degree-of-freedom system. Since VFFIE is point value described, the evaluation of motion for mass point is an independent difference equation when using the central difference, and the resultant has been assembled on the node. Therefore, a mass point can be considered as an SDOF in VFIFE. Nonetheless, to satisfy the stability criterion, the element stiffness must be known according to Eq. (38). Rather than considering the system stiffness matrix as in the traditional FEM, we focus on the axial stiffness term for each node during the internal force evaluation, i.e.
We choose the axial stiffness because it is more dominant than other rotational stiffness. Choosing the largest stiffness for evaluation can avoid divergence if high frequencies are triggered. As shown in Fig. 5, where a node is connected to multiple elements, the calculation of fictitious mass has to evaluate the equivalent stiffness on the connecting node. Instead of assembling the stiffness on the node, the mass of each element can be evaluated separately by Eq. (38).
Since the element mass has been calculated, the nodal mass to be adjusted becomes If there are multiple elements incident to one node, the mass of the node can be summed up as where n e is the number of elements incident to the node of interest. Similarly, the mass moment of inertia of the node is whereÎ # denotes the nodal mass moment of inertia, r # is the radius of gyration and the subscript represents the axis to which the quantity refers. If there are multiple elements incident to one node, the total moment of inertia can also be summed as In summary, the algorithm of the VFIFE associated with kinetic DRM can be developed as shown in Fig. 6.

NUMERICAL E X A MPLES 4.1 Cantilever beam with a transversal load
The large deflection analysis of a planar cantilever beam has been favorably used to test the geometrically nonlinear issue by Mattiasson [31], Barten [32], and Bisshopp and Drucker [33]. As shown in Fig. 7, a cantilever beam of length L is subject to a vertical force P. The loading curve is also shown and its pattern is used throughout all examples. For comparison, the material properties are assumed as E = 10, L = 10 and I = 10, such that the ratios can be identical with those in [31]. The axial deformations are ignored in this problem, so the cross-sectional area is set to a large value in the simulation. There are 11 nodes and 10 elements in the beam, and we   use a constant time step t = 1 × 10 -3 . Assume the density ρ = 0.01 is used for mass damping and kinetic DRM without a control strategy. The kinetic DRM and mass damping are used to simulate the transversal (W) and axial (U) displacements at the free end.
The result is compared with the literature [31]. Figure 8 shows the free-end displacement curve where kinetic damping was applied with the discrete control algorithm of α 2 = 0.49. The results are perfectly in agreement with [31]. Thus, VFIFE with kinetic DRM can handle this problem effectively. Figure 9 shows the comparison of various numerical schemes, excluding mass damping because its value is out of the figure range. The result by kinetic damping with a control strategy, α 2 = 0.49 and α 2 = 0.41, shows a better  convergence rate than without control schemes. We see that the efficiency of the scheme can be improved by properly selecting the coefficient in the kinetic damping with a control strategy. Figure 10a and b, respectively, shows the time history of positions x and y subject to P = 10 for various schemes, which support the result shown in Fig. 9. Figure 11 shows a cantilever 45°bend subject to a concentrated end load where the loading curve is the same as in Fig. 7. The bend has a radius of 100 with a unit square cross-section, and it lies in the x-z plane. The tip load is applied in the y-direction. The material properties are E = 10 7 , J y = J z = 1/12 and G = 0.5 × 10 7 , which are used in the literature [34]. The beam is meshed with nine nodes and eight elements. Assume the density ρ = 1 is used for mass damping and kinetic DRM without control strategy. Various numerical schemes, including kinetic damping with/without control theory and mass damping, are used to obtain tip displacement. Figure 12 shows the results of tip displacements where kinetic damping was applied with the discrete control algorithm of α 2 = 0.49 and t = 1 × 10 -3 . The results show good agreement with those of [34]. Figure 13 compares the number of iterations among the three schemes. Figure 14a shows the comparisons of convergence with different time steps t and numerical schemes, and Fig. 14b shows the variations of kinetic energy for kinetic damping with various control schemes. The results show that the computation diverges as t increases for the scheme without a control strategy. The convergence is tested for different strategies, and the results are listed in

Figure 15
Schematic of space truss (unit of length: mm). Table 2. As shown in the table, the divergence occurs when the time step t = 1 × 10 -2 without a control strategy. Nonetheless, the scheme with a control strategy reveals more robustness in stability than the scheme without a control strategy. Figure 15 shows a space truss whose configuration is taken from [35]. However, here we apply a vertical displacement at the top of the dome as the input rather than a force input at the top in [35]. The truss is made from a material with Young's modulus E = 2.06 × 10 11 N/m 2 , where the area section of each member A = 2.16 × 10 -4 m 2 and ρ = 7800 kg/m 3 . This structure has 73 nodes and 168 elements, where all support conditions are pinned. This example demonstrates that the proposed algorithm can easily handle the structure with multiple elements pinned at one node. The evaluation of internal forces for the truss element is the same as the beam element except that it considers no bending moment in the element and, hence, ignores the mass moment inertia of the node. A constant time step t = 1 × 10 -4 s is used throughout this example unless stated otherwise. Figure 16 shows the force-displacement curve for the top node in the z-direction where kinetic damping was applied with a control coefficient α 2 = 0.49. Figure 17 compares iterations for different numerical schemes. We observe that the convergence rate of kinetic damping with α 2 = 0.43 and kinetic damping without the control method are lower than those of other schemes. Figure 18a shows the response of the force at the top node. We find that the kinetic DRM converges immediately after the displacement reaches the goal (t = 0.1 s). Figure 18b shows the time history of kinetic energy for various schemes, which supports the results shown in Fig. 18a.

Synthesis of a monolithic constant-force compliant mechanism
To further illustrate the efficacy of the algorithm, we have applied the algorithm in a topology synthesis approach for deriving a constant-force compliant mechanism (CM). A constant-force mechanism is a device that can maintain invariant output force under a range of input displacements. For example, as shown in Fig. 19, a CM made of flexible material can exert a constant force at the control node when the position of the control node moves in a certain range of displacement. The topology synthesis method employs genetic algorithm (GA) optimization and VFIFE algorithm to form a functional model for the shape of the CM. The GA progresses by gradually altering discrete parameters for the model and adjusting fitness values from the preceding generations as it searches for an optimal solution corresponding to a predefined objective function. Each individual model as generated in GA is forwarded to VFIFE to compute for output force that corresponds to individual input displacement as listed in a predefined force-displacement relation. During the process of synthesis, the geometrical nonlinearity of the deformed structure is likely to cause numerical divergence where inconsistency of results can exist using the same time step size and density. Therefore, applying the adaptive VFIFE algorithm can avoid those drawbacks during the optimization of CM. In this example, the GA generates 40 individual configurations in each generation, with each individual being generated according to the rules of GA. For each individual, the output forces at 13 positions of the control node are computed via VFIFE. If the calculation in one of positions diverges, this configuration is viewed as an unsuccessful  one and the calculation stops to the next individual. It is intended to reduce the divergence of calculation during the synthesis since the occurrence of divergence may substantially slow down the overall topology synthesis procedure. Figure 20 shows the result of the optimal constant-force CM, where the structure contains 21 nodes and 20 elements. The control node only slides in the y-direction, and the output force was calculated at an interval of 2 mm of the control node input. The material of the structure is Nylon 6, with Young's modulus E = 3 × 10 9 N/m 2 , the area section A = 4 × 10 -6 m 2 , second moment of area of each element J z = 1.33 × 10 -12 m 4 and density ρ = 1150 kg/m 3 . Figure 21 shows the force-displacement curve for the optimized structure as well as the comparison between the target curve and the designed one, where the kinetic damping has been applied with the discrete control algorithm of α 2 = 0.49 and t = 1 × 10 -5 s. The maximum difference is about 2.02% based on |Target force − VFIFE force |/Target force between the range 8 and 24 mm. The accuracy is verified by traditional FEM using ABAQUS (red line). Figure 22 shows the comparison of the number of divergences in the calculation of the force-displacement of the structure with and without adaptive scheme in VFIFE. It can be observed that the algorithm with control scheme shows much fewer divergence cases after the eighth generation. However, the algorithm without control scheme has divergence cases in each generation of GA.

CONCLUSIONS
In this work, we investigated the stability and efficiency of solving nonlinear structural problems using VFIFE with kinetic DRM. We used concepts from the theory of discrete control to derive a method to improve the stability of the time integration scheme for the central difference method. It was demonstrated that by adjusting the coefficients of velocity and acceleration terms in the displacement equation, robustness in the stability of the numerical integration scheme can be achieved. We also showed that by using the schemes derived from kinetic DRM with control theory, the convergence rate can be improved. Numerical results revealed that the proposed procedures are sound and reliable for simulating statics for nonlinear structural problems. In the future, algorithms that can be used in different elements, such as membranes and solids, will be investigated.

ACKNOWLEDGE MENTS
The financial support of Ministry of Science and Technology of Taiwan (MOST 108-2221-E-002-132-) as well as National Taiwan University (NTU) to the first and corresponding authors is acknowledged. Such supports do not constitute an endorsement or implied endorsement by the supporting agencies of the views expressed in the article.

APPENDIX A . RODRIGUES' ROTATION FOR MUL A
If r 1 is a vector in space and s = [s x s y s z ] T is a unit vector about which r 1 rotates by an angle θ , then the rotated vector r 2 can be related to r 1 by r 2 = s cos θ + (s × r 1 ) sin θ + s(s · r 1 )(1 − cos θ ) (A1) or in matrix form r 2 = Rr 1 , (A2) where the rotation matrix R can be expressed as where I is an identity matrix and