## One and Two Dimensional Isoparametric Elements and Gauss Integration: Isoparametric Elements

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 x 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

View Python Code

from sympy import * import sympy as sp xi, x1, x2, x3, EA, p = sp.symbols("xi x_1 x_2 x_3 EA p") N1 = -xi*(1-xi)/sp.sympify('2') N2 = (1-xi)*(1+xi) N3 = xi*(1+xi)/sp.sympify('2') display("Shape Functions:", N1,N2,N3) x = N1*x1+N2*x2+N3*x3 dxxi = x.diff(xi) dxix = 1/dxxi x = x.subs({x1:0,x2:4.5,x3:10}) display("Mapping Function: ", x) B = Matrix([N1.diff(xi), N2.diff(xi), N3.diff(xi)]) K1 = EA * dxix * B.transpose()*B Kb = integrate(K1, (xi,-1,1)) 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 = Matrix([a,b,c]).subs({x1:0,x2:9/sp.sympify('2'),x3:10}) display("Nodal forces: ",ff)

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

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

View Python Code

from sympy import * from mpmath import * import sympy as sp import scipy as sc import matplotlib.pyplot as plt from scipy.integrate import dblquad xi, eta, rb1, rb2 = sp.symbols("xi eta rb_1 rb_2") coordinates = Matrix([[0,0],[1,0.1],[1.2,1.2],[0.2,1]]) t = 1 display("Plane Stress") E = 1 nu = 0 Cc = E/(1+nu)/(1-nu)*Matrix([[1,nu,0],[nu,1,0],[0,0,(1-nu)/2]]) a = coordinates.row_insert(4, Matrix([[0,0]])) fig = plt.figure() ax = fig.add_subplot(111) cp = ax.plot(a[:,0], a[:,1]) ax.grid(True, which='both') plt.xlabel("x") plt.ylabel("y") plt.title("Shape") ax.axhline(y = 0, color = 'k') ax.axvline(x = 0, color = 'k') Shapefun = Matrix([[(1-xi)*(1-eta)],[(1+xi)*(1-eta)], [(1+xi)*(1+eta)],[(1-xi)*(1+eta)]])/4 Nn = Matrix([[0 for x in range(8)] for y in range(2)]) for i in range(4): Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i] x = sum([Shapefun[i]*coordinates[i,0] for i in range(4)]) y = sum([Shapefun[i]*coordinates[i,1] for i in range(4)]) J = Matrix([[x.diff(xi), x.diff(eta)],[y.diff(xi),y.diff(eta)]]) JinvT = J.inv().transpose() mat1 = Matrix([[1,0,0,0],[0,0,0,1],[0,1,1,0]]) mat2 = Matrix([[JinvT[0,0],JinvT[0,1],0,0], [JinvT[1,0],JinvT[1,1],0,0], [0,0,JinvT[0,0],JinvT[0,1]], [0,0,JinvT[1,0],JinvT[1,1]]]) mat3 = Matrix([[0 for x in range(8)] for y in range(4)]) for i in range(4): mat3[0,2*i] = mat3[2,2*i+1] = sp.diff(Shapefun[i],xi) mat3[1,2*i] = mat3[3,2*i+1] = sp.diff(Shapefun[i],eta) B = mat1*mat2*mat3 k1 = B.transpose()*Cc*B K = zeros(8) for i in range(8): for j in range(8): integrand=k1[i,j]*J.det() intlam=lambdify((eta,xi),integrand) #mp.dps is for integral accuracy (# of decimal points) mp.dps = 3 K[i,j] = quad(intlam,[-1,1],[-1,1]) # Scipy method, results in ZeroDivsionError # K = sc.zeros(k1.shape, dtype = float) # for (i,j), expr in sc.ndenumerate(k1): # tmp = sp.lambdify((xi, eta), expr*J.det(), 'math') # K[i,j] = dblquad(tmp, -1, 1, lambda xi: -1, lambda xi:1)[0] display(K) rb = Matrix([rb1,rb2]) fbeforeintegration = Nn.transpose()*rb*J.det() f3 = Matrix([integrate(fbeforeintegration[i],(xi,-1,1),(eta,-1,1))for i in range(8)]) display("Nodal Forces: ", f3)

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

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

View Python Code

from sympy import * from mpmath import * import sympy as sp import matplotlib.pyplot as plt xi, eta, rb1, rb2 = sp.symbols("xi eta rb_1 rb_2") coordinates = Matrix([[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 display("Plane Stress") E = 1 nu = 0.3 Cc = E/(1+nu)/(1-nu)*Matrix([[1,nu,0],[nu,1,0],[0,0,(1-nu)/2]]) a = Matrix([coordinates[0,:], coordinates[4,:], coordinates[1,:], coordinates[5,:], coordinates[2,:], coordinates[6,:], coordinates[3,:], coordinates[7,:],coordinates[0,:]]) fig = plt.figure() ax = fig.add_subplot(111) cp = ax.plot(a[:,0], a[:,1]) ax.grid(True, which='both') plt.xlabel("x") plt.ylabel("y") plt.title("Shape") ax.axhline(y = 0, color = 'k') ax.axvline(x = 0, color = 'k') Shapefun = Matrix([-(1-eta)*(1-xi)*(1+xi+eta)/4, -(1-eta)*(1+xi)*(1-xi+eta)/4, -(1+eta)*(1+xi)*(1-xi-eta)/4, -(1+eta)*(1-xi)*(1+xi-eta)/4, (1-eta)*(1-xi)*(1+xi)/2, (1+xi)*(1-eta)*(1+eta)/2, (1+eta)*(1-xi)*(1+xi)/2, (1-xi)*(1-eta)*(1+eta)/2]) Nn = Matrix([[0 for x in range(16)] for y in range(2)]) for i in range(8): Nn[0,2*i] = Nn[1,2*i+1] = Shapefun[i] x = sum([Shapefun[i]*coordinates[i,0] for i in range(8)]) y = sum([Shapefun[i]*coordinates[i,1] for i in range(8)]) J = Matrix([[x.diff(xi), x.diff(eta)],[y.diff(xi),y.diff(eta)]]) JinvT = J.inv().transpose() mat1 = Matrix([[1,0,0,0],[0,0,0,1],[0,1,1,0]]) mat2 = Matrix([[JinvT[0,0],JinvT[0,1],0,0], [JinvT[1,0],JinvT[1,1],0,0], [0,0,JinvT[0,0],JinvT[0,1]], [0,0,JinvT[1,0],JinvT[1,1]]]) mat3 = Matrix([[0 for x in range(16)] for y in range(4)]) for i in range(8): mat3[0,2*i] = mat3[2,2*i+1] = sp.diff(Shapefun[i],xi) mat3[1,2*i] = mat3[3,2*i+1] = sp.diff(Shapefun[i],eta) B = mat1*mat2*mat3 # This takes a long, long amount of time (The integration) k1 = N(B.transpose()*Cc*B,3) Ka = zeros(16) for i in range(16): for j in range(16): integrand=k1[i,j]*J.det() intlam=lambdify((eta,xi),integrand) #mp.dps is for integral accuracy (# of decimal points) mp.dps = 3 Ka[i,j] = quad(intlam,[-1,1],[-1,1]) print(k1) rb = Matrix([rb1,rb2]) fbeforeintegration = Nn.transpose()*rb*J.det() f3 = Matrix([integrate(fbeforeintegration[i],(xi,-1,1),(eta,-1,1))for i in range(16)]) display("Nodal Forces: ", f3)

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

When using this 4-node bilinear isoparametric formulation, we use the Gauss method to integrate and find the stiffness matrix (-1/sqrt(3), 1/sqrt(3)). So, does this mean that the values of strain and stress that we obtain after are values for this four gauss coordinates?

Correct. But one can output the stress and strain at any point within the element

Hi, just wanted to say that this blog post on isoparametric elements was really helpful in understanding the concept. The examples you provided were easy to follow and the math was explained in a way that made sense. Thanks for sharing your knowledge!

Thank you very much!

Thank you for the great resource. I am curious how non-planar shell elements would be formulated (I.e. 4 node quad in 3D space). Since the element in the natural coordinate system is planar, is the Jacobian still 2×2, and we would simply add additional columns to the last matrix of the [B] matrix for the new shape functions?

What about when we add rotational degrees of freedom to these nodes?