shape2d.h

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

////////////////////////////////////////////////////////////////////////////////
//
//	FN 函数：获取8（4）节点等参元的形函数
//
//  FNA 函数：获取形函数对X的倒数
//
//  FNB 函数：获取形函数对Y的倒数
//
//  Inv_Jaco函数: 获取等参变换中的Jacobi矩阵的逆矩阵
//
//////////////////////////////////////////////////////////////////////////////

#include <fstream.h>
#ifndef __SHAPE2D_
#define __SHAPE2D_
//关于平面二维坐标变换的程序实现
void FN(double A,double B,double *N,int *Ele);
void FNA(double A,double B,double *NA,int *Ele);
void FNB(double A,double B,double *NB,int *Ele);
double Inv_jaco(double A,double B,double *xo,double *yo,
                 double *Nx, double *Ny,int *Ele);

void FN(double A,double B,double *N,int *Ele)
/*参数说明
  A,B-局部坐标ξ,η;
    N-长度为8的数组的指针,存放8个插值基函数值;
  Ele-长度为8的数组的指针,存放插值区域的节点构成;
  如果元素为零.则意味该节点不存在.
*/
{if(Ele[4]) N[4]=0.5*(1-A*A)*(1-B);else N[4]=0.0;
 if(Ele[5]) N[5]=0.5*(1+A)*(1-B*B);else N[5]=0.0;
 if(Ele[6]) N[6]=0.5*(1-A*A)*(1+B);else N[6]=0.0;
 if(Ele[7]) N[7]=0.5*(1-A)*(1-B*B);else N[7]=0.0;

 N[0]=0.25*(1-A)*(1-B)-0.5*N[7]-0.5*N[4];
 N[1]=0.25*(1+A)*(1-B)-0.5*N[4]-0.5*N[5];
 N[2]=0.25*(1+A)*(1+B)-0.5*N[5]-0.5*N[6];
 N[3]=0.25*(1-A)*(1+B)-0.5*N[6]-0.5*N[7];

 return;

}
 void FNA(double A,double B,double *NA,int *Ele)
/*参数说明 
  A,B--局部坐标ξ,η;
   NA--长度为8的数组的指针，存放插值基函数偏导数值；
   Ele--长度为8的数组指针，存放插值区域的节点构成，如果某元素为零，则意味着该节点不存在；
*/
{ if(Ele[4]) NA[4]=-A*(1-B);else NA[4]=0.0;
  if(Ele[5]) NA[5]=0.5*(1-B*B);else NA[5]=0.0;
  if(Ele[6]) NA[6]=-A*(1+B);else NA[6]=0.0;
  if(Ele[7]) NA[7]=-0.5*(1-B*B);else NA[7]=0.0;

  NA[0]=-0.25*(1-B)-0.5*NA[7]-0.5*NA[4];
  NA[1]=0.25*(1-B)-0.5*NA[4]-0.5*NA[5];
  NA[2]=0.25*(1+B)-0.5*NA[5]-0.5*NA[6];
  NA[3]=-0.25*(1+B)-0.5*NA[6]-0.5*NA[7];
  return;
}

  void FNB(double A,double B,double *NB,int *Ele)
 /*参数说明
  A,B-局部坐标ξ,η;
  NB-长度为8的数组指针,存放插值基函数偏导数值;
  Ele-长度为8的数组指针,存放插值区域的节点构成,
       如果元素为零,则意味着该节点不存在;
*/
{if (Ele[4]) NB[4]=-0.5*(1-A*A);else NB[4]=0.0;
 if (Ele[5]) NB[5]=-B*(1+A);else NB[5]=0.0;
 if (Ele[6]) NB[6]=0.5*(1-A*A);else NB[6]=0.0;
 if (Ele[7]) NB[7]=-B*(1-A);else NB[7]=0.0;

 NB[0]=-0.25*(1-A)-0.5*NB[7]-0.5*NB[4];
 NB[1]=-0.25*(1+A)-0.5*NB[4]-0.5*NB[5];
 NB[2]=0.25*(1+A)-0.5*NB[5]-0.5*NB[6];
 NB[3]=0.25*(1-A)-0.5*NB[6]-0.5*NB[7];

 return;
}
 double Inv_jaco(double A,double B,double *xo,double *yo,
                 double *Nx, double *Ny,int *Ele)
/*参数说明
    A,B--局部坐标ξ,η;
     xo--长度为8的数组指针，存放插值节点总体x坐标;
	 yo--长度为8的数组指针，存放插值节点总体y坐标;
	 Nx--长度为8的数组指针，存放插值基函数对x偏导数；
	 Ny--长度为8的数组指针，插值基函数对y偏导数；
    Ele--长度为8的数组指针，存放插值区域的节点构成；如果元素为零，则意味着该节点不存在；
*/
 {    int    i,j; 
      double detj,temp,NA[8],NB[8],jaco[2][2];
      for(i=0;i<2;i++) 
		  for(j=0;j<2;j++) 
			  jaco[i][j]=0.0;
	  FNA(A,B,NA,Ele);  
	  FNB(A,B,NB,Ele);
	for(i=0;i<8;i++) 
		if(Ele[i])
		{   jaco[0][0]+=NA[i]*xo[i]; 
		    jaco[0][1]+=NA[i]*yo[i];
	        jaco[1][0]+=NB[i]*xo[i];  
			jaco[1][1]+=NB[i]*yo[i];
		}
    detj=jaco[0][0]*jaco[1][1]-jaco[0][1]*jaco[1][0];
	if(detj<1.0e-16)
	{    ofstream Out("file.err");  
	    Out<<"detj<=0\n"; 
		Out.close();
		exit(1);
	}
	 temp=jaco[0][0]/detj; 
	 jaco[0][0]=jaco[1][1]/detj; 
	 jaco[1][1]=temp;
	 jaco[0][1]=-jaco[0][1]/detj; 
	 jaco[1][0]=-jaco[1][0]/detj;
	 for(i=0;i<8;i++) 
		 if(Ele[i])
		 {   Nx[i]=jaco[0][0]*NA[i]+jaco[0][1]*NB[i];
	        Ny[i]=jaco[1][0]*NA[i]+jaco[1][1]*NB[i];
		 }
      return(detj);
 }

#endif