cacubf.h

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

//////////////////////////////////////////////////////////////////////////
//
//        CacuBF函数：计算结构的体积力
//   Density--单元密度;
//   tt     --单元厚度,仅用于平面应力情况下;
//   ax     --,若IK==2,则为转动加速度;
//   ay     --在y方向的惯性加速度;
//   xo     --长度为8的数组指针,存放单元节点的x坐标;
//   yo     --长度为8的数组指针,存放单元节点的y坐标;
//   MLoad  --长度为NF+1的数组指针,存放节点的等荷荷载;
//   Ele    --长度为8的数组指针,存放单元节点编号;
//   nGauss --nGauss!2,采用3*3的高斯积分,否则采用2*2的高斯积分.
//////////////////////////////////////////////////////////////

#ifndef __CACUBF_
#define __CACUBF_
#include "shape2d.h"
void CacuBF(double Density,double tt,double ax,double ay,double *xo,
			double *yo,double *MLoad,int *Ele,int IK,int nGauss);
                     
void CacuBF(double Density,double tt,double ax,double ay,double *xo,
			double *yo,double *MLoad,int *Ele,int IK,int nGauss)
{     int    i,ii,jj,ij,je[16];
      double detj,rx,fx,fy,N[8],Nx[8],Ny[8],Pe[16],AB[3],H[3];
     if(nGauss!=2)           //设定高斯积分点坐标及权系数
	 {   AB[0]=-0.7745966692; 
	     AB[1]=0; 
		 AB[2]=-AB[0];
         H[0]=5.0/9.0; 
		 H[1]=8.0/9.0; 
		 H[2]=H[0]; 
		 ij=1;
	 }
          else
		  {  AB[0]=-0.5773302692;
		     AB[2]=-AB[0];
			 H[0]=1;
			 H[2]=1;
			 ij=2;
		  }
    for(i=0;i<8;i++)
      if(Ele[i])
	  {  je[2*i]=2*Ele[i]-1;
	     je[2*i+1]=2*Ele[i];
	  }
       else
	   {  je[2*i]=0;
	      je[2*i+1]=0;
	   }                            // 形成单元自由度的总体编号
  for(i=0;i<16;i++)
  Pe[i]=0.0;
  for(ii=0;ii<3;ii=ii+ij)                                    //临时数组清零
  for(jj=0;jj<3;jj=jj+ij)                                //进行高斯积分
  {  
	  for(i=0;i<8;i++)
	  {   N[i]=0.0;
	     Nx[i]=0.0;
		 Ny[i]=0.0;
	  }                         //清零
    detj=Inv_jaco(AB[ii],AB[jj],xo,yo,Nx,Ny,Ele)*H[ii]*H[jj];
    detj=detj*tt;
    FN(AB[ii],AB[jj],N,Ele);                               // 计算插值函数
     if(IK==2)
	 {  rx=0.0;
       for(i=0;i<8;i++)
         if(Ele[i])
             rx+=xo[i]*N[i];
       fx=Density*ax*ax*rx;
	   detj=detj*2.0*3.1415926*rx;
	 }
      else 
		  fx=Density*ax;                                         // 计算 x方向体积力分量
          fy=Density*ay;                                          // 计算 y方向体积力分量
    for(i=0;i<8;i++)
    if(Ele[i])                                                  // 计算体积力的节点等效载荷  
	{   Pe[2*i]+=N[i]*fx*detj;
     Pe[2*i+1]+=N[i]*fy*detj;
	}
  }                                                      //完成体积力的节点的等效载荷的积分计算
 for(i=0;i<16;i++)
 if(je[i])
 MLoad[je[i]]+=Pe[i];
}

#endif