Open Educational Resources

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 x_1, x_2, and x_3 such that the only requirement is that x_1<x_2<x_3. Another coordinate system \xi can be defined for the element such that the coordinates of the nodes 1, 2, and 3 are \xi_1=-1, \xi_2=0, and \xi_3=1, respectively. A mapping between the general coordinate system x and the coordinate system \xi can be defined (Figure 1) using the same shape functions that define the displacement in the coordinate system \xi as follows:

    \[x = N_1x_1 + N_2x_2 + N_3x_3\]

where the shape functions have the form:

    \[N_1 = \frac{-\xi (1-\xi)}{2}\]

    \[N_2 = (1-\xi)(1+\xi)\]

    \[N_3 = \frac{\xi(1+\xi)}{2}\]

The gradient of this mapping has the following form:

    \[\frac{dx}{d\xi} = \sum_{i=1}^3 \frac{dN_1}{d\xi} x_i\]

Since x_1<x_2<x_3, and since \xi\in[-1,1], then

    \[\frac{dx}{d\xi} > 0\]

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

    \[\frac{d\xi}{dx} = \frac{1}{\frac{dx}{d\xi}} = \frac{1}{\sum_{i=1}^3 \frac{dN_i}{d\xi} x_i}\]

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

    \[u = N_1u_1 + N_2u_2 + N_3u_3\]

The normal strain \varepsilon_{11} is defined in the x coordinate system as follows:

    \[\varepsilon_{11} = \frac{du}{dx} = \frac{du}{d\xi} \frac{d\xi}{dx} = \frac{d\xi}{dx} \hspace{1mm} < \hspace{2mm} \frac{dN_1}{d\xi} \hspace{2mm} \frac{dN_2}{d\xi} \hspace{2mm} \frac{dN_3}{d\xi} \hspace{2mm}> \hspace{1mm} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \end{pmatrix}\]

Figure 1. One Dimensional Quadratic Isoparametric Mapping

The virtual displacement u^* and the associated virtual strain \varepsilon_{11}^* can be chosen as follows:

    \[u^* = N_1u_1^* + N_2u_2^* + N_3u_3^*\]

    \[\varepsilon_{11}^* = \frac{d\xi}{dx} \hspace{1mm} < \hspace{2mm} u_1^* \hspace{2mm} u_2^* \hspace{2mm} u_3^* \hspace{2mm} > \hspace{1mm} \begin{pmatrix} \frac{dN_1}{d\xi} \\ \frac{dN_2}{d\xi} \\ \frac{dN_3}{d\xi} \\ \end{pmatrix}\]

Assuming p 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:

    \[IVW_e = A \int_{x_1}^{x_3} \sigma_{11} \varepsilon_{11}^* dx\]

    \[EVW_e = \int_{x_1}^{x_3} p u^* dx\]

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

    \[\int_{x_1}^{x_3} g(x) dx = \int_{-1}^1 g(x(\xi)) \frac{dx}{d\xi} d\xi\]

Using the constitutive law \sigma_{11}=E\varepsilon_{11}, the virtual work contributions of the element have the forms:

    \[IVW_e = \hspace{1mm} < \hspace{2mm} u_1^* \hspace{2mm} u_2^* \hspace{2mm} u_3^* \hspace{2mm} > \hspace{1mm} EA \int_{-1}^1 \Big(\frac{d\xi}{dx}\Big)^2 \frac{dx}{d\xi} \begin{pmatrix} \frac{dN_1}{d\xi} \\ \frac{dN_2}{d\xi} \\ \frac{dN_3}{d\xi} \\ \end{pmatrix} \hspace{1mm} < \hspace{2mm} \frac{dN_1}{d\xi} \hspace{2mm} \frac{dN_2}{d\xi} \hspace{2mm} \frac{dN_3}{d\xi} \hspace{2mm} > \hspace{1mm} d\xi \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \end{pmatrix}\]

    \[EVW_e = \hspace{1mm} < \hspace{2mm} u_1^* \hspace{2mm} u_2^* \hspace{2mm} u_3^* \hspace{2mm} > \hspace{1mm} \int_{-1}^1 p \begin{pmatrix} N_1 \\ N_2 \\ N_3 \\ \end{pmatrix} \frac{dx}{d\xi} d\xi\]

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

    \[K^e = EA \int_{-1}^1 \Big(\frac{d\xi}{dx}\Big)^2 \frac{dx}{d\xi} \begin{pmatrix} \frac{dN_1}{d\xi} \\ \frac{dN_2}{d\xi} \\ \frac{dN_3}{d\xi} \\ \end{pmatrix} \hspace{1mm} < \hspace{2mm} \frac{dN_1}{d\xi} \hspace{2mm} \frac{dN_2}{d\xi} \hspace{2mm} \frac{dN_3}{d\xi} \hspace{2mm}> d\xi\]

    \[f^e = \int_{-1}^1 p \begin{pmatrix} N_1 \\ N_2 \\ N_3 \\ \end{pmatrix} \frac{dx}{d\xi} d\xi\]

Example

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

Solution

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

    \[x = N_1x_1 + N_2x_2 + N_3x_3 = 4.5 + (5+0.5\xi)\xi\]

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

    \[K^e = EA \begin{pmatrix} 0.2653 & -0.3006 & 0.0353 \\ -0.3006 & 0.5465 & -0.2459 \\ 0.0353 & -0.2459 & 0.21067 \\ \end{pmatrix}\]

    \[f^e = p \begin{pmatrix} 1.333 \\ 6.667 \\ 2.000 \\ \end{pmatrix}\]

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 (x_1,y_1), (x_2,y_2), (x_3,y_3), and (x_4,y_4) be the two dimensional coordinates of nodes 1, 2, 3, and 4, respectively. The element coordinate system \xi and \eta 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:

    \[x = N_1x_1 + N_2x_2 + N_3x_3 + N_4x_4\]

    \[y = N_1y_1 + N_2y_2 + N_3y_3 + N_4y_4\]

where

    \[N_1 = \frac{(1-\xi)(1-\eta)}{4}\]

    \[N_2 = \frac{(1+\xi)(1-\eta)}{4}\]

    \[N_3 = \frac{(1+\xi)(1+\eta)}{4}\]

    \[N_4 = \frac{(1-\xi)(1+\eta)}{4}\]

Figure 2. Isoparametric Mapping in the Two Dimensional Case of a 4-node Bilinear Element

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

    \[J = \begin{pmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \\ \end{pmatrix}\]

    \[J^{-1} = \begin{pmatrix} \frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} \\ \frac{\partial \eta}{\partial x} & \frac{\partial \eta}{\partial y} \\ \end{pmatrix}\]

To demonstrate the role of these matrices, let v_{\xi,\eta}=\{\Delta \xi,\Delta \eta\} represent an infinitesimal vector in the element coordinate system \xi, and \eta. The coordinates of this vector in the general coordinate system x, and y are given by v_{x,y}=\{\Delta x,\Delta y\} such that:

    \[\begin{pmatrix} \Delta x \\ \Delta y \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \\ \end{pmatrix} \begin{pmatrix} \Delta \xi \\ \Delta \eta \\ \end{pmatrix}\]

i.e.,

    \[v_{x,y} = Jv_{\xi, \eta}\]

Similarly,

    \[v_{\xi, \eta} = J^{-1}v_{x,y}\]

The determinant of J can be used to relate areas in the coordinate system defined by x and y and the coordinate system defined by \xi and \eta (See the determinant section). From Figure 2, let d_1=\{\mathrm{d}\xi,0\} and d_2=\{0,\mathrm{d}\eta\} represent two vectors inscribing an area \mathrm{d}a=\mathrm{d}\xi\mathrm{d}\eta. Let dx_1=Jd_1, and dx_2=Jd_2 be the same two vectors in the coordinate system defined by x and y. Let the area measure in this coordinate system be \mathrm{d}A. The relationship between the two area measures is given by:

    \[dA = \det J da = \det Jd\xi d\eta\]

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

    \[\begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 & 0 \\ 0 & N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\ v_4 \\ \end{pmatrix}\]

    \[=Nu_e\]

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 \xi and \eta coordinate system have the following form:

    \[\begin{pmatrix} \frac{\partial u}{\partial \xi} \\ \frac{\partial u}{\partial \eta} \\ \frac{\partial v}{\partial \xi} \\ \frac{\partial v}{\partial \eta} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0 \\ \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} & 0 \\ 0 & \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 & \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\ v_4 \\ \end{pmatrix}\]

The gradients of the displacements functions in the x and y coordinate system can be evaluated as functions of the gradients of the displacements in the \xi and \eta coordinate system as follows:

    \[\begin{pmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} & 0 & 0 \\ \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} & 0 & 0 \\ 0 & 0 & \frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} \\ 0 & 0 & \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} \\ \end{pmatrix} \begin{pmatrix} \frac{\partial u}{\partial \xi} \\ \frac{\partial u}{\partial \eta} \\ \frac{\partial v}{\partial \xi} \\ \frac{\partial v}{\partial \eta} \\ \end{pmatrix}\]

    \[= \begin{pmatrix} \frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} & 0 & 0 \\ \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} & 0 & 0 \\ 0 & 0 & \frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} \\ 0 & 0 & \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} \\ \end{pmatrix} \begin{pmatrix} \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0 \\ \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} & 0 \\ 0 & \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 & \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\ v_4 \\ \end{pmatrix}\]

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

    \[\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \gamma_{12} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 \\ \end{pmatrix} \begin{pmatrix} \frac{\partial u}{\partial x} \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} \\ \frac{\partial v}{\partial y} \\ \end{pmatrix}\]

Setting:

    \[B = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 \\ \end{pmatrix} \begin{pmatrix} \frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} & 0 & 0 \\ \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} & 0 & 0 \\ 0 & 0 & \frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} \\ 0 & 0 & \frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y} \\ \end{pmatrix} \begin{pmatrix} \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0 \\ \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} & 0 \\ 0 & \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 & \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \\ \end{pmatrix}\]

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

    \[\varepsilon = Bu_e\]

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

    \[K^e = t \int_{-1}^1 \int_{-1}^1 B^T CB \det J d\xi d\eta\]

The element nodal forces can be evaluated as follows:

    \[f^e = t\int_{-1}^1 \int_{-1}^1 N^t \rho b \det J d\xi d\eta + \int_{element \hspace{1mm} surface} N^T t_n dS\]

where t is the thickness of the element, C is the plane strain or plane stress material constitutive relationship, \rho b is the body forces vector per unit volume, and t_n is the traction vector per unit area on the surface \mathrm{d}S. Notice that the surface element \mathrm{d}S is given in the x and y coordinate system and it should be transformed accordingly to the \xi and \eta 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 \rho b=\{\rho b_1, \rho b_2\} and a unit thickness. Assume plane stress conditions with E=1 units and \nu=0.

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 \det{J} 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

coordinates = {{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

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 \xi and \eta and the coordinate system x and y uses the quadratic interpolation functions (Figure 3). The interpolation functions in the \xi and \eta coordinate system are given by:

    \[N_1 = -\frac{(1-\eta)(1-\xi)(1+\xi+\eta)}{4}\]

    \[N_2 = -\frac{(1-\eta)(1+\xi)(1-\xi+\eta)}{4}\]

    \[N_3 = -\frac{(1+\eta)(1+\xi)(1-\xi-\eta)}{4}\]

    \[N_4 = -\frac{(1+\eta)(1-\xi)(1+\xi-\eta)}{4}\]

    \[N_5 = \frac{(1-\eta)(1-\xi)(1+\xi)}{2}\]

    \[N_6 = \frac{(1+\xi)(1-\eta)(1+\eta)}{2}\]

    \[N_7 = \frac{(1+\eta)(1-\xi)(1+\xi)}{2}\]

    \[N_8 = \frac{(1-\xi)(1-\eta)(1+\eta)}{2}\]

Figure 3. Isoparametric Mapping in the 8 Node Quadratic Two Dimensional Element

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

coordinates = {{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

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 \xi and \eta and the coordinate system of x and y 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 \xi and \eta coordinate system for the triangular elements are similar to the ones obtained previously.

Figure 4. Two Dimensional Isoparametric Mapping for the 3 node and 6 node Triangle Elements

Video:

Page Comments

  1. Luis Andres MOLLERICON TITIRICO says:
    April 3, 2024 at 1:59 pm

    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?

    Reply
    1. Samer Adeeb says:
      April 3, 2024 at 2:16 pm

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

      Reply
  2. JUMP ASSEMBLE says:
    June 9, 2024 at 11:20 pm

    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!

    Reply
    1. Samer Adeeb says:
      June 10, 2024 at 10:35 pm

      Thank you very much!

      Reply
  3. Nathan Robbins says:
    July 6, 2024 at 12:26 pm

    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?

    Reply

Leave a Reply

Your email address will not be published.