有限元六面体单元 C语言源代码.rar

void Getnumber(double nodeCoo[8][3],double elementK[24][3]) { double Gauss[2]={0.57735026919,-0.57735026919};//高斯点 double WeiG[2]={1.0,1.0};//权重 int local[3][8]={{1,1,-1,-1,1,1,-1-1},{-1,1,1,-1,-1,1,1,-1},{-1,-1,-1,-1,1,1,1,1}};//结点局部坐标 double diffN[3][8];//插值函数关于局部坐标偏导数 double diffNlnverJ[3][8]={0};//插值函数关于整体坐标偏导数 double DJ; double E; double v; double J[3][3]; double *MatFirstadd; double lnverJ[3][3];//存储J的逆矩阵 double B[6][24]={0};//存储B矩阵 double Blnver[24][6];//存储B的转置 double *MatrixinverFirstadd; double BlnverD[24][6]={0};//存储B的转置与D的乘积 int i,j,ksi,eta,zeta,m,n; double Gaussksi,Gausseta,Gausszeta; double G=E/(2*(1+v)),lam=E*v/((1+v)*(1-2*v)); //已知常数E,v double D[6][6]={{lam+2*G,lam,lam,0,0,0},{lam,lam+2*G,lam,0,0,0},{lam,lam,lam+2*G,0,0,0},{0,0,0,G,0,0},{0,0,0,0,G,0},{0,0,0,0,0,G}};//D矩阵 double BTDBJ[24][24]={0}; { for(ksi=0;ksi<2;ksi++) //以ksi方向高斯点计算积分 { Gaussksi=Gauss[ksi]; for(eta=0;eta<2;eta++) //以eta方向高斯点计算积分 { Gausseta=Gauss[eta]; for(zeta=0;zeta<2;zeta++); //以zeta方向高斯点计算积分 { Gausszeta=Gauss[zeta]; for(i=0;i<3;i++) //以坐标进行循环 { if(i==0) { for(j=0;j<8;j++) //以节点进行循环 { diffN[0][j]=(1.0/8)*(1+Gausseta*local[1][j])*(1+Gausszeta*local[2][j]); //求插值函数关于kis的偏函数 } } else if(i==1) { for(j=0;j<8;j++) //以节点进行循环 { diffN[1][j]=(1.0/8)*Gausseta*(1+Gaussksi*local[0][j])*(1+Gausszeta*local[2][j]); //求插值函数关于eta的偏导数 } } else { for(j=0;j<8;j++) //以节点进行循环 { diffN[2][j]=(1.0/8)*Gausszeta*(1+Gaussksi*local[0][j])*(1+Gausseta*local[1][j]); //求插值函数关于zeta的偏导数 } } } for(i=0;i<3;i++) //计算雅可比矩阵 { for(j=0;j<3;j++) { J[i][j]=0; for(m=0;m<8;m++) { J[i][j]=J[i][j]+diffN[i][m]*nodeCoo[m][j]; } } } m=3; DJ=Surplus(J[0],m,m); //计算行列式j的值 MatFirstadd=MatrixOpp(J[0],m,m);//返回矩阵就[3][3]逆矩阵的首地址 for(i=0;i<3;i++) { for(j=0;j<3;j++) { lnverJ[i] [j] =*(MatFirstadd+i*3+j);//利用返回的首地址求出逆矩阵 } } for(i=0;i<3;i++) { for(j=0;j<8;j++) { for(m=0;m<3;m++) { diffNlnverJ[i][j]=diffNlnverJ[i][j]+lnverJ[i][m]*diffN[m][j];//计算插值函数关于整体坐标xyz的偏导数 } } } for(j=0;j<8;j++)//计算B矩阵 { B[0][j*3]=diffNlnverJ [0][j]; B[1][1+j*3]=diffNlnverJ[1][j]; B[2][2+j*3]=diffNlnverJ[2][j]; B[3][3*j]=B[1][1+3*j]; B[3][1+j*3]=B[0][j*3]; B[4][1+j*3]=B[2][2+j*3]; B[4][2+j*3]=B[1][1+j*3]; B[5][j*3]=B[2][2+j*3]; B[5][2+j*3]=B[0][j*3]; } m=6; n=24; MatrixlnverJFirstadd=Matrixlnver(B[0],m,n); for(i=0;i<24;i++) { for(j=0;j<6;j++) { Blnver[i][j]=*(MatrixlnverFirstadd+i*6+j);//利用返回首地址求B矩阵的转置 } } for(i=0;i<24;i++ ) { for(j=0;j<6;j++) { for(m=0;m<0;m++) { BlnverD[i][j]=BlnverD[i][j]+Blnver[i][m]*D[m][j];//计算B的转置与D的乘积 } } } for(i=0;i<24;i++) { for(j=0;j<24;j++) { for(m=0;m<6;m++) { BTDBJ[i][j]=BTDBJ[i][j]+BlnverD[i][m]*B[m][j]*DJ;//计算B的转置，D与B，雅可比矩阵的乘积 } } } for(i=0;i<24;i++) { for(j=0;j<24;j++) { elementK[i][j]=BTDBJ[i][j]+elementK[i][j]; } } } } } } }