In this tutorial, a quadrilateral 4 node isoparametric full integration element with coordinates as node-1: (0, 0), node-2: (5, 0), node-3: (5, 5), and node-4: (0, 5) is considered. Let's assume, Young's modulus, E = 20 GPa, and Poisson's ratio, ν = 0.2. We are interested here to compare the stiffness matrices for the plane strain condition (having unit thickness), which are determined by Mathematica as well as a finite element software, ABAQUS. Let's consider that the node-1 and -2 are fixed, and we are required to calculate the displacement of the other two nodes (3 and 4) under gravitational load. The material density is 2400 kg/m^{3}. We will also find out the stresses and strains at the integration points of the element.

The following Mathematica program is developed to determine the stiffness matrix, stress, strain and displacement manually (with Mathematica) of the nodal points of a quadrilateral 4 node isoparametric full integration element.

__Manual Approach (Mathematica):__

**(*Mathematica Codes to the Solve Problem*)**

**(*Definition of Nodal Coordinates from Problem 1*)**

**coordinates = {{0, 0}, {5,
0}, {5, 5}, {0, 5}};**

**(*Properties in Problem 1*)**

**t = 1;**

**Ee = 20*10^9;**

**Nu = 2/10;**

**Cc = Ee/((1 + Nu)*(1 - 2
Nu))*{{1 - Nu, Nu, 0}, {Nu, 1 - Nu, 0}, {0, **

** 0, 1/2 (1 - 2 Nu)}};**

**(*Defining the Shape Functions for Each Node*)**

**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 =
FullSimplify[Sum[Shapefun[[i]]*coordinates[[i, 1]], {i, 1, 4}]];**

**y =
FullSimplify[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[FullSimplify[mat1.mat2.mat3]];**

**Kbeforeintegration =
FullSimplify[Transpose[B].Cc.B];**

**xiset = {-1/Sqrt[3],
1/Sqrt[3]};**

**etaset = {-1/Sqrt[3],
1/Sqrt[3]};**

**Kk = t*FullSimplify[**

** Sum[(Kbeforeintegration*Det[J] /. {xi ->
xiset[[i]], **

** eta -> etaset[[j]]}), {i, 1, 2}, {j,
1, 2}]];**

**Kk = Round[Kk, 0.0001];**

**Print["The Stiffness
Matrix is,"]**

**ScientificForm[Kk //
MatrixForm, 4]**

**Print[""]**

**(*Determining
the Force Vector*)**

**rb = {0, -9.81*2400};**

**fbeforeintegration =
t*Transpose[Nn].rb*Det[J];**

**Ff = FullSimplify[**

** Sum[(fbeforeintegration /. {xi ->
xiset[[i]], **

** eta -> etaset[[j]]}), {i, 1, 2}, {j,
1, 2}]];**

**Print["The Force Vector
is,"]**

**FF = Transpose[**

** Join[{{"Fx_1 =", "Fy_1
=", "Fx_2 =", "Fy_2 =", "Fx_3 =", "Fy_3
=", **

** "Fx_4 =", "Fy_4 ="}},
{Ff}]];**

**ScientificForm[FF //
MatrixForm, 4]**

**Print[""] **

**(*Determining
the Displacement*)**

**ID = {5, 6, 7, 8};**

**Kkr = Kk[[ID, ID]];**

**Ffr = Ff[[ID]];**

**Uur = Inverse[Kkr].Ffr;**

**Uu = Table[0, {i, 1, 8}];**

**Uu[[ID]] = Uur;**

**Print["The Displacement
Vector is,"]**

**UU = Transpose[**

** Join[{{"u_1 =", "v_1 =",
"u_2 =", "v_2 =", "u_3 =", "v_3 =", **

** "u_4 =", "v_4 ="}},
{Uu}]];**

**ScientificForm[UU //
MatrixForm, 4]**

**Print[""] **

**(*Determining
the Stress and Strain*)**

**Eps = B.Uu;**

**Eps[[3]] = Eps[[3]]/2;**

**Sig = Cc.Eps;**

**EPS = ConstantArray[0, {4,
8}];**

**Do[EPS[[i, 1]] = i, {i, 1,
4}];**

**ii = 1;**

**Do[EPS[[ii, 1]] =
xiset[[i]];**

** EPS[[ii, 2]] = etaset[[j]];**

** EPS[[ii, {3, 4, 5}]] = Eps /. {xi ->
xiset[[i]], eta -> etaset[[j]]};**

** EPS[[ii, {6, 7, 8}]] = Sig /. {xi ->
xiset[[i]], eta -> etaset[[j]]};**

** ii++, {i, 1, 2}, {j, 1, 2}]**

**Print["Stress and Strain
Components are,"]**

**EPS =
Join[{{"x_i", "eta", "e_xx", "e_yy",
"e_xy", "s_xx", "s_yy", **

** "s_xy"}}, EPS];**

**EPS1 = EPS;**

**ScientificForm[EPS //
MatrixForm, 4]**

**Print[""]**

** **

__Program Outputs:__

Now, a model of a quadrilateral 4 node isoparametric full integration element is developed in Abaqus to determine the same stiffness matrix, stress, strain and displacement of the nodal points, which were obtained previously manually.

__Using Abaqus:__

In Abaqus, the geometry is defined by using the given nodal points in assignment. The material is elastic, where Young’s modulus is 20000000000 *pa*, and Poisson’s ratio is 0.2. In meshing, only a single element is considered with full integration option. The following two figures show the variations of displacements in magnitude.

The following figure shows the variations of stress in magnitude.

__The Output File for the Stresses:__

__The Output File for the Strains:__

__Comparison between Stiffness Matrices:__

**Mathematica Calculation**

By comparing all the results (stress, strain, and displacement) between manual calculation (Mathematica) and Abaqus model, we see that except the stress values, all other parameters are close. However, the stresses differ in the Abaqus model.