emap3.c

计算电磁学的３维矢量有限元程序 EMAP-3 is a vector (edge element) code. Vector codes are generally not affected

/*                ***********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     :  Emap3.c (Version 2.02)
*   LAST UPDATE :  February 15, 1995
*   DESCRIPTION :  A 3-D Vector Finite Element Modeling Code for Analyzing
                   Time Varying Complex Electromagnetic Fields.
*   INPUT FILE  :  emap3.in or or 1st argument of emap command
*   OUTPUT FILE :  file(s) specified by 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  19
#define  FRCD_Matrix_MaxColNos 10

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

void Read_Input_Pass_1();
void ParameterInfo();
void Read_Input_Pass_2();
void AssignGlobalCoorDinates();
void AssignHexHedronNodeNumbering();
void AssignHexHedronEdgeNumbering();
void HexHedraSubDiv1();
void HexHedraSubDiv2();
void EdgeSubDiv1();
void EdgeSubDiv2();
void FindTetraHedronVolume();
void ComputeTetraCoFactor();
void S_Matrix();
void T_Matrix();
void TetraSubMatrix();
void FindGlobalMatrix();
void CountHalfBandWidth();
void PartitionGlobalMatrix();
void ComputeRHSVector();
void BandMatrixSolver();
void GlobalEdgeEndsDiv1();
void GlobalEdgeEndsDiv2();
void EfieldatNode();
void Produce_Output();
FILE *InF, *OutF_0, *OutF_1, *OutF_2, *OutF_3;

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

int
IEdge,
edge_end[6][2],
**HexNode,
**HexEdgeNum,
**TetGlobalNodeNum,
**TetEdgeNum,
**GBL_Matrix_ColNos,
*GBL_Matrix_Index,
*ForceStat,
 **xn,
**FORC_Matrix_ColNos,
*InnerEdgeStat,
*BoundEdgeStat,
*ForcdEdgeStat,
**TetGlobalEdgeEnds;
 
double
DivisorX,
DivisorY,
DivisorZ,
CellDimension,
OperateFreq,
WaveNumber,
TetVolume,
A_mat[5][6][6],
S_mat[6][6],
T_mat[6][6],
CoFactor[4][4],
**GBL_Matrix_Data,
**HalfBandMatrix,
**FORC_Matrix_Data,
*RHSVector,
**ForcdValue,
**RELPerm,
** NodeCord,
*Sigma,
*Epsilon,
*EdgeStatus,
*EfieldData;

char Out_FileName0[20];
char Out_FileName1[20];
char Out_FileName2[20];
char Out_FileName3[20];

int 
TotBoundNodeNum,
TotBoundEdgeNum,
TotForcdEdgeNum,
TotTrngleNum,
TotInnerEdgeNum,
start,TotNumHexHedron,
HalfBandWidth,
t,
XdiM,YdiM,ZdiM,
Max_X=0,
Max_Y=0,
Max_Z=0,
Min_X=0,
Min_Y=0,
Min_Z=0,
TotEdgeNum=0,
TotNodeNum=0,
HexNum,
TotTetElement;

/***************************************************************************
*  Utility Functions
****************************************************************************/

 void VTXsub1(Count_i,Count_k,Buff,Buff1)
 int Count_i,Count_k;
 double Buff[4][3],Buff1[3];
 {
 int Count_j;
 for(Count_j=0;Count_j<=2;Count_j++)
 Buff1[Count_j] = Buff[Count_i][Count_j]-Buff[Count_k][Count_j];
 }
 void VTXcross(Buff1,Buff2,Buff)
 double Buff1[3],Buff2[3],Buff[3];
 {
 Buff[0] = Buff1[1]*Buff2[2] - Buff2[1]*Buff1[2];
 Buff[1] = Buff1[2]*Buff2[0] - Buff2[2]*Buff1[0];
 Buff[2] = Buff1[0]*Buff2[1] - Buff2[0]*Buff1[1];
 }
 double Sign(Value)
 int Value;
 {
 double Value1,Value2;
 Value1 = (double)Value ; Value2 = Value1 / fabs(Value1) ;
 return Value2;
 }
 double VTXmag(Buff1,Buff2)
 double Buff1[3],Buff2[3];
 {
 double Value,ValueX,ValueY,ValueZ;
 ValueX = Buff1[0] - Buff2[0];
 ValueY = Buff1[1] - Buff2[1];
 ValueZ = Buff1[2] - Buff2[2];
 Value  = sqrt(ValueX*ValueX + ValueY*ValueY + ValueZ*ValueZ);
 return Value;
 }
   double ff(i,j) int i,j; {
   double Prod;
   Prod = CoFactor[1][i-1] * CoFactor[1][j-1] + CoFactor[2][i-1] * CoFactor[2][j-1]
        + CoFactor[3][i-1] * CoFactor[3][j-1]  ;
              return Prod;}
 int Srch_NZ_Element_Column(RowNum,ColNum)
   {
   int Count_i; {
   for(Count_i=0;Count_i<=GBL_Matrix_Index[RowNum]-1;Count_i++)
    {
    if(ColNum==(GBL_Matrix_ColNos[RowNum][Count_i]))
     {return Count_i;}
     } return -1;       }
   }
   int Locate_NZ_Element_Column(RowNum,ColNum)
     {
      int Count_i; {
      for(Count_i=0;Count_i<=FRCD_Matrix_MaxColNos-1;Count_i++)
       {
       if(ColNum==(FORC_Matrix_ColNos[RowNum][Count_i]))
        {return Count_i;}
       } return -1;        }
     }

/**************************** MAIN PROGRAM *********************************/

main(argc, argv)

int argc;
char *argv[];

{
int con,i,j,k;
char   Ifile[20];
        if (argc > 2) {
        fprintf(stderr,"Usage: Emap3 [input_file] \n"); exit(1);}
        if (argc < 2) {  argv[1] = "emap3.in";   strcpy(Ifile, argv[1]); }
        else                                    strcpy(Ifile, argv[1]);
        fprintf(stdout, "\nEMAP3 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();
HexNode = INT_Matrix(TotNumHexHedron,8);
HexEdgeNum = INT_Matrix(TotNumHexHedron,19);
RELPerm     = FLOAT_Matrix(TotNumHexHedron,2);
NodeCord = FLOAT_Matrix(TotNodeNum,3);
TetGlobalNodeNum  = INT_Matrix(5,4);
TetEdgeNum  = INT_Matrix(5,6);
TetGlobalEdgeEnds = INT_Matrix(TotEdgeNum,2);

Read_Input_Pass_2();
AssignGlobalCoorDinates();
AssignHexHedronNodeNumbering();
AssignHexHedronEdgeNumbering();

GBL_Matrix_Data=FLOAT_Matrix(TotEdgeNum,GBL_Matrix_MaxColNos);
GBL_Matrix_ColNos=INT_Matrix(TotEdgeNum,GBL_Matrix_MaxColNos);
GBL_Matrix_Index=INT_Vector(TotEdgeNum);

for(i=0;i<TotEdgeNum;i++)  { GBL_Matrix_Index[i] = 0;
for(j=0;j<GBL_Matrix_MaxColNos;j++)
{GBL_Matrix_Data[i][j] = 0.0; GBL_Matrix_ColNos[i][j] = 0;}
                           }

 for(k=0;k<=ZdiM-1;k++)
 for(j=0;j<=YdiM-1;j++)
 for(i=0;i<=XdiM-1;i++)   {