emap1.c

一个求解三维时变场的有限元算法,C语言实现,根据MAILE提出的算法实现.

/*                ***********NOTICE**********                              */
/*                                                                         */
/*  The EMAP finite element modeling codes were created at the             */
/*  University of Missouri-Rolla Electromagnetic Compatibility Laboratory. */
/*  They are intended for use in teaching or research.  They may be freely */
/*  copied and distributed PROVIDED THE CODE IS NOT MODIFIED IN ANY MANNER.*/
/*                                                                         */
/*  Suggested modifications or questions about these codes can be          */
/*  directed to Dr. Todd Hubing, Department of Electrical Engineering      */
/*  University of Missouri-Rolla, Rolla, MO  65401.  Principal authors     */
/*  of these codes are Mr. Mohammad Ali and Mr. Girish Bhat of the         */
/*  University of Missouri-Rolla.                                          */
/*                                                                         */
/*  Neither the authors nor the University of Missouri makes any warranty, */
/*  express or implied, or assumes any legal responsibility for the        */
/*  accuracy, completeness or usefulness of these codes or any information */
/*  distributed with these codes.                                          */
/*                                                                         */
/*  2/2/95                                                                 */

/***************************************************************************
*   PROGRAM     :  Emap1.c (Version 2.02)
*   LAST UPDATE :  February 15, 1995
*   DESCRIPTION :  A 3-D Finite Element Modeling Code for Analyzing
                   Time Varying Complex Electromagnetic fields.
                   (a scalar code using a variational approach)
*   INPUT FILE  :  emap1.in or or 1st argument of emap command
*   OUTPUT FILE :  output filenames are specified by the input file
*
****************************************************************************/

/**************************** Include Files ********************************/

#include<stdio.h>
#include<signal.h>
#include<string.h>
#include<stdlib.h>
#include<sys/time.h>
#include<math.h>

/******************** Dynamic Memory Allocation Functions ******************/

int *INT_Vector(nh)
/* Allocates an int vector with range [0 ... nh-1]  */
int nh;
{
int *v;
v=(int *)malloc((unsigned) nh*sizeof(int));
return v;
}
double *FLOAT_Vector(nh)
/* Allocates a double vector with range [0 ... nh-1] */
int nh;
{
double *v;
v=(double *)malloc((unsigned) nh*sizeof(double));
return v;
}

double **FLOAT_Matrix(row,column)
/*Allocates a double matrix with range [0..rw-1][0..cl-1]*/
int row,column;
{
int i;double **m;
m=(double **) malloc((unsigned) row*sizeof(double*));
for(i=0;i<row;i++)
{m[i]=(double *) malloc((unsigned) column*sizeof(double));}
return m;
}
int **INT_Matrix(row,column)
/*Allocates an int matrix with range [0..rw-1][0..cl-1]*/
int row,column;
{
int i;int **m;
m=(int **) malloc((unsigned) row*sizeof(int));
for(i=0;i<row;i++)
{m[i]=(int *) malloc((unsigned) column*sizeof(int));}
return m;
}

/************************* Default Parameters ******************************/

#define  GBL_Matrix_MaxColNos 43
#define  FRCD_Matrix_MaxColNos 10

/*********************** Function Prototypes *******************************/

void ParameterInfo();
void Assign_Global_Coordinates();
void Assign_hexahedron_NodeNumbering();
void Find_CoFactor();
void Find_Volume_Determinant();
void HexaHedra_SubDivision_1();
void HexaHedra_SubDivision_2();
void TetraHedron_SubMatrix();
void HexaHedron_Matrix();
void HexaHedron_link_1();
void HexaHedron_link_2();
void Find_Global_Matrix();
void Read_Input_Pass_1();
void Read_Input_Pass_2();
void Fill_Vectors();
void Partition_Global_Matrix();
void Find_RHS_Vector();
void Band_Matrix_Solver();
void Produce_Output();
void EfieldatNode();
void Count_Half_BandWidth();
void swap();
/***************************************************************************/

double Determinant_TetraMatrix();
int Locate_NZ_Element_Column();
int Srch_NZ_Element_Column();

/*********************** Fixed Matrix Pointers *****************************/

double cofs[4][4];
double Volume_Cord[4][4],VVolume_Cord[4][4];
int   tetrnod[5][4],tetlink[5][4],connect1[8][6];
double Tet_Matrix_Element[5][12][12],Hex_Matrix_Element[24][24];

/************************** Global Pointers ********************************/

double *REL_Permittivity,*FORC_Val,**Dat_Buffer,**GBL_Matrix_Data;
double **MeshNodeCord,**FORC_Matrix_Data,*RHS_Vector,**HalfBand_Matrix;
int **GBL_Matrix_ColNos,*GBL_Matrix_Index;
int **FORC_Matrix_ColNos,*UN_Vector,*FORC_Vector,**HexHed;
int *Loc_UN_Node,*Loc_UN_NodeCord;
FILE *InF, *OutF_0, *OutF_1, *OutF_2, *OutF_3;

/************************** Global Variables *******************************/

int Dimension_X,Dimension_Y,Dimension_Z;
double OperateFreq,CellDimension;
double DivisorX,DivisorY,DivisorZ;
int Min_X=9999,Min_Y=9999,Min_Z=9999,Max_X=-9999,Max_Y=-9999,Max_Z=-9999;
double det=0.0, WaveNumber, freq;
int t, cn=0;
int TotNum_GBL_Nodes,TotNum_HexaHedron;
struct timeval;
int Half_BandWidth=0,TotNum_Unknown_Var=0,TotNum_Known_Var=0;
char Out_FileName0[20],Out_FileName1[20],Out_FileName2[20],Out_FileName3[20];

/************************** Main Function **********************************/

main(argc, argv)

int argc;
char *argv[];

{
int    i,j,GBL_Matrix_RowNos;
char   Ifile[20];
        if (argc > 2) {
        fprintf(stderr,"Usage: Emap1 [input_file] \n"); exit(1);}
        if (argc < 2) {  argv[1] = "emap1.in";   strcpy(Ifile, argv[1]); }
        else                                    strcpy(Ifile, argv[1]);
        fprintf(stdout, "\nEMAP1 Version 2.02\n\n");
        fprintf(stdout, "The input will be read from the file:  %s\n", Ifile);
        InF = fopen(argv[1], "r");
        if (InF == NULL) {
        fprintf(stderr, " Error: Can't open input file.\n"); exit(1); }
        
Read_Input_Pass_1();
ParameterInfo();

/* Allocates Dynamic memory */
HexHed=INT_Matrix(TotNum_HexaHedron,8);
MeshNodeCord=FLOAT_Matrix(TotNum_GBL_Nodes,3);
REL_Permittivity=FLOAT_Vector(TotNum_GBL_Nodes);
Dat_Buffer=FLOAT_Matrix(TotNum_GBL_Nodes,3);

Read_Input_Pass_2();
fprintf(stderr,"TotNum_Unknown_Var=%d\n", TotNum_Unknown_Var);
fprintf(stderr,"TotNum_Known_Var=%d\n", TotNum_Known_Var);

Assign_Global_Coordinates();
Assign_hexahedron_NodeNumbering();

/* Allocates Dynamic memory */
 UN_Vector=INT_Vector(TotNum_Unknown_Var);
 FORC_Vector=INT_Vector(TotNum_Known_Var);
 FORC_Val=FLOAT_Vector(TotNum_Known_Var);
 Loc_UN_Node=INT_Vector(TotNum_Unknown_Var);
 Loc_UN_NodeCord=INT_Vector(TotNum_Unknown_Var);
 FORC_Matrix_Data=FLOAT_Matrix(TotNum_Unknown_Var,FRCD_Matrix_MaxColNos);
 FORC_Matrix_ColNos=INT_Matrix(TotNum_Unknown_Var,FRCD_Matrix_MaxColNos);
 RHS_Vector=FLOAT_Vector(TotNum_Unknown_Var);
 GBL_Matrix_RowNos=TotNum_GBL_Nodes*3;
 GBL_Matrix_Data=FLOAT_Matrix(GBL_Matrix_RowNos,GBL_Matrix_MaxColNos);
 GBL_Matrix_ColNos=INT_Matrix(GBL_Matrix_RowNos,GBL_Matrix_MaxColNos);
 GBL_Matrix_Index=INT_Vector(GBL_Matrix_RowNos);


/*  Initialize Dynamically Allocated Variables  */

for(i=0;i<=3*TotNum_GBL_Nodes-1;i++) {
  GBL_Matrix_Index[i]=0;
  for(j=0;j<=GBL_Matrix_MaxColNos-1;j++) {
    GBL_Matrix_Data[i][j]=0.0;
    GBL_Matrix_ColNos[i][j]=0; }}
   for(i=0;i<=23;i++) for(j=0;j<=23;j++)
   { Hex_Matrix_Element[i][j]=0.0;}
   for(i=0;i<=11;i++)
    for(j=0;j<=11;j++)
      for(t=0;t<=4;t++)
       { Tet_Matrix_Element[t][i][j]=0.0;}

Fill_Vectors();

for(cn=0;cn<=(TotNum_HexaHedron-1);cn++) {      
 if(cn%2 ==1) HexaHedra_SubDivision_2();
 else         HexaHedra_SubDivision_1();

   for(i=0;i<=11;i++)
    for(j=0;j<=11;j++)
      for(t=0;t<=4;t++)
       { Tet_Matrix_Element[t][i][j]=0.0;}
  for(t=0;t<=4;t++)          {
    Find_Volume_Determinant();         
    Find_CoFactor();          
    TetraHedron_SubMatrix(); }

    HexaHedron_Matrix(); 

   if(cn%2 ==1) HexaHedron_link_2();
   else         HexaHedron_link_1();
   Find_Global_Matrix();

    } /* end of cn loop */

Count_Half_BandWidth();

/* Allocates Dynamic memory */
HalfBand_Matrix=FLOAT_Matrix(TotNum_Unknown_Var+1,Half_BandWidth+1);

for(i=1;i<=TotNum_Unknown_Var;i++)
 for(j=1;j<=Half_BandWidth;j++)
 {HalfBand_Matrix[i][j] =0.0;}

Partition_Global_Matrix();
Find_RHS_Vector();

Band_Matrix_Solver();
Produce_Output();
return(0);
}

/*********************End of Main Program **********************************/
/***************************************************************************
*
*   FUNCTION        : ParameterInfo()
*   DESCRIPTION     : Initializes global variables. 
*
****************************************************************************/
void ParameterInfo()
{
double FreeSpaceVel,AbsPermeable,AbsPermitt,WaveLength, UCellLen;
 
 AbsPermeable  = 1.25663706144E-06;
 AbsPermitt    = 8.8542E-12;
 FreeSpaceVel  = 1.0/sqrt(AbsPermeable*AbsPermitt);
 WaveLength = FreeSpaceVel/(OperateFreq*1.0E+06);
 WaveNumber = 2.0*M_PI/WaveLength;

 Dimension_X=Max_X-Min_X; Dimension_Y=Max_Y-Min_Y; Dimension_Z=Max_Z-Min_Z;
 TotNum_HexaHedron=Dimension_X*Dimension_Y*Dimension_Z;
 TotNum_GBL_Nodes=(Dimension_X+1)*(Dimension_Y+1)*(Dimension_Z+1);
 DivisorX = CellDimension;
 DivisorY = CellDimension;
 DivisorZ = CellDimension;
 UCellLen = CellDimension*1000.0;
 
 fprintf(stderr,"Free Space Velocity     =%8.2E m/s\n",FreeSpaceVel );
 fprintf(stderr, "Operating frequency    = %g MHz \n", OperateFreq);
 fprintf(stderr, "Free Space Wave Number = %8.2f\n",WaveNumber);
 fprintf(stderr, "Free Space WaveLength  = %8.2f mm\n",WaveLength*1000.0);
 fprintf(stderr, "Unit cell dimension    = %10.4f  mm\n\n", UCellLen);
 fprintf(stderr, "Dimension of PEC box       = %g mm x %g mm x %g mm \n",
 Dimension_X*UCellLen,Dimension_Y*UCellLen,Dimension_Z*UCellLen);
 fprintf(stderr, "Total Number of cube cell  = %d  \n",TotNum_HexaHedron);
 fprintf(stderr, "Total Number of Nodes      = %d  \n",TotNum_GBL_Nodes);
}

/***************************************************************************
*
*   FUNCTION        : Assign_Global_Coordinates()
*   DESCRIPTION     : Assigns Global coordinates for each meshnode. 
*
****************************************************************************/

void Assign_Global_Coordinates()
 {
  int i,j,k,Xcount,Ycount,Zcount,Count_h=0; 
  double Xtemp,Ytemp,Ztemp;
   {
    Xcount=(int)Dimension_X;Ycount=(int)Dimension_Y;Zcount=(int)Dimension_Z;
    Xtemp=Min_X;  Ytemp=Min_Y;  Ztemp=Min_Z;
     for(k=0;k<=Zcount;k++) {
      for(j=0;j<=Ycount;j++) {
       for(i=0;i<=Xcount;i++) {
          MeshNodeCord[Count_h][0]=Xtemp; 
          MeshNodeCord[Count_h][1]=Ytemp; 
          MeshNodeCord[Count_h][2]=Ztemp; 
          Xtemp=Xtemp+DivisorX; Count_h++; }
          Xtemp=Min_X; Ytemp=Ytemp+DivisorY; }
          Ytemp=Min_Y; Ztemp=Ztemp+DivisorZ; }
    }
 }

/***************************************************************************
*
*   FUNCTION        : Assign_hexahedron_NodeNumbering()
*   DESCRIPTION     : Assigns Global nodenumbering for each meshnode.
*
****************************************************************************/

void Assign_hexahedron_NodeNumbering()
  {
   int i,start=1;
   {
   for(i=1;i<=TotNum_HexaHedron;i++)
     {
      HexHed[i-1][0]=start;
      HexHed[i-1][1]=HexHed[i-1][0]+1;
      HexHed[i-1][2]=HexHed[i-1][0]+Dimension_X+1;
      HexHed[i-1][3]=HexHed[i-1][2]+1;
      HexHed[i-1][4]=HexHed[i-1][0]+(Dimension_X+1)*(Dimension_Y+1);
      HexHed[i-1][5]=HexHed[i-1][4]+1;
      HexHed[i-1][6]=HexHed[i-1][4]+(Dimension_X+1);
      HexHed[i-1][7]=HexHed[i-1][6]+1;

      if((i%(Dimension_Y*Dimension_X)==0))
            { start=start+(Dimension_X+3); }
      else  { start=start+1; }

      if((start%(Dimension_X+1))==0)
           { start++; }

     }
  }
 }

/***************************************************************************
*
*   FUNCTION        : HexaHedra_SubDivision_1()
*   DESCRIPTION     : Divide each Hexahedron into five Tetrahedron.
*                     Find links between them.  (Method 1)
*
****************************************************************************/

void HexaHedra_SubDivision_1()
    {
      tetrnod[0][0]=HexHed[cn][0];   tetlink[0][0]=0;
      tetrnod[0][1]=HexHed[cn][2];   tetlink[0][1]=2;
      tetrnod[0][2]=HexHed[cn][4];   tetlink[0][2]=4;
      tetrnod[0][3]=HexHed[cn][1];   tetlink[0][3]=1;

      tetrnod[1][0]=HexHed[cn][5];   tetlink[1][0]=5;
      tetrnod[1][1]=HexHed[cn][1];   tetlink[1][1]=1;
      tetrnod[1][2]=HexHed[cn][4];   tetlink[1][2]=4;
      tetrnod[1][3]=HexHed[cn][7];   tetlink[1][3]=7;

      tetrnod[2][0]=HexHed[cn][3];   tetlink[2][0]=3;
      tetrnod[2][1]=HexHed[cn][1];   tetlink[2][1]=1;
      tetrnod[2][2]=HexHed[cn][7];   tetlink[2][2]=7;
      tetrnod[2][3]=HexHed[cn][2];   tetlink[2][3]=2;

      tetrnod[3][0]=HexHed[cn][6];   tetlink[3][0]=6;
      tetrnod[3][1]=HexHed[cn][2];   tetlink[3][1]=2;
      tetrnod[3][2]=HexHed[cn][7];   tetlink[3][2]=7;
      tetrnod[3][3]=HexHed[cn][4];   tetlink[3][3]=4;

      tetrnod[4][0]=HexHed[cn][1];   tetlink[4][0]=1;
      tetrnod[4][1]=HexHed[cn][2];   tetlink[4][1]=2;
      tetrnod[4][2]=HexHed[cn][4];   tetlink[4][2]=4;
      tetrnod[4][3]=HexHed[cn][7];   tetlink[4][3]=7;
   }

/***************************************************************************
*
*   FUNCTION        : HexaHedra_SubDivision_2()
*   DESCRIPTION     : Divide each Hexahedron into five Tetrahedron.
*                     Find links between them. (Method 2)
*
****************************************************************************/

void HexaHedra_SubDivision_2()
    {
      tetrnod[0][0]=HexHed[cn][2];   tetlink[0][0]=2;
      tetrnod[0][1]=HexHed[cn][0];   tetlink[0][1]=0;
      tetrnod[0][2]=HexHed[cn][3];   tetlink[0][2]=3;
      tetrnod[0][3]=HexHed[cn][6];   tetlink[0][3]=6;
 
      tetrnod[1][0]=HexHed[cn][7];   tetlink[1][0]=7;
      tetrnod[1][1]=HexHed[cn][6];   tetlink[1][1]=6;
      tetrnod[1][2]=HexHed[cn][3];   tetlink[1][2]=3;
      tetrnod[1][3]=HexHed[cn][5];   tetlink[1][3]=5;