emap4.c

计算电磁学的３维矢量有限元方法通用的原程序

/*                ***********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.                                          */
/*                                                                         */
/* 5/98                                                                    */
/***************************************************************************
*   PROGRAM     :  Emap4.c (Version 1.0)
*   LAST UPDATE :  May 1998 
*   DESCRIPTION :  A 3-D Vector Finite Element Modeling Code for Analyzing
                   Time Varying Complex Electromagnetic Fields.
*   LATEST DEV. :  Inclusion of Isource, modeling of coax cables.
*   INPUT FILE  :  emap4.in or or 1st argument of emap command
*   OUTPUT FILE :  file(s) specified by input file
*   USER'S GUIDE:  http://www.emclab.umr.edu/emap4
*                  (contains all the relevant references)
*
****************************************************************************/


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

#include<stdio.h>
#include<signal.h>
#include<string.h>
#include<stdlib.h>
#include<sys/time.h>
#include<math.h>
#include"complex.h" /* Complex routines */
#include"alloc.h"  /* Memory allocation routines */

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

#define  GBL_Matrix_MaxColNos  19  /* The FEM is sparse and any row has a 
				      maximum of 19 columns */
#define  FRCD_Matrix_MaxColNos 11  
#define  MAX_ITERATIONS 25000      /* The maximum number of iterations for 
				      the CBCG solver */
#define  TOLERANCE  0.00001        /* The amount of accuracy required in the 
				      code */

/************************** 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 ConjugateSolver();
void GlobalEdgeEndsDiv1();
void GlobalEdgeEndsDiv2();
void EfieldatNode();
void Produce_Output();

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

/* The following are the file variables */

FILE *InF, *OutF_0, *OutF_1, *OutF_2, *OutF_3;

int 

TotBoundNodeNum, /* The number of boundary nodes */
TotBoundEdgeNum, /* The number of unforced boundary edges */
TotForcdEdgeNum, /* The number of forced edges */
TotTrngleNum, /* The total number of triangles on the surface (not used) */
TotInnerEdgeNum, /* The number of unforced or "inner" edges */
start, /* Dummy variable used to define the start in edge and 
	  node numbering schemes */
TotNumHexHedron, /* The total number of hexahedra */
t, /* Dummy variable used to specify the tetrahedron number */
XdiM,YdiM,ZdiM, /* The dimensions of the meshed area */
Max_X=0, /* Start and end coordinates of the meshed area */
Max_Y=0, /* The mesh goes from (Min_X,Min_Y,Min_Z) to (Max_X,Max_Y,Max_Z) */
Max_Z=0,
Min_X=0,
Min_Y=0,
Min_Z=0,
TotEdgeNum=0, /* Total number of edges in the mesh */
TotNodeNum=0, /* Total number of nodes in the mesh */
HexNum, /* Index Variable used to specify the hexahedron number */
TotTetElement; /* Total number of tetrahedral elements in the mesh */

int freq_step = 0; /* The variable to indicate the number of the 
		      frequency step in progress */

int

  IEdge, /* Dummy edge index */
  edge_end[6][2], /* temporary variable while calculating the edge-ends */
**HexNode,  /* Node numbers : 1st index is the hexahedron number
	       The second index is the local hexahedron node number */
**HexEdgeNum, /* Edge numbers : 1st index is the hexahedron number
	       The second index is the local hexahedron edge number */
**TetGlobalNodeNum, /* This is a [5:4] matrix. 
		       Node numbers : 1st index is the tetrahedron number [0:4]
		       The second index is the local tetrahedron node number */
**TetEdgeNum, /* Edge numbers : 1st index is the tetrahedron number
	       The second index is the local tetrahedron edge number */
**TetGlobalEdgeEnds, /* Edge Ends : 1st index is the gloal edge number
	       The second index is 0 (beginning) or 1 (end) */
*ForceStat, /* The flag which indicates whether the edge is forced (1) or
	       an unforced edge (0) */
*InnerEdgeStat, /* Contains TotalInnerEdgeNum entries containing the 
		   global edge number of each inner edge */
*BoundEdgeStat, /* Contains TotalBoundEdgeNum entries containing the 
		   global edge number of each unforced boundary edge */
*ForcdEdgeStat, /* Contains TotalForcdEdgeNum entries containing the 
		   global edge number of each forced edge */
*InnerEdgeStat1,/* Contains TotEdgeNum entries containing the 
		   inner edge number of each global edge if applicable
		   otherwise is zero */
*BoundEdgeStat1,/* Contains TotEdgeNum entries containing the 
		   boundary  edge number of each global edge if applicable
		   otherwise is zero */
*ForcdEdgeStat1,/* Contains TotEdgeNum entries containing the 
		   forced edge number of each global edge if applicable
		   otherwise is zero */
/* The above three variables exist for ease and simplicity of the program */

/* The following variables store the Global and the 
   actual matrix to be solved 
   (see chapter on sparse matrices in the user's guide) */
**GBL_Matrix_ColNos,
**LHS_Matrix_ColNos,
**FORC_Matrix_ColNos,
*GBL_Matrix_Index,
*LHS_Matrix_Index;

complex
**GBL_Matrix_Data,
**LHS_Matrix_Data;

complex

A_mat[5][6][6], /* Temp. variable to store the A_mat for each tetrahedron */
S_mat[6][6],/* Temp. variable to store the S_mat for each tetrahedron */
T_mat1[6][6],/* Temp. variable to store the T_mat for each tetrahedron */
*Resist, /* Variable to store the terms arising due to resistors */
*Isource, /* Variable to store the terms arising due to current sources */
*RHSVector, /* The RHS vector in the final equation which is solved */
*ForcdValue, /* The value of the electric fields along the "forced edges */
*RELPerm, /* The relative permittivity of each hexahedron */
*EfieldData; /* The value of the electric field along every edge in the mesh */




int
**xn,
NUM_OF_STEPS=1; /* The number of frequency steps that this code is run */

 
double

*DivisorX, /* This variable stores the value of cell dimension 
	      (length of each brick in the x-direction */
*DivisorY, /* This variable stores the value of cell dimension 
	      (length of each brick in the y-direction */
*DivisorZ, /* This variable stores the value of cell dimension 
	      (length of each brick in the z-direction */
OperateFreq, /* Frequency of operation 
		(the code can run only for one frequency at a time */
FreeSpaceVel, /* "c" */
AbsPermeable, 
AbsPermitt,
WaveLength, /* the free space wavelength */
WaveNumber, /* the free-space wavenumber */
Omega, /* "w"=2*pi*"operatefreq" */
TetVolume, /* Dummy variable to store the volume of each tetrahedron */
CoFactor[4][4], /* Variable which stores the "co-factor matrix" 
		   (the a's, b;s c's and d's.
		   for each tetrahedron (see reference [9] in the 
		   user's guide on the web "http://www.emclab.umr.edu/emap4*/
**NodeCord, /* Stores the node coordinates for each global node 
	       the first index is the global node number, the second index
	       is 0 (x-coord), 1 (y-coord) or 2 (z_coord) */
*Sigma, /* Stores the value of the conductivity of each brick (hexahedron)
	   if applicable, otherwise stores 0 */
*Epsilon, /* Stores the value of the relative permittivity of each brick 
	     (hexahedron) if applicable, otherwise stores 0 */
*EdgeStatus, /* This variable is a vector of length [TotEdgeNum]. 
		If the edge is an inner edge, it stores -1, 
		if it is an unforced boundary edge, it stores -2 
		and if it is a forced edge, it stores the 
		magnitude specified in the input file */
*EdgeStatus1, /* This variable is the same as EdgeStatus and 
		 stores the phase specified in the input file */
FREQ_INC;

complex   /* The variables used in */ 
*X,       /* the conjugate solver  */
*P,       /* routine (see the sparse matrix chapter in the 
	     user's guide for more 
	     information about the variables */
*PP,      
*R,
*RR,
*A_P,
*AH_PP,
Alpha,
Beta,
temp_var;

complex 
*Old_Solution; /* Used in the conjugate solver to feed the previous solution
		  as the starting vector for the next frequency in a 
		  swept frequency simulation */


complex
a_PML,b_PML,c_PML, /* THe values of the PML variables a,b,c inside the specific
		      hexahedron for which the matrix is being computed*/
*a_PML_vec,*b_PML_vec,*c_PML_vec; /* THe values of the PML variables a,b,c 
				     inside all the hexahedra in the PML*/
int *Is_PML; /* Variable which asks the question 
		"Is this hexahedron part of a PML layer?" */ 

char 
Ifile[20], /* Input filename */ 
Out_FileName0[20], /* Output filenames */
Out_FileName1[20], 
Out_FileName2[20], 
Out_FileName3[20];


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

/* Useful functions when dealing with matrices and vectors 
   and manupulate them */ 

 void VTXsub1(Count_i,Count_k,Buff,Buff1)

/* subtracts two vectors of the same matrix from each other */
 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) 

/* Vector cross product */ 

 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) /*sign */
 int Value;
 {
 double Value1,Value2;
 Value1 = (double)Value ; Value2 = Value1 / fabs(Value1) ;
 return Value2;
 }

 double VTXmag(Buff1,Buff2) /* length of vector */
 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;
 }

/* calculates the "F" terms used in the FEM matrix 
   (see appendices A and B in reference [9] or see reference [1])
   These terms are slightly different when PMLs are used */
   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;}

/* ff1 is used when PML matrix terms are calculated. 
   The PMLs multiply an anisotropic dielectric to the matrix dot-product. 
   So you can see that for each term, the multiplying factor is different */ 

   complex ff1(i,j) int i,j; {
   complex Prod;
   Prod = COMplex_Add2( Real_Mul(CoFactor[1][i-1]*CoFactor[1][j-1],a_PML),
	  		Real_Mul(CoFactor[2][i-1]*CoFactor[2][j-1],b_PML),
			Real_Mul(CoFactor[3][i-1]*CoFactor[3][j-1],c_PML))  ;
              return Prod;}

/* Searches for non-zero entries in the row and column specified */

 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(int RowNum,int ColNum)
     {
      int Count_i,test_i,test_j; {
      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 *********************************/