main.cpp

计算雅克比行列式的小程序

#include <stdio.h>
#include <stdlib.h>
#include <iomanip.h>
#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 *NA,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的数组的指针，存放八个插值基函数值；
    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]==1)  NA[4]=-A*(1-B);   else NA[4]=0.0;
  if(Ele[5]==1)  NA[5]=0.5*(1-B*B);else NA[5]=0.0;
  if(Ele[6]==1)  NA[6]=-A*(1+B);   else NA[6]=0.0;
  if(Ele[7]==1)  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]==1)  NB[4]=-0.5*(1-A*A);else NB[4]=0.0;
  if(Ele[5]==1)  NB[5]=-B*(1+A);    else NB[5]=0.0;
  if(Ele[6]==1)  NB[6]=0.5*(1-A*A); else NB[6]=0.0;
  if(Ele[7]==1)  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;


main()
{
	int i,j;
	int Ele[8];
	float a;
    float b;
	double NA[8];
	double NB[8];
	double Nx[8];
	double Ny[8];
	double xo[4]={1,4,7,2};
	double yo[4]={1,1.5,5.5,4.5};
	double detj;
    

	for(i=4;i<=7;i++)
		{
			Ele[i]=-1;
		}
	
	printf("please enter a and b:");
    scanf("%f,%f",&a,&b);
	printf("a=%lf,b=%lf\n",a,b);
	FNA(a,b,NA,Ele); 	
	FNB(a,b,NB,Ele);
	for(i=0;i<=3;i++)
	{
		printf("NA[%d]=%f,NB[%d]=%f\n",i,NA[i],i,NB[i]);
	}
    
    detj=Inv_jaco(a,b,xo,yo,Nx,Ny,Ele);
	printf("detj=%f\n",detj);

    return 0;
}