femtemper.cpp

一个温度场计算的小程序 望各位大虾指点一二

#include <iostream>
#include <stdio.h>
#include <fstream>
#include <iomanip>
#include <cmath>
using namespace std;

const int NN=22;
const int NE=16;

#ifndef ELEMENT_f_MATRIX_H
#define ELEMENT_f_MATRIX_H
   void fmatrix(int NOE,int NGN[][4],float CN[][2],int sbound[][3],double f[],double p[4],
	            float E,float v,float aa);
#endif

#ifndef GLOBLE_STIFFNESS_MATRIX_H
#define GLOBLE_STIFFNESS_MATRIX_H
        void GlbMatrix( int NGN[][4], float CN[][2],
                                 float T0, float E,float v,float aa,
                                 int sbound[][3],
								 double H[][NN], double R[][NN], double F[NN],double P[NN]);
#endif

#ifndef ELEMENT_h_MATRIX_H
#define ELEMENT_h_MATRIX_H
   void hmatrix(int NOE,int NGN[][4],float CN[][2],double h[][4],
	            double g[][4],int sbound[][3],float E,float v,float aa);
#endif

#ifndef ELEMENT_r_MATRIX_H
#define ELEMENT_r_MATRIX_H
   void rmatrix(int NOE,int NGN[][4],float CN[][2],
	             double r[][4],float E,float v,float aa);
#endif

#ifndef  TEMPERATURE_H
#define  TEMPERATURE_H
              void adjust(double GK[][NN],int fbound[NN],float Tb,double FR[NN],double temper[NN]);
              void gauss(double GK[][NN],double temper[NN]);
              void temperature(double H[][NN],double R[][NN],double F[NN],double P[NN],
				 float T0,float Ta,float Tb,int fbound[NN],float time,
				 int d,double temper[NN]);
#endif

int main(int argc, char* argv[])
{
	const int NN=22;
    const int NE=16;
	int i,j;
    int d=0;//NE为单元数，NN为结点数
    float T0=0.0,time=0.0,E=0.0,v=0.0,aa=0.0;//E为导热系数，v为表面放热系数,a导温系数
    
	 int n=80; //每行最多读入80 个字符
   char *newlin=new char[n];
   FILE *fp ;
   fp = fopen("Input.txt", "r") ;
   fscanf(fp,"%s",&newlin);
   fscanf(fp, "%f", &T0) ;
   fscanf(fp, "%f", &time) ;
   fscanf(fp, "%d", &d) ;
   fscanf(fp, "%f", &E) ;
   fscanf(fp, "%f", &v) ;
   fscanf(fp, "%f", &aa) ;
   int (*NGN)[4]=new int [NE][4];
   float (*CN)[2]=new float [NN][2];
   fscanf(fp,"%s",&newlin); //遇到字符串，跳行；可以跳过.input.txt.中的注释
   
   for(i=0;i<NE;i++)
      for(j=0;j<4;j++)
	  {
		  fscanf(fp, "%d", &NGN[i][j]) ;
	  }
   fscanf(fp,"%s",&newlin);
   for(i=0;i<NN;i++)
      for(j=0;j<2;j++)
	  {
		  fscanf(fp, "%f", &CN[i][j]) ;
	  }
     
	  float Tb=0;
	  int *fbound=new int[NN];//一类边界
	  fscanf(fp,"%s",&newlin);
	  fscanf(fp, "%f", &Tb) ;
      for(i=0;i<NN;i++)
      {
		 fscanf(fp, "%d", &fbound[i]) ;
	  }
      
	  float Ta=0.0;
	  int(*sbound)[3]=new int[NN][3];//二类边界
	  fscanf(fp,"%s",&newlin);
      fscanf(fp, "%f", &Ta) ;
	  for(i=0;i<NE;i++)
		  for(j=0;j<3;j++)
		  {
               fscanf(fp, "%d", &sbound[i][j]) ;
		  }
	  fclose(fp);
      
	   double H[NN][NN];
	   double R[NN][NN];
	   double F[NN];
	   double P[NN];
	   double temper[NN];
       int *count=new int[NN];
	   for(i=0;i<NN;i++)
	        for(j=0;j<NN;j++)
		          H[i][j]=0.0;
	   for(i=0;i<NN;i++)
	        for(j=0;j<NN;j++)
		          R[i][j]=0.0;
       for(i=0;i<NN;i++)
	   {
		    F[i]=0.0;
	        P[i]=0.0;
	   }
	   for(i=0;i<NN;i++)
		    temper[i]=0.0;
	    printf("**********************************************************\n");
        printf("*****Welcome to Temperature Finite Element Analysis!***********\n");
        printf("**********************************************************\n");
           
        GlbMatrix(NGN,CN,T0,E,v,aa,sbound,H,R,F,P);
           
        temperature(H,R,F,P,T0,Ta,Tb,fbound,time,d,temper);
        return 0;
}

void fmatrix(int NOE,int NGN[][4],float CN[][2],int sbound[][3],double f[],double p[],
	            float E,float v,float aa)
{
     int i=0,j=0,k=0,m=0,n=0;
	 int ngn[4];
         for(i=0;i<4;i++)ngn[i]=NGN[NOE-1][i];
     float cn[4][2];
         for(i=0;i<4;i++)
             for(j=0;j<2;j++)
                   cn[i][j]=CN[ngn[i]-1][j];
	 float a[3]={0,0,0};
	 float b[3]={0,0,0};
	 a[0]=((-1)*cn[0][0]+cn[1][0]+cn[2][0]-cn[3][0])/4;
     a[1]=((-1)*cn[0][0]-cn[1][0]+cn[2][0]+cn[3][0])/4;
     a[2]=(cn[0][0]-cn[1][0]+cn[2][0]-cn[3][0])/4;
     b[0]=((-1)*cn[0][1]+cn[1][1]+cn[2][1]-cn[3][1])/4;
     b[1]=((-1)*cn[0][1]-cn[1][1]+cn[2][1]+cn[3][1])/4;
     b[2]=(cn[0][1]-cn[1][1]+cn[2][1]-cn[3][1])/4;
	 double gusa[2]={-0.5773503,0.5773503};
     double gusb[2]={-0.5773503,0.5773503};
     double jac=0.0;
	 float u[4]={-1,1,1,-1};
	 float w[4]={-1,-1,1,1}; 
	 for(m=0;m<2;m++)
		 for(n=0;n<2;n++)
		 {
              jac=(a[0]+a[2]*gusb[n])*(b[1]+b[2]*gusa[m])-(a[1]+a[2]*gusa[m])*(b[0]+b[2]*gusb[n]);
              for(i=0;i<4;i++)
			  {
				  f[i]+=((1+u[i]*gusa[m])*(1+w[i]*gusb[n])*jac)/(4*aa);
			  }
		 }
     double s=0.0;
	 if(sbound[NOE-1][0]==0)
	 {
		 for(i=0;i<4;i++)
		 {
			 p[i]=0.0;
		 }
	 }
	 else
	 {
        s=(v/E)*sqrt(pow((CN[sbound[NOE-1][1]-1][0]-CN[sbound[NOE-1][2]-1][0]),2)+pow((CN[sbound[NOE-1][1]-1][1]-CN[sbound[NOE-1][2]-1][1]),2)); 
 
		for(i=0;i<4;i++)
		{
			if(ngn[i]==sbound[NOE-1][1])
			{
				p[i]=1.0/2.0;
			}
			else
			{
				if(ngn[i]==sbound[NOE-1][2])
				{
					p[i]=1.0/2.0;
				}
			}
		}
	 }
	 for(i=0;i<4;i++)
	 {
		 p[i]=p[i]*s;
	 }
}

void GlbMatrix( int NGN[][4], float CN[][2],
               float T0, float E,float v,float aa,  int sbound[][3],
			   double H[][NN], double R[][NN],double F[NN],double P[NN])
{
    int i,j;
       int NOE=1;
       double h[4][4],g[4][4],k[4][4],r[4][4],f[4],p[4];
	   for(i=0;i<4;i++)
		   for(j=0;j<4;j++)
		   {
			   h[i][j]=0.0;
			   g[i][j]=0.0;
			   k[i][j]=0.0;
			   r[i][j]=0.0;
		   }
	   for(i=0;i<4;i++)
	   {
		   f[i]=0.0;
	       p[i]=0.0;
	   }
       int gr=0,gc=0;
       int ngn[4];
       double temp=0.0;
	   for(NOE=1;NOE<=NE;NOE++)
	   {
           for(i=0;i<4;i++)ngn[i]=NGN[NOE-1][i]-1;
                  hmatrix(NOE,NGN,CN,h,g,sbound,E,v,aa);
           for(i=0;i<4;i++)
		   {
                 gr=ngn[i];
                 for(j=0;j<4;j++)
				 {
                       gc=ngn[j];
                       temp=h[i][j]+g[i][j];
                       H[gr][gc]+=temp;
				 }
		   }
	   }
       for(NOE=1;NOE<=NE;NOE++)
	   {
           for(i=0;i<4;i++)ngn[i]=NGN[NOE-1][i]-1;
                  rmatrix(NOE,NGN,CN,r,E,v,aa);
           for(i=0;i<4;i++)
		   {
                 gr=ngn[i];
                 for(j=0;j<4;j++)
				 {
                       gc=ngn[j];
                       R[gr][gc]+=r[i][j];
				 }
		   }
	   }
	   for(NOE=1;NOE<=NE;NOE++)
	   {
               for(i=0;i<4;i++)
			   {
			       ngn[i]=NGN[NOE-1][i]-1;
			   }
               fmatrix(NOE,NGN,CN,sbound,f,p,E,v,aa); 
			   for(i=0;i<4;i++)
			   {
				   gr=ngn[i];
				   F[gr]+=f[i];
				   P[gr]+=p[i];
			   }
	   }
}

void hmatrix(int NOE,int NGN[][4],float CN[][2],double h[][4],
	            double g[][4],int sbound[][3],float E,float v,float aa)
{
     int i=0,j=0,k=0,m=0,n=0;
	 int ngn[4];
         for(i=0;i<4;i++)ngn[i]=NGN[NOE-1][i];
     float cn[4][2];
         for(i=0;i<4;i++)
             for(j=0;j<2;j++)
                   cn[i][j]=CN[ngn[i]-1][j];
	 float a[3]={0,0,0};
	 float b[3]={0,0,0};
	 for( k=0;k<4;k++)
          for(i=0;i<4;i++)
		  {
			  h[k][i]=0;
		  }
	 a[0]=((-1)*cn[0][0]+cn[1][0]+cn[2][0]-cn[3][0])/4;
     a[1]=((-1)*cn[0][0]-cn[1][0]+cn[2][0]+cn[3][0])/4;
     a[2]=(cn[0][0]-cn[1][0]+cn[2][0]-cn[3][0])/4;
     b[0]=((-1)*cn[0][1]+cn[1][1]+cn[2][1]-cn[3][1])/4;
     b[1]=((-1)*cn[0][1]-cn[1][1]+cn[2][1]+cn[3][1])/4;
     b[2]=(cn[0][1]-cn[1][1]+cn[2][1]-cn[3][1])/4;
	 double gusa[2]={-0.5773503,0.5773503};
     double gusb[2]={-0.5773503,0.5773503};
     double jac=0;
	 float u[4]={-1,1,1,-1};
	 float w[4]={-1,-1,1,1};
	 for(m=0;m<2;m++)
	     for(n=0;n<2;n++)
		 {
			 jac=(a[0]+a[2]*gusb[n])*(b[1]+b[2]*gusa[m])-(a[1]+a[2]*gusa[m])*(b[0]+b[2]*gusb[n]);
			 for(i=0;i<4;i++)
			       for(j=0;j<4;j++)
				   {
                         h[i][j]+=((u[i]*(b[1]+b[2]*gusa[m])*(1+w[i]*gusb[n])
		                           -w[i]*(b[0]+b[2]*gusb[n])*(1+u[i]*gusa[m]))
						          *(u[j]*(b[1]+b[2]*gusa[m])*(1+w[j]*gusb[n])
		                           -w[j]*(b[0]+b[2]*gusb[n])*(1+u[j]*gusa[m]))
								  +(w[i]*(a[0]+a[2]*gusb[n])*(1+u[i]*gusa[m])
								   -u[i]*(b[1]+b[2]*gusa[m])*(1+w[i]*gusb[n]))
								  *(w[j]*(a[0]+a[2]*gusb[n])*(1+u[j]*gusa[m])
								   -u[j]*(b[1]+b[2]*gusa[m])*(1+w[j]*gusb[n])))/(16*jac);
				   }
		 }
     double s=0.0;
     if(sbound[NOE-1][0]==0)
	 {
		 for(i=0;i<4;i++)
               for(j=0;j<4;j++)
			   {
				   g[i][j]=0;
			   }
	 }
	 else
	 {
         s=(v/E)*sqrt(pow((CN[sbound[NOE-1][1]-1][0]-CN[sbound[NOE-1][2]-1][0]),2)+pow((CN[sbound[NOE-1][1]-1][1]-CN[sbound[NOE-1][2]-1][1]),2));
		 for(i=0;i<4;i++)
		 {
			 if(ngn[i]==sbound[NOE-1][1])
			 {
				 g[i][i]=1.0/3.0;
			 }
			 else
			 {
				 if(ngn[i]==sbound[NOE-1][2])
				 {
					 g[i][i]=1.0/3.0;
				 }
				 else
				 {
					 g[i][i]=0.0;
				 }
			 }
		 }
         for(i=0;i<4;i++)
		 {
		     if(g[i][i]!=0)
			 {
			     for(j=0;j<4;j++)
				    g[j][i]=1.0/6.0;
				    g[i][j]=1.0/6.0;
				    g[i][i]=1.0/3.0;
			 }
		 }
	 }
	 for(i=0;i<4;i++)
		 for(j=0;j<4;j++)
		 {
			 g[i][j]=s*g[i][j];
		 }
}

void rmatrix(int NOE,int NGN[][4],float CN[][2],
	             double r[][4],float E,float v,float aa)
				
{
     int i=0,j=0,k=0,m=0,n=0;
	 int ngn[4];
         for(i=0;i<4;i++)ngn[i]=NGN[NOE-1][i];
     float cn[4][2];
         for(i=0;i<4;i++)
             for(j=0;j<2;j++)
                   cn[i][j]=CN[ngn[i]-1][j];
	 float a[3]={0,0,0};
	 float b[3]={0,0,0};
	 for( k=0;k<4;k++)
		  for(i=0;i<4;i++)
		  {
			  r[k][i]=0;
		  }
	 a[0]=((-1)*cn[0][0]+cn[1][0]+cn[2][0]-cn[3][0])/4;
     a[1]=((-1)*cn[0][0]-cn[1][0]+cn[2][0]+cn[3][0])/4;
     a[2]=(cn[0][0]-cn[1][0]+cn[2][0]-cn[3][0])/4;
     b[0]=((-1)*cn[0][1]+cn[1][1]+cn[2][1]-cn[3][1])/4;
     b[1]=((-1)*cn[0][1]-cn[1][1]+cn[2][1]+cn[3][1])/4;
     b[2]=(cn[0][1]-cn[1][1]+cn[2][1]-cn[3][1])/4;
	 double gusa[2]={-0.5773503,0.5773503};
     double gusb[2]={-0.5773503,0.5773503};
     double jac=0;
	 float u[4]={-1,1,1,-1};
	 float w[4]={-1,-1,1,1};
     for(m=0;m<2;m++)
	     for(n=0;n<2;n++) 
		 {
			 jac=(a[0]+a[2]*gusb[n])*(b[1]+b[2]*gusa[m])-(a[1]+a[2]*gusa[m])*(b[0]+b[2]*gusb[n]);
			 for(i=0;i<4;i++)
			       for(j=0;j<4;j++)
				   {
					   r[i][j]+=((1+u[i]*gusa[m])*(1+w[i]*gusb[n])*(1+u[j]*gusa[m])*(1+w[j]*gusb[n])*jac)/(aa*16);
				   }
         }
}

void adjust(double GK[][NN],int fbound[NN],float Tb,double FR[NN],double temper[NN])
{
     int i;
	 for(i=0;i<NN;i++)
	 {
		 temper[i]=FR[i];
	 }
	 for(i=0;i<NN;i++)
	 {
		 if(fbound[i]==1)
		 {
			 GK[i][i]=1.0E10;//大数法
			 temper[i]=Tb*1.0E10;
		 }
	 }
}
void gauss(double GK[NN][NN],double temper[NN])
{
      int i,j,k;
      double temp=0.0;
      double temp1=0.0;
      for(i=0;i<NN;i++)
      {
           temp=GK[i][i];
           for(j=0;j<NN;j++)
           {   GK[i][j]/=temp;
           }
           temper[i]/=temp;
           for(j=i+1;j<NN;j++)
           {
               temp1=GK[j][i];
               for(k=i;k<NN;k++)
               {
                    GK[j][k]-=temp1*GK[i][k];
               }
               temper[j]-=temper[i]*temp1;
           }
      }
      for(i=NN-1;i>=0;i--)
      {
            for(j=i-1;j>=0;j--)
            {
                  temper[j]=temper[j]-GK[j][i]*temper[i];
                  GK[j][i]=0;
            }
      }
}
void temperature(double H[][NN],double R[][NN],double F[NN],double P[NN],
				 float T0,float Ta,float Tb,int fbound[NN],float time,
				 int d,double temper[NN])
{
	int i,j,k;
	double GK[NN][NN];
	for(i=0;i<NN;i++)
		for(j=0;j<NN;j++)
		{
			GK[i][j]=0.0;
		}
	double FR[NN];
	for(i=0;i<NN;i++)
	{
		FR[i]=0.0;
	}
	float dt=0;
    double Tn[NN];
	for(i=0;i<NN;i++)
		{
			Tn[i]=T0;
		}
	 dt=time/d;
     for(i=1;i<=d;i++)
	 {
		 for(j=0;j<NN;j++)
			 for(k=0;k<NN;k++)
			 {
				 GK[j][k]=H[j][k]+R[j][k]/dt;
			 }
		 for(j=0;j<NN;j++)
		 {
			 FR[j]=(F[j]*30/pow(1+dt*i,2)+P[j]*Ta);
			 for(k=0;k<NN;k++)
			 {
                   FR[j]+=R[j][k]*Tn[k]/dt;
			 }
		 }
         adjust(GK,fbound,Tb,FR,temper);
         gauss(GK,temper);
		 for(j=0;j<NN;j++)
		 {
			 Tn[j]=temper[j];
		 }
	 }
	 ofstream outfile;
	 outfile.open("temperature.txt");
	 outfile<<fixed<<showpoint<<setprecision(2);
	 for(i=0;i<NN;i++)
	 {
		 outfile<<temper[i]<<endl;
	 }
	 outfile.close();
}