Weighted Integration Route to Stiffness Matrix of Quadrilaterals for Speed, Accuracy and Functionally Graded Material Application

Corresponding Author: P.V. Jeyakarthikeyan Department of Mechanical Engineering, SRM University, India Email: jeyapv@gmail.com Abstract: A weighted integration route with robust one-point integration (hourglass-controlled) is proposed as efficient, time saving alternative to Gauss quadrature for stiffness matrix of bilinear quadrilaterals. One-point rule relies on sampling at the center of the element to linearize the geometric transformation and average the material property over it. This enables, for a given element, explicit integration of stiffness matrix yielding a first approximation. For a second and better approximation, this procedure is applied independently to each of the four sub-squares of the mapped 2-square of the element and the matrices are assembled. A weighted addition of the two approximations produces a stiffness matrix as accurate as from 3-point Gauss-quadrature (G9P). Whereas, due to explicit integrations, obtaining stiffness matrix in this way demands less than a third of the time needed for 2-point Gauss-quadrature (G4P). On both counts (speed and accuracy) this approach outperforms Gaussquadrature. Sampling (material and geometry) at 5-points makes this element superior to G4P for Functionally Graded Material (FGM) applications. Bench mark examples by this approach are validated with Gauss quadrature and analytical solutions.


Introduction
Even to this day developing stiffness matrices for plane quadrilateral elements keeps enticing researchers (Jeyakarthikeyan et al., 2017). The various attempts include obtaining closed-form solutions as well as reducing the time needed for numerical integration. Simple exactly integrated numerical universal matrices like those for triangles and tetrahedrons (Subramanian and Bose, 1982;1983) also are possible only for parallelograms. But, unlike triangles, just parallelogram and rectangular elements alone are not sufficient to discretize arbitrary planar domains without including general convex quadrilaterals or triangles. This justifies developing stiffness matrices for quadrilaterals.
Numerical integrations have come to stay when developing stiffness matrices for finite elements. Attempts to integrate the functions exactly result many a time in complex logarithmic terms (serving mere academic interests) and their numerical evaluations are potentially time taking. Yet, there are instances when the shapes of elements lend themselves transformed linearly into standard shapes enabling exact integrations without such bothersome terms, saving computing time substantially (Jeyakarthikeyan et al., 2016;Zhou and Vecchio, 2006). Example: Triangulation of a planar domain with straight edged triangles (Jeyachandrabose and Kirkhope, 1984).
Finite element method itself being based on approximations, exact integrations do not mean much beyond a point. One prefers approximate numerical integrations to standardize computing procedures and save servicing costs. In FEM, when generating stiffness matrices, a lot of effort is put into identifying as small a number of sampling points as possible to minimize computational effort. From the very beginning Gauss quadrature is identified as the unchallenged option and adopted without question for all FEM applications (Irons and Ahmad, 1984). Commercial software routinely implement Gauss quadrature due to its ability to fit higher order polynomials with fewer sampling points improving accuracies and reducing computing time relative to conventional numerical methods (Newton-Cotes (Hoffman and Frankel, 2001;Zienkiewicz et al., 2005), Romberg etc., for example). It does not discriminate usual elements from special ones. Reduced Integration is quite tempting from time-saving angle but it could result in loss of stability and destroy the order of accuracy unless handled judiciously.
Numerical integrations are approximate and the degree of approximation depends not only on the order of integration but also on the distortion of the parent element in real plane. Efforts are taken to see that the shapes of finite elements are as regular as possible when generating finite element grids. Heavily distorted elements usually are avoided to reduce their harmful influence on the approximations. And so, there are routine warnings and elements with shapes not satisfying "close-to-regular" norms are flagged by commercial software. With these observations we proceed to describe a novel, fast and accurate alternate route to obtain stiffness matrices for plane quadrilateral elements. Hansbo (1998) presents an attractive one-point Gauss quadrature with stabilizing hourglass control for onedegree-of-freedom-per-node (1DPN) bilinear quadrilateral. An almost similar idea is enunciated by Mizukami (1986). The trick lies in integrating the basis functions leading to the stiffness matrix exactly with only Jacobian and material matrices sampled at the origin of the transform plane (mid-point rule). This enables explicitly integrated (approximate) stiffness matrix at the element level incorporating hourglassinstability control. The stiffness matrix is stable and computationally cheap but not accurate enough and so, several elements are needed to discretize a given domain and achieve a decent accuracy in the solution. We examine the element closely as under to see how this drawback can be overcome. Mizukami (1986) observes that this approximate stiffness matrix for a quadrilateral ( Fig. 1) is exact for a matching parallelogram. This parallelogram is readily identified in the following manner. Applying mid-point rule to Jacobian modifies the nonlinear geometric transformation into a linear one with a set of (leastsquare) substitute shape functions. This means that (ξ, η) to (x, y) transformation is 'explicitly' reversible. Thus, these linear substitute shape functions return the corner nodes of the 2-square in the transform plane as a set of corner nodes for a parallelogram instead of the original quadrilateral (Subramanian and Bose, 1982;1983) for a similar observation). It is verified further that the midpoints of the sides of this parallelogram and the quadrilateral are the same (Fig. 2). This means that the stiffness matrix obtained by this method is exact for the parallelogram and approximate for the quadrilateral.
It is easily verified that the following matrix 'QP' (quadrilateral to parallelogram) can transform (by post multiplication of the vector of corner coordinates of) any given quadrilateral into (corresponding corner coordinates of) a matching parallelogram in a least-square sense: Though, Mizukami (1986) and Hansbo (1998) (MH, for short), examine 1DPN quadrilaterals to develop stable one-point quadrature (for parallelograms!) as shown in Fig. (2), their idea is applicable equally to 2DPN quadrilaterals employed in structural mechanics. This forms the basis for the efficient computation of stiffness matrices for four-node plane quadrilaterals meriting FGM applications.
Computationally cheap and explicit 'one-point' (Reese, 2005;Jacquotte and Oden, 1984;Seleson, 2014;Cardoso and Yoon, 2005) MH stiffness matrix is yet approximate and for a prescribed level of accuracy many elements are needed to discretize effectively a given domain when, at the same time, fewer elements with relatively costly G4P (2-point Gauss quadrature) would be adequate. To strike a balance between these two extremes of speeding up computation and employment of fewer elements, the following is proposed. For a given quadrilateral the basic 'one-point' matrix is obtained explicitly as a first approximation. For a second approximation a cluster of four stable one-point quadrilaterals are carved out of the given quadrilateral and assembled into an (18×18) matrix (in principle) and reduced to (8×8) super element by enforcing the bilinear nature of the displacement description over the element (This procedure is not a static condensation, a procedure associated usually with standard super element formation, but is akin to obtaining transition elements from higher order elements: Zienkiewicz et al. (2005;Jeyachandrabose and Kirkhope, 1984)). These two approximations (with h-squared error) are weighted and extrapolated to give a vastly improved stiffness matrix that is as good as employing G9P. This new matrix is obtainable explicitly reducing cost of computation. Further, it is safely concluded that, since the arithmetic operations involve, independently stable matrices of the constituent elements, the resulting matrix also is stable. Thus, in a sense, every four quadrilaterals in the original discretization employing one-point integration can be replaced conveniently with larger, superior, computationally efficient and far more accurate, single, stable quadrilateral (Reese, 2005;Jacquotte and Oden, 1984;Seleson, 2014;Cardoso and Yoon, 2005;Zienkiewicz et al., 2005).

Explicit Stiffness Matrices for Plane Parallelogram Elements
Conventionally stiffness matrices for bilinear quadrilaterals are generated using G4P integration schemes that integrate cubic polynomials exactly. Any day, the terms to be integrated for obtaining stiffness matrix are never simple polynomials except when geometric transformation is linear. In quadrilaterals this is possible only for parallelograms (which include rectangles). For other shapes the terms to be integrated are of the form f(ξ, η) /g(ξ, η) where f and g are polynomials and Gauss quadrature turns approximate. Yet, (Mizukami, 1986, Hansbo, 1998 conclude independently that this approximation ('mid-point rule' for geometry and material) is stable and satisfactory and the stiffness matrix can be obtained explicitly.
For the entire element a first approximation for the stiffness matrix [K 1 ], is obtained using 'one-point integration' ('mid-point rule'). This matrix has an order of error h-squared. Following a procedure similar to that for triangles the stiffness matrix is conveniently generated in terms of numerical universal matrices for the matching parallelogram (as the mid-point rule for a quadrilateral renders the geometric transformation linear). Further, without loss of generality, origin is taken as a corner of this parallelogram in the (x, y) plane (Appendix: A).
If the MH quadrilateral is developed for a functionally graded material, it samples the values for the material constants only at its origin and the stiffness matrix is generated explicitly. In this form this is a basic (one-point) MH-Homogeneous Graded Material (HGM) element.
For the second approximation, the quadrilateral is divided into sub-quadrilaterals corresponding to the four sub-squares of its mapped 2-square in the (ξ, η) plane. One point integrations are carried out over each subsquare independently. In principle, triple matrix products are generated for all the sub-square elements and added. This gives a better estimate of the stiffness matrix for the given quadrilateral and the error is of the order of (h/2)-squared.
Thus there are two approximations for the stiffness matrix with different orders of error associated with them and a weighted addition produces a stiffness matrix possessing stability and greater accuracy. In a different context very recently the idea of obtaining stiffness matrices in this way is proposed by the authors for the first time. Actual calculations would follow the procedure outlined here to optimize computing effort.

Procedure
Considering the sub-element (1-5-9-8) of the quadrilateral (Fig. 1), one point integrated stiffness matrix K sub I takes the same form given earlier for the complete quadrilateral except that the coefficients of Gmatrix (Appendix: A) are calculated in terms of the coordinates of (the matching parallelogram for) the subelement and its material property associated with its midpoint. Also, the matrix is associated with displacements local to this element. To translate the contribution of this sub-element matrix to that of the global quadrilateral, a (8×8) transformation matrix C I must be operated upon K sub I in a triple-matrix-product form to yield C I T K sub I C I = K N I ; K N I is the individual contribution of this element to the stiffness matrix of the given quadrilateral (1-2-3-4). Similar contributions from the other sub-elements are added to obtain the complete stiffness matrix of the quadrilateral (1-2-3-4). The R matrix below in C I corresponds to the nodal order (1-5-9-8) of the first subelement in Fig. (1): Instead if the universal matrices A, B and C (Appendix A) are modified as A′= R T A R,B′ = R T B R and C′ = E T B R the stiffness matrix K N 1 for the subelement (1-5-9-8) is given by: In the above: Also therein, all are (4×4) numerically constant submatrices. To save time it is prudent to calculate in advance the (4×4) A′, B′ and C′ and generate K sub I , K sub II , K sub III and K sub IV . It should be noted that the R matrix and thus the A′, B′ and C′ based on the first sub-element nodal order (1-5-9-8), would remain the same also for the remaining three sub-elements of the quadrilateral only when their nodes are ordered cyclically as (2-6-9-5), (3-7-9-6) and (4-8-9-7) respectively. Accordingly the coefficients of G-matrix are calculated in terms of the co-ordinates local to the sub-element in the prescribed order. But this approach renders, for example, elements of K sub II match nodes in the order (2-3-4-1) of the given quadrilateral. So do elements of K sub III and K sub IV match with nodes (3-4-1-2) and (4-1-2-3) respectively. This requirement is taken care of easily using simple program logic when transferring the individual contributions of the sub-elements to the whole quadrilateral. Thus the total contribution from these four sub-elements to the stiffness matrix of the quadrilateral (1-2-3-4) gives [K 2 ].
One point sampling (only for geometric and material interpolation functions) being cheap computationally, easily saves time compared to Gauss Quadrature. The sampling further ensures that material gradation, if any, in the element is reflected in the final matrix, enabling Functionally Graded Material (FGM) element formulation effortlessly. In illustration, examples with FGM elements also are presented in this study.
Importantly four one-point quadrilaterals could be combined into one super quadrilateral for an added advantage. Once again it is not hard to see that reduction in time ensues from sampling at the Transform-Plane (TP) origin both for the first approximation and later for the second approximation (where each sub-element quadrilateral is sampled independently at its own TP origin) with only the constant terms of geometric function and material function contributing to the integral estimates. But then the field variable functions are integrated once and for all to obtain constant universal matrices A, B and C.
Though individually both estimates for the stiffness matrix are approximate, their weighted sum yields one with far greater accuracy. We realize that [K 1 ] and [K 2 ], are one-division and two-division approximations with errors of O(h 2 )  for the quadrilateral and may successfully be combined with weights. Taking β = (-1/3), the final K matrix is obtained as below: The weights are easily identified with those of wellknown Richardson's extrapolation. In passing, it is pointed out that Richardson's extrapolation by itself has a stabilizing effect Zlatev et al., 2010;2014a;2014b;Faragó et al., 2013). Despite the apparently elaborate exercise to find the stiffness matrix, the entire process still requires a third of the time needed for G4P. This is detailed after the numerical examples.

Numerical Examples
Example 1

A Cantilever Beam under Shear Load (Linear Elasticity)
A cantilever beam under shear load (often studied in the past by several researchers (Lee and Hobbs, 1998;Ramm et al., 2001;Chang-Chun et al., 1987;Yew et al., 1995;Della Croce and Venini, 2004)) is examined in the present context, Fig. (3a). This problem is analyzed with five distorted quadrilateral elements as in Fig. (3b). Table (1) compares the nodal displacements using elements obtained by the proposed method with those using MH quadrilaterals. Results of displacements using Gauss quadrature for these elements also are given. Clearly the results of the proposed method match closely with displacements using G9P, while MH quadrilateralbased results are far away. It is verified that even when each of the five elements is discretized further into four elements (i.e., a total of 20 elements as in Fig. 3c) the resulting nodal displacements show no perceptible improvement. The results using G9P (being a higher order quadrature), is taken as the benchmark to compare. The quadrilaterals of the proposed method and Gauss quadrature outperform MH quadrilaterals when their shapes are far from rectangle.
The convergence rate for this example is shown in Fig. (4), using the proposed method and conventional G4P. The convergence rate in the energy norm error indicator (€ es ) both by the proposed method and G4P is given in Fig. (5) (h = 1/5, h = 1/10, h = 1/20, h = 1/40 and h = 1/80). We plot log 10 (h) on the horizontal axis and log 10 (€ es ), on the vertical axis (Hansbo, 1998;Bui et al., 2014). The results from both the proposed method and G4P converge with reducing nodal ratio (h) almost at the same rate. It is clear that the proposed method produces results as good as Gauss quadrature. As for the accuracy of the method Table (1) compares the results for the cantilever bending example with both G4P and G9P favoring proposed method whose results are closer to the more accurate G9P. The energy norm is estimated with the following formula: Example 2

Functionally Graded Material (FGM) and Homogenous Graded Material (HGM) Quadrilateral Element Examples
Nowadays, one finds increasing use of functionally graded materials to alleviate effects of high stress concentration and withstand temperature gradients in structural applications (Kim and Paulino, 2002;Kubair and Bhanu-Chandar, 2008;Enab, 2014). It is shown here how the proposed element can be successfully employed in FGM related applications.
Consider a plate under plane stress conditions, assuming constant Poisson's ratio ν and Young's modulus E(x) varying in the following fashion (Kim and Paulino, 2002): The results from using FGM and HGM elements for different loading conditions are shown in Fig. 6a and 7a. The sampling points are, the center of any given element and centers of all sub-squares of the element in its mapped 2-square. In tune with the examples presented by Kim and Paulino (2002), Della Croce and Venini (2004), for comparison by the proposed method, the finite element mesh is created as 9×9, 4-node quadrilaterals (either MH-HGM or FGM). Figure 6b and 7b show FGM plate with material variation in the Cartesian x-direction for different boundary conditions (Kim and Paulino, 2002). In this section analytical and G4P solutions are compared with proposed solutions. For FGM plate, exponential material variations along x-direction are examined for the following loading conditions: • Bending • Tension The appropriate stress values obtained by G4P and proposed methods using FGM and MH-HGM elements are successfully compared with analytical solution. Figures (6b and 7b) show that the stresses obtained with homogeneous elements (MH) result in piecewise constant approximation as these elements sample just a single value each for material property. So it is clear that homogeneous elements estimate stress values only at their centroids. Also, the values for the stress jumps at the nodes for MH elements progressively increase with increasing values for x-coordinate for the loading conditions studied (Bending and Tension); Fig. 6b and  7b. These examples show that the proposed method for obtaining stiffness matrices of quadrilaterals is efficient and readily accommodate FGM applications (Enab, 2014;Mechee and Kadhim, 2016;Raffaella et al., 2016).    Figure 6c, 6d, 7c and 7d show the (node-averaged) stress contours by the proposed method, G4P and MH-HGM elements. The stress distributions using proposed method is the same as by G4P and (so only one figure is presented for both cases and) is seen to do better than MH-HGM elements. It is evident that the maximum and minimum stress values of MH-HGM elements are higher even though the aspect ratio of all the elements is unity. The scales for all the stress plots compared are kept the same to appreciate their relative merits.

Time Comparison: Flops and CPU Time by Proposed Method and G4P
Until recently, Floating point operation count (FLOPS) has been a standard to assess the execution speeds of algorithms. But, in modern computers, after introducing high speed computing engines with memory references and cache memory, FLOPS no longer is the norm in estimating execution speeds (Qian, 2015). We find that the FLOPS for obtaining stiffness matrix of a typical element by the proposed method is 10685. But, of these, only 2633 is required to find the stiffness matrices of all sub-elements for first and second approximations. The remaining 8052 is taken up for assembly operations. On the other hand G4P requires 4477 and G9P, 10062 FLOPS. These figures are quite misleading because when time calculations are performed on the same platform [MATLAB 13] the G4Pitself takes more than three times the time needed by the proposed method. This apparent anomaly when examined closely shows that the major amount of time taken for calculations by the proposed method is for sub-element stiffness matrix generation only. The assembly operations involving Equation (4) (with many zeros and ones) detailing 8052 flops practically takes no time. Table 2a presents the comparison of time between calculations performed separately for sub-matrices only against time for sub-element assembling only. It shows assembling is carried out in a flash.  This is easily understood as one-point sampling at (0,0) deals most of the time with integers whereas G4P works with float variables needing more time. The Table  2b, below summarizes this view and sets at rest the time efficiency question in favor of proposed method. The computation time needed for the four node quadrilateral element using G4P and proposed method are compared using the same platform (PC-Intel-i5, 2.3GHz, 8GB RAM) for four different numbers of elements (5×10 3 , 5×10 4 , 5×10 5 and 5×10 6 ) and is shown in Table 2b. Only the time required for the computation of the element stiffness matrix is accounted.

Conclusion
The foregoing analysis convincingly underlines that weighted integration route with robust one-point integration (hourglass-controlled) can wholly replace the traditional Gauss quadrature in terms of efficiency, speed and accuracy to obtain stiffness matrices of bilinear quadrilaterals. Sampling at the center of the quadrilateral (mid-point rule), by linearizing the geometric transformation and averaging the material property over an element, ensures obtaining stiffness matrix quickly and explicitly through constant universal matrices. (It is noted that the inferences elsewhere on the value of universal matrices in reducing computation time in connection with triangular elements (Jeyakarthikeyan et al., 2016) are equally valid here). Independently, the two approximations for the stiffness matrix for the quadrilateral are by themselves quite approximate but their weighted addition produces a stiffness matrix that is almost as accurate as G9P. This happens due to the application of the Richardson's extrapolation which plays on the error orders and improves the results many fold. Richardson's extrapolation also has the distinction (Jeyakarthikeyan et al., 2017) for stabilizing the combined results. Most satisfying is that the entire exercise needs less than a third of the time needed for G4P as confirmed by time calculations. In this instance one-point sampling and explicit integrations do the trick and save time.
Instead of sampling at four points as in G4P, here, in two steps, the sampling is done at five points creating relatively a better description of the material gradation over the element in the final matrix. Accordingly, FGM examples presented using this element confirm this. As an alternative, if needed, one could combine, stiffness matrix-wise, four one-point MH quadrilaterals into one large, far more accurate quadrilateral. All these go to show the effectiveness of the weighted integration path underlined in this study.
One does see that the method itself is not restricted to bilinear quadrilaterals and is readily extended to all members of the quadrilateral family. The idea is immediately applicable to three dimensional brick elements where 8 Gauss points can be replaced with 9 one point sampling and preliminary studies show that the resulting matrices are accurate and obtained in less than a third of the time needed for compatible Gauss quadrature scheme. This study is underway and will be presented in future. Strain: With D as the material matrix, strain energy equation for quadrilateral element is given by: On substituting the equation (A1) in the above equation we get:

In the above [G] = [P]T[D]
[P] is a (4×4) symmetric matrix of constants for the plane stress and plane strain problems and the values in it depend only on the material and geometry of the element. The material matrix D is a function of Young's modulus and Poisson's ratio for non-homogeneous elements (FGM). Young's modulus and Poisson's ratio are known functions of the spatial coordinates. When the material is considered homogeneous the matrix D is invariant with respect to the spatial position. In this study, if the element is non-homogenous, Lagrange polynomials are taken to define material variation and sampling is done at the centroid of the element.