This C++ Code calculates the stiffness matrix for a given problem. It has two options: you can choose either Q4 element or Q8 element. It uses the 3 Gauss Points. The x and y coordinates for the 2D rectangular element should be inputted manually. The calculation we be given in an "output.data" file. Here is brief description of Q4 and Q8
Q4 element:
Q4 element is a quadrilateral element that has four nodes. Its displacement field has four terms. It is called Bilinear Quadratic.
Q8 element:
Q8 element is a quadrilateral element that has eight nodes. Its displacement field has eight terms. It is called Quadratic Quadrilateral.
User is considered to have a basic knowledge of Finite Element Method.
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <fstream>
#include <time.h>
#include <conio.h>
#include <string>
using namespace std;
void main()
{
string output;
ofstream out;
output = "output.data";
out.open(output.c_str());
int i, j, a,b,c,f,g,k,h, l,m,p, q, nelnodes;
double dxdxi[2][2], IJ[2][2], coord[8][2], dNdxi[2][8], Jdeterm, dNdx[2][8];
double xi[2][9],Kel[16][16], D[3][3], B[3][16], BT[16][3], BTD[16][3],w[3][3]
,w1D[3], BTDB[16][16];
const double E = 30000000., nu = 0.25;
const int ncoord =2; const int npoints = 9;
cout<<"\n===============================================";
cout<<"\n ME862 HW06 Habib Ahmadi";
cout<<"\n Element Stiffness Matrix, Q4 and Q8 with 3GP";
cout<<"\n===============================================";
cout<<endl;
//input for the number of element node element node and global coordinate
cout<< "Please enter number of element nodes"<<endl;
cin>>nelnodes;
cout<<"Enter coordinate values, all local X value first, then Y values"<<endl;
if(nelnodes ==4)
{
for(b=0;b<ncoord; b++)
{
for(a=0;a<nelnodes;a++)
{
cin>>coord[a][b];
}
}
}
else if (nelnodes == 8)
{
for(b=0;b<ncoord; b++)
{
for(a=0;a<nelnodes;a++)
{
cin>>coord[a][b];
}
}
}
for(p=0;p<2*nelnodes;p++)
{
for(q=0;q<2*nelnodes;q++)
{
Kel[p][q]= 0.;
}
}
for(i=0;i<3;i++)
{
for(j=0;j<3;j++)
{
//====================== No. integration points =============================
//
// Defines the number of integration points:be used for
// each element type
// if (nelnodes == 4)
// n = 4;}
//else if (nelnodes == 8)
//{ n = 9;}
//====================== INTEGRATION POINTS ==================================
// Defines positions of integration points
// 2D Rectangular elements
// ( nelnodes==4 || nelnodes==8 )
xi[0][0]= -0.7745966692;
xi[1][0] = xi[0][0];
xi[0][1] = 0.0;
xi[1][1] = xi[0][0];
xi[0][2] = -xi[0][0];
xi[1][2] = xi[0][0];
xi[0][3] = xi[0][0];
xi[1][3] = 0.0;
xi[0][4] = 0.0;
xi[1][4] = 0.0;
xi[0][5] = -xi[0][0];
xi[1][5] = 0.0;
xi[0][6] = xi[0][0];
xi[1][6] = -xi[0][0];
xi[0][7] = 0.;
xi[1][7] = -xi[0][0];
xi[0][8] = -xi[0][0];
xi[1][8] = -xi[0][0];
//================= INTEGRATION WEIGHTS ==================================
//
// Defines integration weights w_i, done for 3 points
// 2D Rectangular elements
// ( nelnodes==4 || nelnodes==8 )
w1D[0] = 0.555555555;
w1D[1]= 0.888888888;
w1D[2] = 0.55555555555;
w[j][i] = w1D[j]*w1D[i];
//================= SHAPE FUNCTION DERIVATIVES ======================
// 2D Rectangular elements
if ( nelnodes == 4 )
{
dNdxi[0][0] = -0.25*(1.-xi[1][3*i+j]);
dNdxi[1][0] = -0.25*(1.-xi[0][3*i+j]);
dNdxi[0][1] = 0.25*(1.-xi[1][3*i+j]);
dNdxi[1][1] = -0.25*(1.+xi[0][3*i+j]);
dNdxi[0][2] = 0.25*(1.+xi[1][3*i+j]);
dNdxi[1][2] = 0.25*(1.+xi[0][3*i+j]);
dNdxi[0][3] = -0.25*(1.+xi[1][3*i+j]);
dNdxi[1][3] = 0.25*(1.-xi[0][3*i+j]);
}
else if (nelnodes == 8)
{
dNdxi[0][0] = 0.25*(1.-xi[1][3*i+j])*(2.*xi[0][3*i+j]+xi[1][3*i+j]);
dNdxi[1][0] = 0.25*(1.-xi[0][3*i+j])*(xi[0][3*i+j]+2.*xi[1][3*i+j]);
dNdxi[0][1]= 0.25*(1.-xi[1][3*i+j])*(2.*xi[0][3*i+j]-xi[1][3*i+j]);
dNdxi[1][1] = 0.25*(1.+xi[0][3*i+j])*(2.*xi[1][3*i+j]-xi[0][3*i+j]);
dNdxi[0][2] = 0.25*(1.+xi[1][3*i+j])*(2.*xi[0][3*i+j]+xi[1][3*i+j]);
dNdxi[1][2] = 0.25*(1.+xi[0][3*i+j])*(2.*xi[1][3*i+j]+xi[0][3*i+j]);
dNdxi[0][3] = 0.25*(1.+xi[1][3*i+j])*(2.*xi[0][3*i+j]-xi[1][3*i+j]);
dNdxi[1][3] = 0.25*(1.-xi[0][3*i+j])*(2.*xi[1][3*i+j]-xi[0][3*i+j]);
dNdxi[0][4]= -xi[0][3*i+j]*(1.-xi[1][3*i+j]);
dNdxi[1][4] = -0.5*(1.-xi[0][3*i+j]*xi[0][3*i+j]);
dNdxi[0][5] = 0.5*(1.-xi[1][3*i+j]*xi[1][3*i+j]);
dNdxi[1][5] = -(1.+xi[0][3*i+j])*xi[1][3*i+j];
dNdxi[0][6] = -xi[0][3*i+j]*(1.+xi[1][3*i+j]);
dNdxi[1][6] = 0.5*(1.-xi[0][3*i+j]*xi[0][3*i+j]);
dNdxi[0][7] = -0.5*(1.-xi[1][3*i+j]*xi[1][3*i+j]);
dNdxi[1][7] = -(1.-xi[0][3*i+j])*xi[1][3*i+j];
}
// Compute the jacobian matrix && its determinant && jacobian inverse
for (f = 0;f<ncoord;f++)
{
for (c = 0;c<ncoord; c++)
{
dxdxi[f][c] = 0.;
for (a = 0;a<nelnodes;a++)
{
dxdxi[f][c]+= dNdxi[f][a]*coord[a][c];
}
}
}
Jdeterm = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
IJ[0][0]=dxdxi[1][1]/Jdeterm;
IJ[0][1]=-dxdxi[0][1]/Jdeterm;
IJ[1][0]=-dxdxi[1][0]/Jdeterm;
IJ[1][1]=dxdxi[0][0]/Jdeterm;
// Convert shape function derivatives:derivatives wrt global coords
for (f = 0;f<ncoord; f++)
{
for (a = 0;a<nelnodes; a++)
{
dNdx[f][a] = 0.;
for (c = 0;c<ncoord;c++)
{
dNdx[f][a] += IJ[f][c]*dNdxi[c][a];
}
}
}
B[0][0] = dNdx[0][0];
B[0][1] = 0.;
B[0][2] = dNdx[0][1];
B[0][3] = 0;
B[0][4] = dNdx[0][2];
B[0][5] = 0.;
B[0][6] = dNdx[0][3];
B[0][7] = 0.;
B[0][8] = dNdx[0][4];
B[0][9] = 0.;
B[0][10] = dNdx[0][5];
B[0][11] = 0;
B[0][12] = dNdx[0][6];
B[0][13] = 0.;
B[0][14] = dNdx[0][7];
B[0][15] = 0.;
B[1][0] = 0.;
B[1][1] = dNdx[1][0];
B[1][2] = 0;
B[1][3] = dNdx[1][1];
B[1][4] = 0.;
B[1][5] = dNdx[1][2];
B[1][6] = 0.;
B[1][7] = dNdx[1][3];
B[1][8] = 0.;
B[1][9] = dNdx[1][4];
B[1][10] = 0;
B[1][11] = dNdx[1][5];
B[1][12] = 0.;
B[1][13] = dNdx[1][6];
B[1][14] = 0.;
B[1][15] = dNdx[1][7];
B[2][0] = dNdx[1][0];
B[2][1] = dNdx[0][0];
B[2][2] = dNdx[1][1];
B[2][3] = dNdx[0][1];
B[2][4] = dNdx[1][2];
B[2][5] = dNdx[0][2];
B[2][6] = dNdx[1][3];
B[2][7] = dNdx[0][3];
B[2][8] = dNdx[1][4];
B[2][9] = dNdx[0][4];
B[2][10] = dNdx[1][5];
B[2][11] = dNdx[0][5];
B[2][12] = dNdx[1][6];
B[2][13] = dNdx[0][6];
B[2][14] = dNdx[1][7];
B[2][15] = dNdx[0][7];
BT[0][0] = dNdx[0][0];
BT[1][0] = 0.;
BT[2][0] = dNdx[0][1];
BT[3][0] = 0;
BT[4][0] = dNdx[0][2];
BT[5][0] = 0.;
BT[6][0] = dNdx[0][3];
BT[7][0] = 0.;
BT[8][0] = dNdx[0][4];
BT[9][0] = 0.;
BT[10][0] = dNdx[0][5];
BT[11][0] = 0;
BT[12][0] = dNdx[0][6];
BT[13][0] = 0.;
BT[14][0] = dNdx[0][7];
BT[15][0] = 0.;
BT[0][1] = 0.;
BT[1][1] = dNdx[1][0];
BT[2][1] = 0;
BT[3][1] = dNdx[1][1];
BT[4][1] = 0.;
BT[5][1] = dNdx[1][2];
BT[6][1] = 0.;
BT[7][1] = dNdx[1][3];
BT[8][1] = 0.;
BT[9][1] = dNdx[1][4];
BT[10][1] = 0;
BT[11][1] = dNdx[1][5];
BT[12][1] = 0.;
BT[13][1] = dNdx[1][6];
BT[14][1] = 0.;
BT[15][1] = dNdx[1][7];
BT[0][2] = dNdx[1][0];
BT[1][2] = dNdx[0][0];
BT[2][2] = dNdx[1][1];
BT[3][2] = dNdx[0][1];
BT[4][2] = dNdx[1][2];
BT[5][2] = dNdx[0][2];
BT[6][2] = dNdx[1][3];
BT[7][2] = dNdx[0][3];
BT[8][2] = dNdx[1][4];
BT[9][2] = dNdx[0][4];
BT[10][2] = dNdx[1][5];
BT[11][2] = dNdx[0][5];
BT[12][2] = dNdx[1][6];
BT[13][2] = dNdx[0][6];
BT[14][2] = dNdx[1][7];
BT[15][2] = dNdx[0][7];
//
// //================= SHAPE FUNCTIONS ==================================
// //
// // Calculates shape functions for various element types. No need here to
// calculate. We only need derivative of it.
// //
// // 2D Rectangular element
// if ( nelnodes == 4 )
// {
// N[1] = 0.25*(1.-xi[1])*(1.-xi[2]);
// N[2] = 0.25*(1.+xi[1])*(1.-xi[2]);
// N[3] = 0.25*(1.+xi[1])*(1.+xi[2]);
// N[4] = 0.25*(1.-xi[1])*(1.+xi[2]);
// }
// else if(nelnodes == 8)
// {
// N[1] = -0.25*(1.-xi[1])*(1.-xi[2])*(1.+xi[1]+xi[2]);
// N[2] = 0.25*(1.+xi[1])*(1.-xi[2])*(xi[1]-xi[2]-1.);
// N[3] = 0.25*(1.+xi[1])*(1.+xi[2])*(xi[1]+xi[2]-1.);
// N[4] = 0.25*(1.-xi[1])*(1.+xi[2])*(xi[2]-xi[1]-1.);
// N[5] = 0.5*(1.-xi[1]*xi[1])*(1.-xi[2]);
// N[6] = 0.5*(1.+xi[1])*(1.-xi[2]*xi[2]);
// N[7] = 0.5*(1.-xi[1]*xi[1])*(1.+xi[2]);
// N[8] = 0.5*(1.-xi[1])*(1.-xi[2]*xi[2]);
// }
//
//================= Material Stiffness ==================================
D[0][0]= E/(1-nu*nu);
D[0][1]= E*nu/(1-nu*nu);
D[0][2]= 0.;
D[1][0]= D[0][1];
D[1][1]=E/(1-nu*nu);
D[1][2]=0.;
D[2][0]=D[0][2];
D[2][1]=D[1][2];
D[2][2]=0.5*E*(1-nu)/(1-nu*nu);
for (g = 0;g<2*nelnodes; g++)
{
for (c = 0;c<3; c++)
{
BTD[g][c] = 0.;
for (k = 0;k<3;k++)
{
BTD[g][c] += BT[g][k]*D[k][c];
}
}
}
//================= ELEMENT STIFFNESS MATRIX ================================
for(l=0;l<2*nelnodes;l++)
{
for(h=0;h<2*nelnodes;h++)
{
BTDB[l][h]=0.;
for(m=0;m<3;m++)
{
BTDB[l][h]+= w1D[i]*w1D[j]*BTD[l][m]*B[m][h]*Jdeterm;
}
}
}
for(p=0;p<2*nelnodes;p++)
{
for(q=0;q<2*nelnodes;q++)
{
Kel[p][q]+= BTDB[p][q];
}
}
}
}
cout<<"\n\n\n";
cout<< " Please look at the output.data file to see the stiffness matrix"<<endl;
out<< " The Stiffness Matrix is:"<<endl<<endl;
for (i=0;i<2*nelnodes;i++)
{
for(j=0;j<2*nelnodes;j++)
out<<Kel[i][j]<<"\t"<<"\t";
out<<endl;
}
out.close();
}
