cacuss.h

这是等参单元的有限元程序

//////////////////////////////////////////////////////////
//    Cass函数是计算单元内节点应力
//
////////////////////////////////////////////////////////
#include "shape2d.h"
#include "eld.h"
void Cass(double E,double u,double *xo,double *yo,double *uo,double *vo,
			double (*SS)[5],int *Ele,int IK);

void Cass(double E,double u,double *xo,double *yo,double *uo,double *vo,
			double (*SS)[5],int *Ele,int IK)

{    int    i,j,LL,ii;
     double Nx[8],Ny[8],D[4][4],Strain[4];
     double A[8]={-1.0,1.0,1.0,-1.0,0.0,1.0,0.0,-1.0};
	 double B[8]={-1.0,-1.0,1.0,1.0,-1.0,0.0,1.0,0.0};  //设定各节点局部坐标　
     if(IK==2) LL=4; 
        else LL=3;
     eld(E,u,D,IK);   //计算D矩阵
     for(ii=0;ii<8;ii++) 
	 if(Ele[ii])             //计算节点应力
	 {
		 for(i=0;i<8;i++)
		 {  Nx[i]=0.0;
		    Ny[i]=0.0;
		 }          //清零
        Inv_jaco(A[ii],B[ii],xo,yo,Nx,Ny,Ele);         //坐标变换
       for(i=0;i<4;i++) Strain[i]=0.0;
       for(i=0;i<8;i++) Strain[0]+=Nx[i]*uo[i];             //x向应变
       for(i=0;i<8;i++) Strain[1]+=Ny[i]*vo[i];             //y向应变
       for(i=0;i<8;i++) Strain[2]+=Nx[i]*vo[i]+Ny[i]*uo[i];    //剪应变
        if(IK==2)
			if(fabs(xo[ii])>1.0e-6) 
				Strain[3]=uo[ii]/xo[ii];    //周向应变
               else Strain[3]=Strain[0];
       for(i=0;i<LL;i++)
		   for( j=0;j<LL;j++)
              SS[Ele[ii]][i+1]+=D[i][j]*Strain[j];          //计算总体编号为Ele[ii]的节点应力,并将各分量叠加
         SS[Ele[ii]][0]+=1;       //将Ele[ii]节点的引用次数叠加
	 }
	 return;
}