cacusf.h

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

///////////////////////////////////////////////////////////////////////////////
//
//qritz 函数：
//
//		功能：处理分布载荷
//
//      变量：xx0,yy0：分布力作用的单元表面三个或两个节点的坐标
//
//            q0: 分布力的大小（作用在节点上的值）
//
//            pload：载荷向量
//
//            tt：单元厚度
//
//            Node:分布力作用的单元表面节点号
//
//            IK: 指示是平面问题还是轴对称问题，＝2时是轴对称问题											
//
///////////////////////////////////////////////////////////////////////////////
#ifndef __CACUSF_
#define __CACUSF_
#include "shape1d.h"
  void CacuSF(double tt,double *qo,double *xo,double *yo,double *pload,
			int *Node,int IK,int nGauss);
  
  void CacuSF(double tt,double *qo,double *xo,double *yo,double *pload,
			int *Node,int IK,int nGauss)
{ int i,ii,ij,je[6];
  double detj,rx,xa,ya,P,qx,qy,N[3],Na[3],Pe[6],A[3],H[3];
  if(nGauss==3)//设置高斯积分点及权系数;
  {
    A[0]=-0.7745966692; A[1]=0;  A[2]=-A[0];
	H[0]=5.0/9.0;       H[1]=8.0/9.0; H[2]=H[0]; ij=1;
  }
  else
  {
	A[0]=-0.5773502962; A[2]=-A[0]; H[0]=1.0;H[2]=1.0;ij=2;
  }
  for(i=0;i<3;i++)  if(Node[i])  { je[2*i]=2*Node[i]-1;  je[2*i+1]=2*Node[i];}
  else
  { je[2*i]=0;  je[2*i+1]=0;}//形成单元边界自由度的总体编号;
  for(i=0;i<6;i++)  Pe[i]=0.0;
  for(ii=0;ii<3;ii=ii+ij)//进行高斯积分;
  {
  shape(A[ii],N,Node);//进行一维插值基函数;
  shape_a(A[ii],Na,Node);//计算一维插值基函数局部坐标的导数;
  P=0;xa=0;ya=0;for(i=0;i<3;i++)  if(Node[i])
  {P+=qo[i]*N[i];xa+=xo[i]*Na[i];ya+=yo[i]*Na[i];}
  if(IK==2)
  { rx=0.0;for(i=0;i<3;i++)  if(Node[i]) rx+=xo[i]*N[i];
    detj=2.0*3.1415926*rx*H[ii];
   }
  else detj=H[ii]*tt;//计算积分加权因子;
  qx=-P*ya;qy=P*xa;//计算压力的x分量和y分量;
  for(i=0;i<3;i++)  if(Node[i])
  {Pe[2*i]+=N[i]*qx*detj;Pe[2*i+1]+=N[i]*qy*detj;}
   }
  for(i=0;i<6;i++)  if(je[i]) pload[je[i]]+=Pe[i];
}

#endif