## Finite Element Analysis: One and Two Dimensional Isoparametric Elements and Gauss Integration

In the previous section, the basic elements that are used for discretizing a plane domain were introduced. These elements can be used to mesh material bodies with regular geometries. For irregular geometries, quadrilateral and triangle elements will not necessarily follow the shapes of the basic elements described (a quadrilateral element will not necessarily be a square or a rectangle, and the triangle element will not necessarily have a right angle, as described before). Isoparametric elements enable meshing an irregular domain with triangular or quadratic elements that do not maintain orthogonality between the sides of the element. The purpose of the isoparametric formulation is to create shape functions that would ensure the compatibility of the displacement between neighboring elements while maintaining the requirements for shape functions mentioned in the previous section. When the elementsâ€™ sides do not necessarily have right angles, the isoparametric formulation is used to map those irregular elements into elements with regular shapes. The mapping functions used are the same shape functions that are used for the displacements, and hence, the name isoparametric (same parameters). After the mapping, integration is performed in the coordinate system of the element with the regular shape. The integration is performed numerically rather than using exact integrals. Numerical integration is preferred over exact integration for two reasons. The first is to reduce the computational resources required, and the second because the irregular geometries render radical (fraction) quantities that could be problematic for exact integrals. Many numerical integration techniques are available but the technique that is widely used for finite element analysis is the Gauss numerical integration technique.

### Isoparametric Elements

#### One Dimensional 3-Node Quadratic Isoparametric Element

the one dimensional quadratic element described previously was defined with the middle node equidistant from the end nodes. However, the isoparametric formulation allows defining a general one dimensional quadratic element with nodal coordinates , , and such that the only requirement is that . Another coordinate system can be defined for the element such that the coordinates of the nodes 1, 2, and 3 are , , and , respectively. A mapping between the general coordinate system and the coordinate system can be defined (Figure 1)using the same shape functions that define the displacement in the coordinate system as follows:

where the shape functions have the form:

The gradient of this mapping has the following form:

Since , and since , then

Therefore, the inverse of the gradient is well defined as:

If , , and are the horizontal displacements of the nodes 1, 2, and 3, respectively, then, the horizontal displacement function is defined in the coordinate system as follows:

The normal strain is defined in the coordinate system as follows:

The virtual displacement and the associated virtual strain can be chosen as follows:

Assuming as a distributed load per unit length of the one dimensional element, then, the contribution to the virtual work equation of this element are as follows:

To utilize the numerical Gauss integration scheme (as described later) Integration is performed in the coordinate system. In general, for a function we have:

Using the constitutive law , the virtual work contributions of the element have the forms:

Therefore, the stiffness matrix and the nodal forces vector for this element have the form:

##### Example

Let a nonlinear isoparametric one dimensional element be placed such that the nodal coordinates , , and are equal to 0, 4.5, and 10 units, respectively. Assume that a constant load per unit length is applied in the direction of . Find the stiffness matrix and the nodal forces for that element.

##### Solution

The mapping function between the two coordinate systems has the form:

Applying the above equations, the stiffness matrix and the nodal forces for the element have the form:

The following Mathematica code was utilized.

View Mathematica Code

Clear[x1, x2, x3, x, xi, p] N1 = -xi (1 - xi)/2; N2 = (1 - xi) (1 + xi); N3 = xi (1 + xi)/2; x = N1*x1 + N2*x2 + N3*x3; dxxi = FullSimplify[D[x, xi]]; dxix = 1/dxxi; x1 = 0; x2 = 4.5; x3 = 10; FullSimplify[x] B = FullSimplify[{{D[N1, xi], D[N2, xi], D[N3, xi]}}] K1 = FullSimplify[EA*dxix*Transpose[B].B]; Kb = Integrate[K1, {xi, -1, 1}]; Kb // MatrixForm a = Integrate[p*N1*dxxi, {xi, -1, 1}]; b = Integrate[p*N2*dxxi, {xi, -1., 1}]; c = Integrate[p*N3*dxxi, {xi, -1, 1}]; ff = {a, b, c}; Chop[ff] // MatrixForm

#### Two Dimensional 4-node Bilinear Isoparametric Element

Similar to the one dimensional case, an irregular quadrilateral can be defined with four nodes such that the only requirement is that the node numbering is counterclockwise and that the element area is positive. Let , , , and be the two dimensional coordinates of nodes 1, 2, 3, and 4, respectively. The element coordinate system and is chosen such that the coordinates of nodes 1, 2, 3, and 4 are: (-1,-1), (1,-1), (1,1), and (-1,1), respectively (Figure 2). The mapping between the two coordinate systems is performed using the same shape functions for the bilinear quadrilateral as follows:

where

The requirement that the element area is positive and that the node numbering is counterclockwise ensures that the Jacobian of the coordinate transformation is invertible at every point within the element. and its inverse have the forms:

To demonstrate the role of these matrices, let represent an infinitesimal vector in the element coordinate system , and . The coordinates of this vector in the general coordinate system , and are given by such that:

i.e.,

Similarly,

The determinant of can be used to relate areas in the coordinate system defined by and and the coordinate system defined by and (See the determinant section). From Figure 2, let and represent two vectors inscribing an area . Let , and be the same two vectors in the coordinate system defined by and . Let the area measure in this coordinate system be . The relationship between the two area measures is given by:

This formula can be used to transform integration over the and coordinate system to integration over the , and coordinate system. The displacement function in and coordinate system is given by:

To evaluate the strains, the gradients of the displacement functions in the two different coordinate systems need to be evaluated. The gradients of the displacement functions in the and coordinate system have the following form:

The gradients of the displacements functions in the and coordinate system can be evaluated as functions of the gradients of the displacements in the and coordinate system as follows:

The above relations allow the calculation of the two dimensional strains as follows:

Setting:

Therefore, the strains can be written in terms of the degrees of freedom as follows:

Then, the stiffness matrix of the element can be evaluated as follows:

The element nodal forces can be evaluated as follows:

where is the thickness of the element, is the plane strain or plane stress material constitutive relationship, is the body forces vector per unit volume, and is the traction vector per unit area on the surface . Notice that the surface element is given in the and coordinate system and it should be transformed accordingly to the and coordinate system.

Note that the results of this section can be extended in a straightforward manner to the three dimensional 8 node isoparametric brick element.

##### Example

Find the stiffness matrix and the nodal forces vector for a 4 node isoparametric plane element whose coordinates are (0,0), (1,0.1), (1.2,1.2), and (0.2,1) with a constant body forces vector and a unit thickness. Assume plane stress conditions with units and .

##### Solution

The following Mathematica code can be utilized to obtain the stiffness matrix and the nodal forces due to a constant body-force vector per unit volume for any 4 node isoparametric element as a function of the coordinates of the 4 nodes of the element. Study the code carefully. In particular, notice the following: If the nodes of the element are changed such that they are not defined in a counterclockwise fashion, then becomes negative and the stiffness matrix has negative diagonal entries (The stiffness matrix is not positive definite in that case!). Most finite element analysis software would give a warning message that such element is not adequately defined since the area of the element is negative. Also, notice that the stiffness matrix is obtained using the built-in numerical integration function NIntegrate in Mathematica. Try switching to the builtin exact integration function Integrate to notice the difference in the computational resources required. The numerical integration is a much faster technique.

View Mathematica Codecoordinates = {{0, 0}, {1, 0.1}, {1.2, 1.2}, {0.2, 1}}; t = 1; (*Plane Stress*) Ee = 1; nu = 0; Cc = Ee/(1 + nu)/(1 - nu)*{{1, nu, 0}, {nu, 1, 0}, {0, 0, (1 - nu)/2}}; a = Polygon[{coordinates[[1]], coordinates[[2]], coordinates[[3]], coordinates[[4]], coordinates[[1]]}]; Graphics[a] Shapefun = Table[0, {i, 1, 4}]; Shapefun[[1]] = (1 - xi) (1 - eta)/4; Shapefun[[2]] = (1 + xi) (1 - eta)/4; Shapefun[[3]] = (1 + xi) (1 + eta)/4; Shapefun[[4]] = (1 - xi) (1 + eta)/4; Nn = Table[0, {i, 1, 2}, {j, 1, 8}]; Do[Nn[[1, 2 i - 1]] = Nn[[2, 2 i]] = Shapefun[[i]], {i, 1, 4}]; x = Sum[Shapefun[[i]]*coordinates[[i, 1]], {i, 1, 4}]; y = Sum[Shapefun[[i]]*coordinates[[i, 2]], {i, 1, 4}]; J = Table[0, {i, 1, 2}, {j, 1, 2}]; J[[1, 1]] = D[x, xi]; J[[1, 2]] = D[x, eta]; J[[2, 1]] = D[y, xi]; J[[2, 2]] = D[y, eta]; JinvT = Transpose[Inverse[J]]; mat1 = {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 1, 1, 0}}; mat2 = {{JinvT[[1, 1]], JinvT[[1, 2]], 0, 0}, {JinvT[[2, 1]], JinvT[[2, 2]], 0, 0}, {0, 0, JinvT[[1, 1]], JinvT[[1, 2]]}, {0, 0,JinvT[[2, 1]], JinvT[[2, 2]]}}; mat3 = Table[0, {i, 1, 4}, {j, 1, 8}]; Do[mat3[[1, 2 i - 1]] = mat3[[3, 2 i]] = D[Shapefun[[i]], xi]; mat3[[2, 2 i - 1]] = mat3[[4, 2 i]] = D[Shapefun[[i]], eta], {i, 1, 4}]; B = Chop[Simplify[mat1.mat2.mat3]]; k1 = Simplify[Transpose[B].Cc.B]; K = t*NIntegrate[k1*Det[J], {xi, -1, 1}, {eta, -1, 1}]; rb = {rb1, rb2}; f3 = Integrate[Transpose[Nn].rb*Det[J], {xi, -1, 1}, {eta, -1, 1}]; f3 // MatrixForm K // MatrixForm

#### Two Dimensional 8-node Quadratic Isoparametric Element

The 8-node quadratic isoparametric plane element is better suited for curved surfaces since the mapping function between the coordinate system of and and the coordinate system and uses the quadratic interpolation functions (Figure 3). The interpolation functions in the and coordinate system are given by:

The development of the equations is similar to the previous section and is left to the reader. The only difference is the number of degrees of freedom of the element. The extension to the three dimensional case (the 20-node brick element) is also straightforward. The following is the Mathematica code that can be utilized to produce the stiffness matrix and the nodal forces due to a constant body forces vector for the 8-node isoparametric element as a function of the coordinates. Notice that the numerical integration built-in function in Mathematica has a higher accuracy than the traditional Gauss integration scheme (which will be presented later) but requires higher computational resources.

View Mathematica Codecoordinates = {{0, 0}, {1, 0.1}, {1.2, 1.2}, {0.2, 1}, {0.5, 0}, {1.1,0.65}, {0.7, 1.2}, {0.1, 0.5}}; t = 1; Ee = 1; nu = 0.3; (*Plane Stress*) Cc=Ee/(1+nu)/(1-nu)*{{1,nu,0},{nu,1,0},{0,0,(1-nu)/2}}; a=Polygon[{coordinates[[1]],coordinates[[5]],coordinates[[2]],coordinates[[6]],coordinates[[3]],coordinates[[7]],coordinates[[4]],coordinates[[8]],coordinates[[1]]}]; Graphics[a] Shapefun=Table[0,{i,1,8}]; Shapefun[[1]]=-(1-eta)(1-xi)(1+xi+eta)/4; Shapefun[[2]]=-(1-eta)(1+xi)(1-xi+eta)/4; Shapefun[[3]]=-(1+eta)(1+xi)(1-xi-eta)/4; Shapefun[[4]]=-(1+eta)(1-xi)(1+xi-eta)/4; Shapefun[[5]]=(1-eta)(1-xi)(1+xi)/2; Shapefun[[6]]=(1+xi)(1-eta)(1+eta)/2; Shapefun[[7]]=(1+eta)(1-xi)(1+xi)/2; Shapefun[[8]]=(1-xi)(1-eta)(1+eta)/2; Nn=Table[0,{i,1,2},{j,1,16}]; Do[Nn[[1,2i-1]]=Nn[[2,2i]]=Shapefun[[i]],{i,1,8}]; x=Sum[Shapefun[[i]]*coordinates[[i,1]],{i,1,8}]; y=Sum[Shapefun[[i]]*coordinates[[i,2]],{i,1,8}]; J=Table[0,{i,1,2},{j,1,2}]; J[[1,1]]=D[x,xi]; J[[1,2]]=D[x,eta]; J[[2,1]]=D[y,xi]; J[[2,2]]=D[y,eta]; JinvT=Transpose[Inverse[J]]; mat1={{1,0,0,0},{0,0,0,1},{0,1,1,0}}; mat2={{JinvT[[1,1]],JinvT[[1,2]],0,0},{JinvT[[2,1]],JinvT[[2,2]],0,0},{0,0,JinvT[[1,1]],JinvT[[1,2]]},{0,0,JinvT[[2,1]],JinvT[[2,2]]}}; mat3=Table[0,{i,1,4},{j,1,16}]; Do[mat3[[1,2i-1]]=mat3[[3,2i]]=D[Shapefun[[i]],xi];mat3[[2,2i-1]]=mat3[[4,2i]]=D[Shapefun[[i]],eta],{i,1,8}]; B=Chop[FullSimplify[mat1.mat2.mat3]]; k1=FullSimplify[Transpose[B].Cc.B]; Ka=t*NIntegrate[k1*Det[J],{xi,-1,1},{eta,-1,1},AccuracyGoal->8]; Ka//MatrixForm rb={rb1,rb2}; f3 = Chop[NIntegrate[Transpose[Nn]*Det[J], {xi, -1, 1}, {eta, -1, 1}, AccuracyGoal -> 8].rb]; f3//MatrixForm

#### Two Dimensional 3 node and 6 node Triangle Isoparametric Elements

The development of the isoparametric triangle elements is also similar to the previous section. The isoparametric mapping between the coordinate system of and and the coordinate system of and is shown in Figure 4 for the two types of triangular elements studied (Linear and quadratic). The development of the equations follows the previous section. The shape functions in the and coordinate system for the triangular elements are similar to the ones obtained previously.

### Gauss Integration

The Gauss integration scheme is a very efficient method to perform numerical integration over regular domains. In fact, if the function to be integrated is a polynomial of an appropriate degree, then the Gauss integration scheme produces exact results. The Gauss integration scheme has been implemented in almost every finite element analysis software due to its simplicity and computational efficiency. This section outlines the basic principles behind the Gauss integration scheme, along with its application in the finite element analysis method. In particular, this section will discuss the development of the reduced-integration elements and the benefits and drawbacks of utilizing them in finite element analysis models.

#### Motivation

The isoparametric elements developed above imply that all the integrated functions are defined such that . In general, numerical integration of such a function has the form:

For the Gauss integration method, is called an integration point and is called the associated weight. Now, if is always affine, i.e., then:

So, for functions that are very close to being affine, a numerical integration scheme with 1 integration point that is with an associated weight of 2 can be employed.

#### Gauss Integration over One Dimensional Domains

Let . Then, the definite integral of over the domain can be approximated by the following formula:

where and is a weight factor associated with each . are termed the integration points. The number of integration points required to obtain sufficient accuracy depend on the shape or form of the function . If is a linear function, then one integration point is enough to exactly integrate . If is a cubic polynomial function, then two integration points are sufficient to obtain the exact solution. Figure 5 shows a graphical representation of the Gauss numerical integration scheme for integrating polynomials up to the fifth degree. Following is a table of the number of integration points required to obtain the exact solution for a polynomial function with different polynomial degrees:

##### Example

Verify that the exact integral of a general cubic polynomial on the interval can be obtained by the two integration points and with the two weight factors .

##### Solution

Let . The exact integral is equal to:

On the other hand, Gauss integration with two integration points is equal to:

Equating the Gauss integration with the exact solution yields the following four equations:

These four equations are obtained by noticing that the exact solution and the Gauss integration solution should be equal given any values for the coefficients , , , and . The solution to the above four equations yields , , . The Mathematica code to solve the four equations is:

View Mathematica CodeSolve[{w1 + w2 == 2, w1*x1^3 + w2*x2^3 == 0, w1*x1^2 + w2*x2^2 == 2/3,w1*x1 + w2*x2 == 0}, {w1, w2, x1, x2}]

Note that this example can be extended to show that Gauss integration in one dimension with points can exactly integrate a polynomial of order .

#### Gauss Integration Over Two Dimensional Domains

Let . Then, the integral of over the domain can be approximated by the following formula:

where is an integration point while and are weight factors associated with and respectively. It can be easily seen that this is a direct consequence of the one dimensional Gauss integration formula as follows:

#### Gauss Integration for Two Dimensional Quadrilateral Isoparametric Elements

As described for the two dimensional linear quadrilateral element, the strains of the 4 node elements are linear, and thus, the entries in the matrix are quadratic functions of position. Thus, a numerical integration scheme is required to produce accurate integrals for the isoparametric 4-node elements (Figure 6). On the other hand, the strains in the two dimensional quadratic quadrilateral are quadratic functions of position, and thus, the entries in the matrix are fourth-order polynomials of position. Thus, a numerical integration scheme is required to produce accurate integrals for the 8 node isoparametric element (Figure 6). It should be noted that this can extended the three dimensional case in a straightforward manner. The 8 node brick element has 8 integration points, while the 20 node brick element has 27 integration points.

#### Reduced Gauss Integration (Under Integration) for Two Dimensional Quadrilateral Isoparametric Elements

Elements that are integrated as per the previous section are termed “full integration” elements. For example, to integrate the stiffness matrix of the 4-node plane isoparametric quadrilateral element, each entry will be evaluated at the four locations () of the integration points in the isoparametric element. Thus, the number of computations (excluding the addition) for the integration of the element are computations. The number of computations for the stiffness matrix of the 8-node plane isoparametric quadrilateral element is 2,304 computations since it uses a Gauss integration points scheme! To save computational resources, a “reduced integration” option for each of those elements is available in finite element analysis software. The reduced integration version of the 4-node plane isoparametric quadrilateral element uses only one integration point in the center of the element (Figure 7), and the number of computations is reduced from 256 to 64 computations. However, this is accompanied by a reduction in accuracy. Figure 7 shows two nonlinear functions to be integrated using one integration point vs. two integration points. In the first case, the function has a small variation along the domain and the integration seems to be accurate. However, in the second case, the function has a large variation and is equal to zero in the centre of the domain. The reduced integration will predict zero stiffness!

For the 8-node plane isoparametric quadrilateral, a Gauss integration points scheme is used instead of the Gauss integration points scheme, which reduces the number of computations from 2,304 to 1,024 computations (Figure 8). The benefit of using reduced integration is not only in saving computational resources but also in balancing out the over-stiffness introduced by assuming a certain deformation field within an element. In fact, reduced integration can lead to better performance in most cases. However, while reduced integration can dramatically reduce the computational time, it can lead to a structure that is too flexible. In fact, the displacements are over predicted when a structure is discretized using reduced integration elements. The reason, as will be shown in the following section, is because some information is lost during sampling of the stiffness at the fewer integration points of the reduced integration elements.

##### Spurious Modes Associated with Reduced Integration Elements

Spurious modes, or sometimes called zero energy modes, are combinations of displacements that in general produce strains; however, such strains cannot be detected at the integration or sampling points. So, the internal energy calculated at the integration points is equal to zero for such mode. In other words, a very small amount of load would be enough to produce a very large displacement. If a spurious mode exists for the whole structure, it becomes unstable. Usually spurious modes exist in under-integrated element (reduced integration). This is because the reduced number of integration points cannot detect a spurious strain mode. The two very common spurious modes are those for the bilinear quadrilateral (4-node plane element) and the quadratic displacement quadrilateral (8-node plane element). In general, any bending mode is a spurious mode in a 4-node reduced-integration linear element. The strain modes shown cannot be detected at the integration point; the strain at the center of the element (at the location of the single integration point) is equal to zero (Figure 9). On the other hand, the spurious mode of an 8 node reduced integration quadratic element is called “the hour glass” mode (Figure 10) and is produced with the following configuration of nodal displacements:

The spurious modes will be illustrated using the following problem.

##### Example

Find the internal energy associated with the shown displacement for the 4-node plane stress quadrilateral element shown below, assuming a bilinear displacement interpolation function.

Use exact integration, full integration ( integration points), and reduced integration (1 integration point).

##### Solution

The shape functions of the shown element are:

The displacement functions and in the directions and , respectively, are:

Therefore, the strains have the form:

The stress strain relationship for plane stress is given by:

The stiffness matrix can be evaluated using the following integral where is the element thickness:

Notice that the system has only one degree of freedom (the unknown displacement variable ), and therefore, the stiffness matrix has the dimensions .

**Using Exact Integration:**

The stiffness matrix evaluated using exact integration is:

The strain energy associated with this displacement mode is equal to:

**Using Full Integration:**

The stiffness matrix evaluated using full integration is:

where and which gives the same results as the exact integration.

** Using Reduced Integration:**

The stiffness matrix evaluated using reduced integration is:

where and . The strain energy associated with this displacement mode is equal to zero. Thus, the reduced integration scheme predicts that the displacement field shown does not require any external forces applied to the element. Thus, this element is unstable if loads produce the displacement mode shown, i.e., bending! The following is the Mathematica code used for the above calculations:

View Mathematica CodeN1=(1-xi)(1-eta)/4; N2=(1+xi)(1-eta)/4; N3=(1+xi)(1+eta)/4; N4=(1-xi)(1+eta)/4; B={D[N1-N2+N3-N4,xi],0,D[N1-N2+N3-N4,eta]}; Cc=Ee/(1+nu)/(1-nu)*{{1,nu,0},{nu,1,0},{0,0,(1-nu)/2}}; Kbeforeintegration=FullSimplify[B.Cc.B] (*Using Exact Integration*) Kexact=Integrate[Kbeforeintegration,{eta,-1,1},{xi,-1,1}] Energyexact=us*Kexact*us/2 (*Using Full Integration*) etaset={1/Sqrt[3],-1/Sqrt[3]}; xiset={1/Sqrt[3],-1/Sqrt[3]}; Kfull=FullSimplify[Sum[(Kbeforeintegration/.{xi->xiset[[i]],eta- >etaset[[j]]}),{i,1,2},{j,1,2}]] Energyfull=us*Kfull*us/2 (*Using reduced Integration*) etaset=xiset={0}; Kreduced=Sum[4(Kbeforeintegration/.{xi->xiset[[i]],eta- >etaset[[j]]}),{i,1,1},{j,1,1}] Energyreduced=us*Kreduced*us/2

#### Gauss Integration for Two Dimensional Triangle Isoparametric Elements

The strains of the 3 node triangle element are constant, and thus, the entries in the matrix are constant. Thus, a 1 point numerical integration scheme is required to produce accurate integrals for the isoparametric 3-node elements (Figure 11). For the 6 node triangle elements, a 3 points numerical integration scheme is required to produce accurate integrals (Figure 11). It should be noted that the shown locations for the integration points and their number can be extended to the three-dimensional case in a straightforward manner.

### Final Notes

The use of isoparametric elements and numerical integration dramatically increases the robustness of the finite element analysis method. The techniques described in the previous sections allow meshing of irregular domains with compatible elements, i.e., elements that ensure the continuity of the displacement field across boundaries (How?). However, to ensure the continuity of the displacement field, the elements of the mesh should be of the same type; otherwise, the user has to apply mesh constraints to ensure the compatibility of the displacement between elements. For example, a 4-node bilinear displacement element when connected to an 8-node quadratic displacement element will have incompatible displacements at the side of the connection. This can be overcome by constraining the mid-side node of the quadratic element to always lie on the straight side of the bilinear element. Otherwise, the results at the location of discontinuous displacement fields would be erroneous.

The stiffness matrix and the nodal loads of the structure can be calculated using the Gauss numerical integration scheme, which speeds up the computational time required. However, the price of using the Gauss numerical integration is the loss of accuracy. Gauss numerical integration produces accurate results only when the integrals are exact polynomials. The integrals are exact polynomials only when the elements of the mesh are exactly aligned with the isoparametric mapping into the coordinates of and . When the shape of the elements in the mesh deviates from the shape of the element in the coordinate system of and , the functions to be integrated are no longer exact polynomials but rather radicals; thus, the numerical integration scheme is no longer accurate. Most finite element analysis software give warnings when elements are highly distorted since this will lead to inaccurate computations of the stiffness matrix during the numerical integration process. In addition, for nonlinear problems when the stiffness changes with deformation, it is possible that the element will reach a highly distorted shape that the numerical integration scheme is no longer valid.

Hi,

Please could you explain what is on the y axis of Figure 7? Thank you!

It is the function to be integrated. In this case, you can think of it as the strain energy ($E epsilon^2/2$).