1.h

有限单元法中笛卡尔坐标系下六面体单元刚度阵

#include <math.h>
#include <iostream.h>
#define IPERM(i)    ((i==3)?1:(i+1))//为计算行列式时，余子式所对应的顺序
const double XK[4]={0,-0.77459667,0.0,0.77459667};//高斯点坐标
const double WGT[4]={0,0.55555556,0.88888889,0.55555556};//高斯点系数加权值
class calculate_hexahedron_element
{private:
 double XX[9][4],U,S[25][25],E;
 double Ni[9],//形函数插值函数  
        Pi[4][9];//对局部坐标的导数值
 int i,j,l,lx,ly,lz;
 double E1,E2,E3;
 double A[4][4];//雅可比行列式
 double B[4][4];//伴随矩阵
 double SA[9][4];
 double FAC,CC1,CC2,CC3,C,DET,GT,G,C1,C2,C3,UI,VI,WI,UJ,VJ,WJ,UU,VV,WW,UV,VU,UW,WU,WV,VW;
 int K3,K2,K1,L1,L2,L3;
 int I,J,K;
public:
void  input(double xx[9][4]/*八节点在整体坐标系下的坐标分量*/,
			    //double S[25][25]/*八节点单元的单刚由此数组存储*/，
				double u/*泊松比*/,double e/*弹性模量*/)
{for(i=1;i<=8;i++) 
     for(j=1;j<=3;j++)   
	 {XX[i][j]=xx[i][j];}//cout<<XX[i][j]<<endl;
 U=u;E=e;
 for(i=1;i<=25;i++)
 {for(j=1;j<=24;j++)  S[i][j]=0.0;
 }
}
void DERIV()//计算形函数插值函数和对局部坐标的导数值 
{Ni[1]=(1+E1)*(1+E2)*(1+E3)/8;
 Ni[2]=(1+E1)*(1-E2)*(1+E3)/8;
 Ni[3]=(1+E1)*(1-E2)*(1-E3)/8;
 Ni[4]=(1+E1)*(1+E2)*(1-E3)/8;
 Ni[5]=(1-E1)*(1+E2)*(1+E3)/8;
 Ni[6]=(1-E1)*(1-E2)*(1+E3)/8;
 Ni[7]=(1-E1)*(1-E2)*(1-E3)/8;
 Ni[8]=(1-E1)*(1+E2)*(1-E3)/8;

 Pi[1][1]=(1+E2)*(1+E3)/8;
 Pi[1][2]=(1-E2)*(1+E3)/8;
 Pi[1][3]=(1-E2)*(1-E3)/8;
 Pi[1][4]=(1+E2)*(1-E3)/8;
 Pi[1][5]=-(1+E2)*(1+E3)/8;
 Pi[1][6]=-(1-E2)*(1+E3)/8;
 Pi[1][7]=-(1-E2)*(1-E3)/8;
 Pi[1][8]=-(1+E2)*(1-E3)/8;

 Pi[2][1]=(1+E1)*(1+E3)/8;
 Pi[2][2]=-(1+E1)*(1+E3)/8;
 Pi[2][3]=-(1+E1)*(1-E3)/8;
 Pi[2][4]=(1+E1)*(1-E3)/8;
 Pi[2][5]=(1-E1)*(1+E3)/8;
 Pi[2][6]=-(1-E1)*(1+E3)/8;
 Pi[2][7]=-(1-E1)*(1-E3)/8;
 Pi[2][8]=(1-E1)*(1-E3)/8;

 Pi[3][1]=(1+E1)*(1+E2)/8;
 Pi[3][2]=(1+E1)*(1-E2)/8;
 Pi[3][3]=-(1+E1)*(1-E2)/8;
 Pi[3][4]=-(1+E1)*(1+E2)/8;
 Pi[3][5]=(1-E1)*(1+E2)/8;
 Pi[3][6]=(1-E1)*(1-E2)/8;
 Pi[3][7]=-(1-E1)*(1-E2)/8;
 Pi[3][8]=-(1-E1)*(1+E2)/8;

for(i=1;i<=3;i++)
    for(j=1;j<=3;j++)
	{C=0.0;
	    for(l=1;l<=8;l++)
		{C+=Pi[i][l]*XX[l][j];}
	 A[i][j]=C;
	}
for(I=1;I<=3;I++)
{J=IPERM(I);K=IPERM(J);
 B[I][I]=A[J][J]*A[K][K]-A[K][J]*A[J][K];
 B[I][J]=A[K][J]*A[I][K]-A[I][J]*A[K][K];
 B[J][I]=A[J][K]*A[K][I]-A[J][I]*A[K][K];
}
/*/************************
for(i=1;i<=3;i++)
{for(j=1;j<=8;j++)
 cout<<Pi[i][j]<<"  ";
cout<<endl;
}
*///************************
DET=A[1][1]*B[1][1]+A[1][2]*B[2][1]+A[1][3]*B[3][1];
for(I=1;I<=3;I++)
{for(J=1;J<=8;J++)
	{C=0.0;
	    for(K=1;K<=3;K++)
		{C+=B[I][K]*Pi[K][J];}
	 SA[J][I]=C/DET;
//cout<<SA[J][I]<<"  "<<DET<<endl;
	}
}
//cout<<DET<<endl;
/*/8************************************
for(i=1;i<=3;i++) 
{for(j=1;j<=3;j++)
{cout<<Ni[i]<<" ";}
cout<<endl;
}
*///**************************************
}


void	STIFST()//计算八节点单元的单刚
{FAC=E/((1-2*U)*(1+U));
 CC1=1-U;CC2=U;CC3=0.5-U;
 int aaa=0;
 for(lx=1;lx<=3;lx++)  
 {E1=XK[lx];
     for(ly=1;ly<=3;ly++)
	 {E2=XK[ly];
         for(lz=1;lz<=3;lz++)
		 {E3=XK[lz];
		 DERIV();
         GT=WGT[lx]*WGT[ly]*WGT[lz];
		 G=GT*FAC;C1=G*CC1;C2=G*CC2;C3=G*CC3;
		      for(I=1;I<=8;I++)
			  {K3=3*I;K2=K3-1;K1=K2-1;
		       UI=SA[I][1];VI=SA[I][2];WI=SA[I][3];
		         for(J=I;J<=8;J++)
//原程序为	  for(J=1;J<=8;J++)		
				 {L3=3*J;L2=L3-1;L1=L2-1;
				  UJ=SA[J][1];VJ=SA[J][2];WJ=SA[J][3];
				  UU=UI*UJ;VV=VI*VJ;WW=WI*WJ;
				  UV=UI*VJ;VU=VI*UJ;
				  UW=UI*WJ;WU=WI*UJ;
				  VW=VI*WJ;WV=WI*VJ;
                  S[K1][L1]=S[K1][L1]+C1*UU+C3*(VV+WW);
				  S[K2][L2]=S[K2][L2]+C1*VV+C3*(UU+WW);
				  S[K3][L3]=S[K3][L3]+C1*WW+C3*(VV+UU);
				  //*********************************8
				 // if(I==8&&J==8) {aaa++;cout<<K3<<" "<<L3<<" "<<aaa<<" "<<C1<<" "<<WW<<"  "<<C3<<"  "<<VV<<"  "<<UU<<endl;}
				  //**********************************
				  S[K1][L2]=S[K1][L2]+C2*UV+C3*VU;
				  S[K1][L3]=S[K1][L3]+C2*UW+C3*WU;
				  S[K2][L3]=S[K2][L3]+C2*VW+C3*WV;
				  if(I==J) continue;
				  S[K2][L1]=S[K2][L1]+C2*VU+C3*UV;
				  S[K3][L1]=S[K3][L1]+C2*WU+C3*UW;
				  S[K3][L2]=S[K3][L2]+C2*WV+C3*VW;
				 }
			  }
		 }
	 }

 }
for(I=1;I<=24;I++)
    for(J=I;J<=24;J++)
	{S[J][I]=S[I][J];}
}

void output()
{for(I=1;I<=24;I++)
{    for(J=1;J<=24;J++) {cout<<S[I][J]<<endl;}
   cout<<I<<endl;cout<<endl;
}
}
};