{'fraction_non_alphanumeric': 0.03432147784245769, 'fraction_numerical': 0.01605033816750227, 'mean_word_length': 4.534450651769087, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Isogeometric analysis of thin Reissner–Mindlin shells: locking phenomena and B-bar method

We propose a local type of B-bar formulation, addressing locking in degenerated Reissner–Mindlin shell formulation in the context of isogeometric analysis. Parasitic strain components are projected onto the physical space locally, i.e. at the element level, using a least-squares approach. The formulation allows the ﬂexible utilization of basis functions of different orders as the projection bases. The introduced formulation is much cheaper computationally than the classical ¯ B method. We show the numerical consistency of the scheme through numerical examples, moreover they show that the proposed formulation alleviates locking and yields good accuracy even for slenderness ratios of 10 5 , and has the ability to capture deformations of thin shells using relatively coarse meshes. In addition it can be opined that the proposed method is less sensitive to locking with irregular meshes.


Introduction
The conventional finite element method (FEM) employs polynomial basis functions to represent the geometry and the unknown fields. The commonly employed approximation functions are Lagrangian polynomials. However, these Lagrange polynomials are usually built upon a mesh structure which needs to be generated, from a CAD (computer-aided design) model provided for the domain of interest. This mesh generation leads to the loss of certain geometrical features: e.g. a circle becomes a polyhedral domain. Moreover, Lagrange polynomials lead to low order continuity between neighboring elements, which is inadequate in applications described by high order partial differential equations. The introduction of isogeometric analysis (IGA) [1] provides a general theoretical framework for the concept of "design-through-analysis" which has attracted considerable attention. The key idea of IGA consists in a direct link between CAD and the simulation, by utilizing the same functions to approximate the unknown field variables as those used to describe the geometry of the domain under consideration, similar to the idea proposed in [2]. Moreover, it also provides a systematic construction of high-order basis functions [3]. Note that, more recently, a generalization of the isogeometric concept was proposed, whereby the geometry continues to be described by NURBS functions as in the CAD, but the unknown field variables are allowed to live in different (spline) spaces. This leads to the concept of sub-and super-geometric analysis, also known as Geometry Independent Field approximaTion (GIFT), described within a boundary element framework in [4] and proposed in [5,6] and later refined in [7]. Related ideas, aiming at the construction of tailored spline spaces for local refinement were proposed recently in [8].
In the literature, the IGA has been applied to study the response of plate and shell structures, involving two main theories, viz., the Kirchhoff-Love theory and the Reissner-Mindlin theory. Thanks to the C 1 -continuity of the adopted NURBS basis functions, Kiendl et al. [9] developed an isogeometric shell element based on Kirchhoff-Love shell theory. The isogeometric Kirchhoff-Love shell element for large deformations was presented in [10]. The isogeometric Reissner-Mindlin shell element was implemented in [11], including linear elastic and nonlinear elasto-plastic constitutive behavior. The blended shell formulation was proposed to glue a Kirchhoff-Love shell with a Reissner-Mindlin shell in [12]. In addition, the isogeometric Reissner-Mindlin shell formulation derived from the continuum theory was presented in [13], in which the exact director vectors were used to improve accuracy. The solid shell element was developed in [14], in this formulation the NURBS basis functions were used to construct the mid-surface and a linear Lagrange basis function was used to interpolate the thickness field.
The Kirchhoff-Love type formulation is rotation-free, but it is only valid for thin structures. Moreover due to the absence of rotational degrees of freedom (DoFs), special techniques are required to deal with the rotational boundary conditions [9,15,16] and multi-patch connection [17]. Theoretically, the Reissner-Mindlin theory is valid for both thick and thin shell structures. However it is observed from the literature [11,18] that, when the shell kinematics is represented by Reissner-Mindlin theory, both FEM and IGA approaches suffer from locking for thin structures, especially for low order elements with coarse meshes. This has been attracting engineers and mathematicians to develop robust elements that could alleviate this pathology. Adam et.al. proposed a family of concise and effective selective reduced integration (SRI [19]) rules within the IGA framework for beams [20], plates/shells [21], and non-linear shells using T-splines [22]. Elguedj et al. [23] presentedB method andF method to handle nearly incompressible linear and non-linear problems. TheB method has been successfully applied to shear locking problems in curved beams [24], 2D solid shells [25], 3D solid shells [26] and later in nonlinear solid shell formulation [27]. Echter and Bischoff [18,28] proposed the discrete shear gap (DSG) method [29] to alleviate the shear locking phenomena effectively. Other approaches include twist Kirchhoff theory [30], virtual element method [31], collocation method [32,33], simple first order shear deformation theory [34], and single variable method [35,36]. The above approaches have been employed with Lagrangian elements or IGA framework with varying orders of success.
The approaches proposed by Bouclier, Elguedj and Combescure [25][26][27] focus on solid-shells in IGA and achieve good un-locking performance with theB method, thus it is promising to implement and test this method for degenerated Reissner-Mindlin shell elements in the context of IGA. This paper is the extension of our previous paper [37] on Timoshenko beams. In this research, in order to alleviate the locking phenomena, the locking strains are projected onto bi-dimensional low-order physical space by the least square method. The novel idea behind the formulation is to use multiple sets of basis functions to project the locking strains locally i.e. element-wise, instead of projecting globally i.e. all over the patch, thereby reducing the computational effort significantly. More importantly, the local projection allows one to use different orders of projecting basis functions, which could further improve the un-locking performance.
The outline of this paper is as follows: Sect. 2 gives an overview of the degenerated Reissner-Mindlin shell formulation. In Sect. 3, we present the novel approach, the localB method to alleviate the locking (both shear and membrane) problems encountered in thin structures whilst employing Reissner-Mindlin formulation. The accuracy, robustness and the convergence properties are demonstrated by some benchmark examples in Sect. 4, followed by concluding remarks in the last Sect. 5. Some implementation details are provided in the "Appendix". Figure 1 represents the mid-surface of a shell structure in the parametric and physical spaces. For simplicity, we consider a shell of constant thickness h, and assume the linear elastic material to be homogeneous and isotropic, which is governed by Young's modulus E and Poisson's ratio ν. The main differences between the Reissner-Mindlin shell theory and the Kirchhoff-Love shell theory lie in the assumptions on the deformation behavior of the section and further in the resulting independent kinematic quantities attached to the mid-surface. According to the Reissner-Mindlin theory, along the thickness direction a first order kinematic description is adopted for the transverse shear deformation. Assuming a Cartesian coordinate system, an arbitrary point P belongs to the 3D shell structure can be found by

Reissner-Mindlin shell model
where x stands for the shell mid-surface, modeled by a function of (ξ, η), as shown in Fig. 1. Let parameter ζ ∈ [−1, 1] denotes the shell thickness. On mid-surface x, the normal vector is calculated by and we denote Note that (t 1 , t 2 , n) form a local frame on the mid-surface, and the components of these vectors are indicated by subscript x, y, z respectively. Assuming small deformations, the displacement of the point P is recovered by in which u and θ are the displacement vector and the rotation vector at the mid-surface point projected from P. Then we have the linearized strain tensor

Isogeometric approach
Following the degenerated type formulation, we start our simulation from the mid-surface geometry of a  [40] and references therein. The geometry of the undeformed mid-surface is described by where x A define the locations of control points. The discrete displacement field on the mid-surface is interpolated as specifically u A = (u, v, w) T and θ A = (θ x , θ y , θ z ) T . Once the mid-surface is modeled according to Eq. (6) and (7), an arbitrary point P in the shell body can be traced by the following discrete forms and Regarding the exact normal vector field n, we construct an approximated field by collocated points, explained in "Appendix A". Upon employing the Galerkin framework and using the following discrete spaces for the displacement field the variational equation reads: find u ∈ S such that in which the bi-linear term is with D g denoting the global constitutive matrix. More details are provided in "Appendix B". When the displacements and the rotations are approximated by polynomials of the same order, the discretized framework experiences locking (shear and membrane) when the shell thickness becomes very small. The numerical procedure fails to satisfy the Kirchhoff limit as the shear strain does not vanish with the thickness approaching zero. One popular explanation of shear locking is that variables involved are not compatible [35], which is also known as field inconsistency. In the next section, we introduce theB method, trying to address the locking syndrome from the above point of view.

B-bar method for Reissner-Mindlin shells
In this section, we firstly illustrate the main idea behind thē B method. After introducing the classicalB method, in order to reduce the computational cost of the formulation, a local form of theB formulation is proposed, followed by a novel projecting scheme.

The idea behind the B-bar method
Assuming an unit length Timoshenko beam, which is discreted by one element with 2 nodes. The element length parameter x ∈ [0, 1]. By the help of shape functions the deflection field w and the rotation field θ are constructed as in which w i and θ i are nodal variables, i = 1, 2. For thin beams, the Timoshenko beam theory tends to degenerate into the Euler beam theory, and the shear strain (in Timoshenko beam theory) tends to vanish: and in discrete form denoted by its superscript h: The field inconsistency phenomenon occurs in the above equation, that is the order (wrt. x) of θ i is higher than the order of w i , which makes the shear strain difficult to vanish. TheB method modifies the field to be order-consistency. Firstly we build a pseudo shear strain as whereγĀ h is the corresponding pseudo DoF. Since the order of the initial shape functions N i is one, we define our one order lower shape function(s) asNĀ = 1 withĀ = 1. Next we perform a least square (L 2 ) projection and obtain a system of linear equation(s): and then Solve the above equation forγ 1 h thus the pseudo shear strain we previously built becomes Compare the pseudo shear strain in Eq. (23) with the original one in Eq. (17), it is noticed that the shape functions for w i remain unchanged, but the order of shape functions for θ i is reduced to be the same as w i . The pseudo shear strain is order consistency, it is believed that this strain could achieve good un-locking performance. From the strain matrix perspective, theB method modifies the original strain matrix which explains the name of this method.

Classical B-bar method in IGA
The novel idea behind theB method [23,24] is to use a projected strain instead of the original one to formulate the bi-linear term The projected strainε ε ε and the original strain ε ε ε are equivalent in the sense of the L 2 projection. Denoting the approximation space for displacement field by Q p,q , a common way is using the projection space of one order lower, i.e. Qp ,q = Q p−1,q−1 (see Fig. 2a). Lower order B-spline basis functionsNĀ are built from lower order knot vectors, with all the weights given as wĀ = 1,Ā = {1, . . . ,nm}. The L 2 projection process is performed on the physical domain as in which the discretized form of the projected strain is whereεĀ h means the projection of ε ε ε h ontoNĀ. Finally, we in which MĀB is the inner product matrix The Gram matrix M is inverseable since the involved bases are obviously linear-independent. In addition, the Hu-Washizu principle could be utilized to prove the variational consistency of theB method [23][24][25]41], and [24,26] proved that theB method is equivalent to a mixed formulation.

Local B-bar formulation for Reissner-Mindlin shells
The motivation of the localB concept comes from the observation that one needs to calculate the inverse of the matrix M in Eq. (29), which could be computationally expensive if the projection is applied patch-wise. Thus, from a practical point of view, it is highly recommended to project the strains locally [26,37], i.e. element-wise on Ω e : with This kind of treatment can be also found from the work of local least square method [38,39]. It is worthy pointing out that Eq. 31 implies the element strain matrix is modified as comparing with the original one B in Appendix Eq. B.8. Moreover, with the assumption that projecting the original strains onto a lower order space can help to reduce locking phenomena, it is expected that using the possibly lowest order for each element maximizes the un-locking performance. Therefore, instead of projecting the strains simply onto Q p−1,q−1 , different sets of projection spaces are adopted in this work. The projecting spaces need to be chosen carefully to avoid ill-condition or rand deficiency, readers interested in the corresponding mathematical theory is recommended to see [42], which is about volume locking (nearly-incompressible) problems. Here, the projection scheme in our previous work [37] is extended to bi-dimensional cases. Specifically, for space Q 2,2 , we adopt Q 1,1 for corner elements, Q 0,1 and Q 1,0 for boundary elements, and Q 0,0 for inner elements, as shown in Fig. 2b. This generalized strategy is similar with the SRI strategy in Fig. 2c, however the SRI strategy strongly depends on continuity conditions across elements [21]. The shape of the projecting functions are plotted in Fig. 3 for an intuitive understanding. The degenerated shell model is represented by its 2D midsurface, while the strain formula is expressed in 3D involving the thickness parameter ζ , in this case if the whole strain is projected, rank deficiency appears and the formulation yields inaccurate results. Following the similar approach as outlined in solid-shell element [25,26], in this work only the average strain through the thickness is projected: where MID(ε ε ε h e ) is defined to be the average strain through the thickness within a single element   The modified bi-linear form is written as The above calculation process is illustrated in Algorithm 1.
A primary comparison of the computational efforts is performed in Table 1, in which the meaning of the Step No. can be found in Algorithm 1. For classical IGA, the process is rather simple but it suffers from locking. In the classical B method, the calculation of the inverse of matrix M = (NĀ,NB) Ω requires large amounts of memory. Reviewing the literature of theB method since its appearance [43], one valuable contribution is the introduction of the localB concept [26,37], which saves lots of computational effort with acceptable accuracy loss. Compared with the globalB method, the localB method shows promising advantages in terms of computational efficiency. Based on the localB method, in the present research we choose different orders of basis functions as projection bases in order to achieve an improved un-locking performance.

Numerical examples
In this section, we demonstrate the performance of the proposed localB method for Reissner-Mindlin plates/shells by solving several standard benchmark problems. The proposed formulation is implemented within the open source C++ IGA framework Gismo 1 [44].   Fig. 2 Step2. Determine projection basesN according to current element position (see Fig. 2 Step3. Coefficient matrix MĀB = (NĀ,NB ) Ωe and its inverse Step4. Strain matrix B Step5. Calculate average strain through the thickness Step6. Least square projection Step7. Element stiffness matrix Step8 The following conventions are employed whilst discussing the results: -IGA 2: Classical IGA without any treatment of locking, approximation orders p = q = 2. -B-bar 2:B method, global projection, p = q = 2. -LB 2: LocalB method [26], the projection is applied element by element as shown in Eq. (31) and Fig. 2a, -SRI 2: Selective reduced integration scheme introduced in [21], p = q = 2 and we restrict C 1 continuity between elements. -GLB 2: LocalB method with present projection strategy as shown in Fig. 2b, p = q = 2.

Simply supported square plate
Consider a simply supported square plate with thickness h and edge length L = 1 m. Owing to symmetry, only one quarter of the plate, i.e, a = L/2 is modeled as shown in Fig. 4. The plate is assumed to be made up of homogeneous isotropic material with Young's modulus E = 200 GPa and Poisson's ratio ν = 0.3, and subjected to uniform pressure p.
The control points and the corresponding weights are given in Table 2.
The analytical out-of-plane displacement field for a thin plate with simply supported edges is given by where D = Eh 3 /(12(1 − ν 2 )). The transverse displacement is constant when the applied load is proportional to h 3 , thus the numerical solution is independent of the plate thickness. Due to the fact that w A indicates the maximum displacement of the plate, only the errors of w A are studied instead of the field error. The reference value is w A = −2.21804×10 −6 m. Figure 5 shows the normalized displacement as a function of mesh refinement. For thickness h = 10 −3 m, although the conventional IGA yields inaccurate results for coarse meshes, the results tend to improve upon refinement. However, IGA suffers from severe shear locking syndrome for h = 10 −5 m. In this case the results remain almost identical when increasing the number of elements, indicating that the elements are extremely locked. ForB method slight rank deficiency occurs for coarse meshes. GLB gets more accurate results than LB when h = 10 −5 m. It is inferred that the proposed formula- tion alleviates shear locking phenomena and yields accurate results even for extremely thin plates of slenderness ratio 10 5 . Figure 6 compares the convergence behavior by a severe locking case, the plate thickness is fixed as h = 10 −5 m. In this study, IGA of order 2 is locked for even more than 1000 elements, the convergence rate is nearly zero. Elements by IGA of order 3 start with smaller error, but suffer from locking until they are sufficiently refined. LB behaves the same as GLB for mesh 2 × 2. GLB achieves a good convergence rate in general, but at the well-refined steps the errors become unreasonably increased, which is also observed for IGA of order 3. The reason is that the converging results mismatch the reference value slightly, as pointed out in the partial enlarged Fig. 6b.
In Fig. 7 we analysis the computational time of the global/localB methods and IGA. Here we use a ×64 computer with Inter CPU i5-4590 (3.30 GHz), 8 GB RAM, then run the code only by a single thread. Overall theB methods add slightly computational effort to IGA, for instance about 2.64% in case of 64 × 64 elements, thus for a given accuracy theB methods are more cost-effective than IGA. For present method, the time is spend mostly on calculating the average strain (Step No. 5) and the least square projection (Step No. 6). For each element the dimension of the coefficient matrix M is one or two, therefore the computational time is not increased significantly. However in classicalB method, inverting M is more time consuming as the dimension of M is proportional to the number of elements (Fig. 7).

Scordelis-Lo roof
Here we consider the Scordelis-Lo roof problem to further demonstrate the capability of the proposed formulation in    Owing to symmetry only one quarter of the roof is modeled as shown in Fig. 8. The roof is modeled with control points (x, y, z) and weights w, given in Table 3 For the geometry considered here, R/h = 100 and L/h = 200, the structure experiences membrane locking as the transverse shear strain is negligible. The roof is dominated by membrane and bending deformations. The results of the normalized vertical displacement with mesh refinement is shown in Fig. 9. The results from the proposed formulation is compared with selective reduced integration technique [21]. It can be seen that accurate results are achieved by the present method. For coarse meshes theB methods, including classi-calB method and localB method as well as present method, lead to rank insufficient matrices, however for other problems slightly rank deficiency is also possible as recorded in the strain smoothing method by a single subcell [46]. In addition, the contour plot of the deflection w B by IGA and GLB are given in Fig. 10. It is can be noted that GLB captures the deformation quite well even for coarse meshes.

Pinched cylinder
The examples in Sects. 4.1 and 4.2 show that the proposed formulation yields accurate results when the structure experience either shear locking or membrane locking. We present here the pinched cylinder problem with the aim of demonstrating the robustness of the proposed formulation when the structure experiences both shear and membrane locking. Again, due to symmetry only one quarter of the cylinder is modeled as shown in Fig. 11. The corresponding control points and weights are given in Table 4  Results of w B for the Scordelis-Lo roof problem (Fig. 8) case to evaluate the performance of the method when the structure is dominated by bending behaviour. The convergence of the vertical displacement with mesh refinement is shown in Fig. 12 and it is evident that the proposed formulation yields more accurate results than the conventional bi-quadratic IGA. The contour plot of w C in Fig. 13 indicates Table 3 Control points and weights of the Scordelis-Lo roof problem (Fig. 8 . 10 Deflection field w × 10 2 for the Scordelis-Lo roof (Fig. 8) by IGA of order 2 (a, b), and by the proposed local B method (c, d)

Pinched hemisphere with hole
Consider a pinched hemisphere with 18 • hole subjected to equal and opposite concentrated forces applied at the four cardinal points. The hemisphere is modeled with Young's modulus E = 68.25 MPa, Poisson's ratio ν = 0.3 and concentrated force, P = 1 N. As before, owing to symmetry, only one quadrant of the hemisphere is modeled as shown in Fig. 14. The control points and weights employed to model the hemisphere are given in Table 5. This structure experiences severe membrane and shear locking, and the discretization experiences severe mesh distortion. The mesh distortion further enhances the locking phenomenon.
To evaluate the convergence properties, the horizontal displacement u ref D = 0.0940 m [21] is taken as the reference solution. The results form the proposed formulation are plot- Table 4 Control points and weights of the pinched cylinder problem (Fig. 11 12 Results of w C of the pinched cylinder (Fig. 11) ted in Fig. 15. As the elements are refined, GLB behaves similarly to SRI, but the results obtained by GLB is slightly accurate. Moreover, the field of displacement u x is plotted in Fig. 16, it seems that the proposed method allows to capture the correct structural response.

Square plate with irregular meshes
Throughout Sect. 4.1, it is pointed out that in the cases of plates with high slenderness ratio, e.g. 10 3 or 10 5 , the classical IGA elements suffer from locking, and presentB method has still achieved satisfactory results. However, the studies are performed in somewhat ideal conditions by using structured meshes. In order to explore the generality and applicability of the present method, we focus on the influence of two kinds of unstructured meshes, i.e. control mesh distortion and nonuniform knots.
The influence of the control mesh distortion on the performance of the proposed formulation is investigated by considering two different kinds of distortions: control mesh expansion and rotation, as shown in Fig. 17. Thickness h = 10 −5 m is considered here as the case when both locking and mesh distortion appear, and h = 10 −1 m is chosen as the reference case when only mesh distortion appears. The influence of mesh distortion on the normalized center displacement is pictured in Fig. 18. When locking and mesh distortion occur at the same time, accuracy of IGA decreases substantially, while GLB keeps a good degree of accuracy even for severe distortions. It is inferred that when locking happens, the results with classical IGA deteriorates with control mesh distortion, while the proposed formulation is less sensitive to control mesh distortion.
The influence of the non-uniformly distributed knots on the performance of the proposed formulation is investigated by considering a chosen knot ξ i moving from its left neighbor ξ i−1 to its right neighbor ξ i+1 . In the following example we build the plate geometry by and let ξ 5 move from ξ 4 to ξ 6 , η 5 move from η 4 to η 6 simultaneously, as illustrated in Fig. 19, while keeping the control mesh normally distributed. Figure 20 shows the results in terms of knots distribution, once again the classical IGA leads to poor results for unstructured knots distribution, but the proposedB formulation keeps good accuracy. It should be noted that, when ξ 5 = ξ 4 and η 5 = η 4 , or when ξ 5 = ξ 6 and η 5 = η 6 , the special cases of repeated knots (C 0 continuity) are also included. For presentB formulation, there is no need to consider inter knot multiplicity explicitly since the projection procedure is applied element-wise. Thus, lower order knot vectors can be used without inner repeated knots for projection, in spite of element continuity. A similar and detailed conclusion is drawn in [37].

Conclusion remarks
The localB method is adopted to unlock the degenerated Reissner-Mindlin shell elements within the framework of isogeometric analysis. The shell mid-surface and the unknown field are described with non-uniform rational Bsplines. The proposed method uses multiple sets of lower order B-spline basis functions as projection bases, by which the locking strains are modified in the sense of L 2 projection, in this way field-consistent strains are obtained. The salient features of the proposed localB method are: (a) has less computational effort than the classical, global B method; (b) yields better accuracy than IGA of order 2 especially in cases of coarse meshes or mesh distortions; (c) suppresses both shear and membrane locking commonly encountered when lower order elements are employed and suitable for both thick and thin models.
Future work includes extending the approach to large deformations, large deflections and large rotations as well as investigating the behaviour of the stabilization technique for enriched approximations such as those encountered in partition of unity methods [46][47][48].    . 16 Field of u × 10 2 of the pinched hemisphere (Fig. 14) by IGA of order 2 (a, b), and by present localB method (c, d)   Fig. 4) with nonuniformly knots distribution (Fig. 19). h stands for thickness surface point associated with the control point x A . Although no approximation error occurs, it is time-consuming to solve for the mid-surface projecting points of x A . In this paper the normal vectors at the Greville abscissae (ξ,η) are adopted to construct a normal vector field approximately [21], the field is re-written as R A (ξ, η)ñ A . Following Eq. 2, these normal vectorsñ A at the Greville abscissae points can be easily calculated. Thus we have and

Appendix B: Detailed formulation of degenerated Reissner-Mindlin shell element
In the following we introduce some necessary aspects to obtain the element stiffness matrix: (a) Using Voigt notation, the relation between the strain vector and the stress vector is expressed as where D g is the global constitutive matrix, and here D l is a given local constitutive matrix. To fulfill the plane stress state σ 33 = 0, the local constitutive matrix is given by in which κ is the shear correction factor and it is set to be 5/6 in this research. To fit the physical shape of the mid-surface, a transformation is employed for D l , specifically t 2x t 2y t 2y t 2z t 2z t 2x n x n x n y n y n z n z n x n y n y n z n z n x 2t 1x t 2x 2t 1y t 2y 2t 1z t 2z t 1x t 2y + t 2x t 1y t 1y t 2z + t 2y t 1z t 1z t 2x + t 2z t 1x 2t 2x n x 2t 2y n y 2t 2z n z t 2x n y + n x t 2y t 2y n z + n y t 2z t 2z n x + n z t 2x 2t 1x n x 2t 1y n y 2t 1z n z t 1x n y + n x t 1y t 1y n z + n y t 1z t 1z n x + n z t 1x (B.4) At each Gauss quadrature point, vectors in Eqs. 2 and 3 are calculated, then the transformation in Eq. B.2 is performed.
(b) Based on Eq. A.1, the Jacobi matrix is calculated by then the physical gradient of the shape functions are ⎡