/*                ***********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++)   {
 HexNum = k*YdiM*XdiM + j*XdiM + i;
 con = pow(-1.0,(double)i)*pow(-1.0,(double)j)*pow(-1.0,(double)k);

 if(con==1) {HexHedraSubDiv1(); EdgeSubDiv1(); GlobalEdgeEndsDiv1();}
 else       {HexHedraSubDiv2(); EdgeSubDiv2(); GlobalEdgeEndsDiv2();}

 for(t=0;t<=4;t++)
 {
  FindTetraHedronVolume();
  ComputeTetraCoFactor();
  S_Matrix();
  T_Matrix();
  TetraSubMatrix();
  FindGlobalMatrix();
 }

                          }
CountHalfBandWidth();
HalfBandMatrix = FLOAT_Matrix(TotInnerEdgeNum+1,HalfBandWidth+1);
FORC_Matrix_Data=FLOAT_Matrix(TotInnerEdgeNum,FRCD_Matrix_MaxColNos);
FORC_Matrix_ColNos=INT_Matrix(TotInnerEdgeNum,FRCD_Matrix_MaxColNos);
RHSVector = FLOAT_Vector(TotInnerEdgeNum);
for(i=0;i<=TotInnerEdgeNum;i++)
for(j=0;j<=HalfBandWidth;j++) {HalfBandMatrix[i][j] =0.0;}
PartitionGlobalMatrix();
ComputeRHSVector();
BandMatrixSolver();
Produce_Output();
return(0);
}

/***************************************************************************/

void ParameterInfo()
{
double FreeSpaceVel,AbsPermeable,AbsPermitt,WaveLength; 

 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;
 printf("OperateFrequency=%g MHz\n",OperateFreq);
 printf("WaveLength=%f m\n",WaveLength );
 printf("FreeSpaceVelocity=%E m/s\n",FreeSpaceVel );
 printf("WaveNumber=%g\n",WaveNumber);

 TotEdgeNum      = 6*XdiM*YdiM*ZdiM + (XdiM+YdiM+ZdiM)
                 + 3*(XdiM*YdiM+YdiM*ZdiM+ZdiM*XdiM);
 TotNodeNum      = (XdiM+1)*(YdiM+1)*(ZdiM+1);
 TotNumHexHedron = XdiM*YdiM*ZdiM;
 TotTetElement = 5 * TotNumHexHedron;
 TotTrngleNum    = 4*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM);
 TotBoundEdgeNum = 6*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM);
 TotBoundNodeNum = 2*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM+1);
 DivisorX = CellDimension; 
 DivisorY = CellDimension; 
 DivisorZ = CellDimension; 

 printf("CellDimension = %f m \n",CellDimension);
 printf("Total Edge Number = %d \n",TotEdgeNum);
 printf("Total Node Number = %d \n",TotNodeNum);
 printf("Total Cube Number = %d \n",TotNumHexHedron);
 printf("Total Tetrahedron Number = %d \n",TotTetElement);
 printf("Total Boundary Edge Number = %d \n",TotBoundEdgeNum);
 printf("Total Boundary Node Number = %d \n",TotBoundNodeNum);

  /* Local edge orientation of tetrahedron edges */
  edge_end[0][0] = 0; edge_end[0][1] = 1;
  edge_end[1][0] = 0; edge_end[1][1] = 2;
  edge_end[2][0] = 0; edge_end[2][1] = 3;
  edge_end[3][0] = 1; edge_end[3][1] = 2;
  edge_end[4][0] = 3; edge_end[4][1] = 1;
  edge_end[5][0] = 2; edge_end[5][1] = 3;

 }

/***************************************************************************
*   FUNCTION        : AssignGlobalCoorDinates()
****************************************************************************/

void AssignGlobalCoorDinates()
 {
  int i,j,k,NodeNum; double Xtemp,Ytemp,Ztemp;
   {
     NodeNum=0; Xtemp=0.0; Ytemp=0.0; Ztemp=0.0;
     for(k=0;k<=ZdiM;k++) {
       for(j=0;j<=YdiM;j++) {
        for(i=0;i<=XdiM;i++) {
          NodeCord[NodeNum][0] = Xtemp;
          NodeCord[NodeNum][1] = Ytemp;
          NodeCord[NodeNum][2] = Ztemp;
          NodeNum++;   Xtemp = Xtemp + DivisorX; }
          Xtemp = 0.0; Ytemp = Ytemp + DivisorY; }
          Ytemp = 0.0; Ztemp = Ztemp + DivisorZ; }
    }
 }

/***************************************************************************
*   FUNCTION        : AssignHexHedronNodeNumbering()
****************************************************************************/

void AssignHexHedronNodeNumbering()
{
 int i;
  {
   start=1;
   for(i=0;i<=TotNumHexHedron-1;i++)
     {
      HexNode[i][0]=start;
      HexNode[i][1]=HexNode[i][0]+1;
      HexNode[i][2]=HexNode[i][0]+XdiM+1;
      HexNode[i][3]=HexNode[i][2]+1;
      HexNode[i][4]=HexNode[i][0]+(XdiM+1)*(YdiM+1);
      HexNode[i][5]=HexNode[i][4]+1;
      HexNode[i][6]=HexNode[i][4]+(XdiM+1);
      HexNode[i][7]=HexNode[i][6]+1;
      if(((i+1)%((YdiM*XdiM))==0)) {start=start+(XdiM+3);}
      else                     {start=start+1;}
      if((start%(XdiM+1))==0)  { start++; }
     }
  }
}

/***************************************************************************
*   FUNCTION        : AssignHexHedronEdgeNumbering()
****************************************************************************/

void AssignHexHedronEdgeNumbering()
  {
  int EdgeNum,xx1,xx2,xx3,kk,dy1,dy2,divz,divz1,i,j,k;
  {
        for(k=0;k<=ZdiM-1;k++) 
        for(j=0;j<=YdiM-1;j++)
        for(i=0;i<=XdiM-1;i++) {  
        EdgeNum = k*YdiM*XdiM + j*XdiM + i;
        for(kk=0;kk<=17;kk++)   {HexEdgeNum[EdgeNum][kk] = 0;}
        divz=(6*XdiM*YdiM+3*XdiM+3*YdiM+1)*k;
        divz1=  6*XdiM*YdiM + 3*XdiM + 3*YdiM + 1 ;
        dy1 = (3*XdiM+1)*j; dy2 = (3*XdiM+2)*j;
        xx1 = 1 + (YdiM+1)*XdiM + (2*XdiM+1)*YdiM;
        xx2 = xx1 + (2*XdiM+1);
        xx3 = xx2 + (XdiM + 1);
        HexEdgeNum[EdgeNum][0] = 1 + i + dy1 + divz;
        HexEdgeNum[EdgeNum][1] = 1 + XdiM + 2*i + dy1 + divz;
        HexEdgeNum[EdgeNum][2] = HexEdgeNum[EdgeNum][1] + 1;  
        HexEdgeNum[EdgeNum][3] = HexEdgeNum[EdgeNum][2] + 1;
        HexEdgeNum[EdgeNum][4] = 1 + 3*XdiM + 1 + i + dy1 + divz;
        HexEdgeNum[EdgeNum][5] = xx1 + 2*i + dy2 + divz;
        HexEdgeNum[EdgeNum][6] = HexEdgeNum[EdgeNum][5] + 1; 
        HexEdgeNum[EdgeNum][7] = HexEdgeNum[EdgeNum][6] + 1;
        HexEdgeNum[EdgeNum][8] = xx2 + i + dy2 + divz;
        HexEdgeNum[EdgeNum][9] = HexEdgeNum[EdgeNum][8] + 1;
        HexEdgeNum[EdgeNum][10]= xx3 + 2*i + dy2 + divz;
        HexEdgeNum[EdgeNum][11]= HexEdgeNum[EdgeNum][10] + 1;
        HexEdgeNum[EdgeNum][12]= HexEdgeNum[EdgeNum][11] + 1; 
        HexEdgeNum[EdgeNum][13]= HexEdgeNum[EdgeNum][0] + divz1;
        HexEdgeNum[EdgeNum][14]= HexEdgeNum[EdgeNum][1] + divz1; 
        HexEdgeNum[EdgeNum][15]= HexEdgeNum[EdgeNum][2] + divz1;
        HexEdgeNum[EdgeNum][16]= HexEdgeNum[EdgeNum][3] + divz1; 
        HexEdgeNum[EdgeNum][17]= HexEdgeNum[EdgeNum][4] + divz1;
                               }
        
 }
}

/***************************************************************************
*   FUNCTION        : HexHedraSubDiv1()
****************************************************************************/

void HexHedraSubDiv1()
{
TetGlobalNodeNum[0][0]=HexNode[HexNum][0];TetGlobalNodeNum[1][0]=HexNode[HexNum][3];
TetGlobalNodeNum[0][1]=HexNode[HexNum][3];TetGlobalNodeNum[1][1]=HexNode[HexNum][5];
TetGlobalNodeNum[0][2]=HexNode[HexNum][2];TetGlobalNodeNum[1][2]=HexNode[HexNum][7];
TetGlobalNodeNum[0][3]=HexNode[HexNum][6];TetGlobalNodeNum[1][3]=HexNode[HexNum][6];
 
TetGlobalNodeNum[2][0]=HexNode[HexNum][0];TetGlobalNodeNum[3][0]=HexNode[HexNum][0];
TetGlobalNodeNum[2][1]=HexNode[HexNum][1];TetGlobalNodeNum[3][1]=HexNode[HexNum][6];
TetGlobalNodeNum[2][2]=HexNode[HexNum][3];TetGlobalNodeNum[3][2]=HexNode[HexNum][4];
TetGlobalNodeNum[2][3]=HexNode[HexNum][5];TetGlobalNodeNum[3][3]=HexNode[HexNum][5];
 
TetGlobalNodeNum[4][0]=HexNode[HexNum][0];TetGlobalNodeNum[4][2]=HexNode[HexNum][3];
TetGlobalNodeNum[4][1]=HexNode[HexNum][5];TetGlobalNodeNum[4][3]=HexNode[HexNum][6];
}

/***************************************************************************
*   FUNCTION        : HexHedraSubDiv2()
****************************************************************************/

void HexHedraSubDiv2()
{
TetGlobalNodeNum[0][0]=HexNode[HexNum][0];TetGlobalNodeNum[1][0]=HexNode[HexNum][1];
TetGlobalNodeNum[0][1]=HexNode[HexNum][1];TetGlobalNodeNum[1][1]=HexNode[HexNum][7];
TetGlobalNodeNum[0][2]=HexNode[HexNum][2];TetGlobalNodeNum[1][2]=HexNode[HexNum][4];
TetGlobalNodeNum[0][3]=HexNode[HexNum][4];TetGlobalNodeNum[1][3]=HexNode[HexNum][5];
 
TetGlobalNodeNum[2][0]=HexNode[HexNum][1];TetGlobalNodeNum[3][0]=HexNode[HexNum][2];
TetGlobalNodeNum[2][1]=HexNode[HexNum][7];TetGlobalNodeNum[3][1]=HexNode[HexNum][6];
TetGlobalNodeNum[2][2]=HexNode[HexNum][3];TetGlobalNodeNum[3][2]=HexNode[HexNum][4];
TetGlobalNodeNum[2][3]=HexNode[HexNum][2];TetGlobalNodeNum[3][3]=HexNode[HexNum][7];

TetGlobalNodeNum[4][0]=HexNode[HexNum][1];TetGlobalNodeNum[4][2]=HexNode[HexNum][4];
TetGlobalNodeNum[4][1]=HexNode[HexNum][2];TetGlobalNodeNum[4][3]=HexNode[HexNum][7];
}

/***************************************************************************
*   FUNCTION        : EdgeSubDiv1()
****************************************************************************/

void EdgeSubDiv1()
{
TetEdgeNum[0][0] = HexEdgeNum[HexNum][2];  
TetEdgeNum[0][1] = HexEdgeNum[HexNum][1]; 
TetEdgeNum[0][2] = HexEdgeNum[HexNum][8];    
TetEdgeNum[0][3] =-HexEdgeNum[HexNum][4];  
TetEdgeNum[0][4] =-HexEdgeNum[HexNum][11];  
TetEdgeNum[0][5] = HexEdgeNum[HexNum][10];   

TetEdgeNum[1][0] = HexEdgeNum[HexNum][9];   
TetEdgeNum[1][1] = HexEdgeNum[HexNum][12];  
TetEdgeNum[1][2] = HexEdgeNum[HexNum][11]; 
TetEdgeNum[1][3] = HexEdgeNum[HexNum][16];  
TetEdgeNum[1][4] =-HexEdgeNum[HexNum][15];  
TetEdgeNum[1][5] =-HexEdgeNum[HexNum][17];   
 
TetEdgeNum[2][0] = HexEdgeNum[HexNum][0];   
TetEdgeNum[2][1] = HexEdgeNum[HexNum][2];   
TetEdgeNum[2][2] = HexEdgeNum[HexNum][6]; 
TetEdgeNum[2][3] = HexEdgeNum[HexNum][3];  
TetEdgeNum[2][4] =-HexEdgeNum[HexNum][7];  
TetEdgeNum[2][5] = HexEdgeNum[HexNum][9]; 
 
TetEdgeNum[3][0] = HexEdgeNum[HexNum][8];   
TetEdgeNum[3][1] = HexEdgeNum[HexNum][5];  
TetEdgeNum[3][2] = HexEdgeNum[HexNum][6]; 
TetEdgeNum[3][3] =-HexEdgeNum[HexNum][14];  
TetEdgeNum[3][4] = HexEdgeNum[HexNum][15];  
TetEdgeNum[3][5] = HexEdgeNum[HexNum][13]; 
 
TetEdgeNum[4][0] = HexEdgeNum[HexNum][6];   
TetEdgeNum[4][1] = HexEdgeNum[HexNum][2];   
TetEdgeNum[4][2] = HexEdgeNum[HexNum][8];  
TetEdgeNum[4][3] =-HexEdgeNum[HexNum][9]; 
TetEdgeNum[4][4] =-HexEdgeNum[HexNum][15];  
TetEdgeNum[4][5] = HexEdgeNum[HexNum][11];  
}

/***************************************************************************
*   FUNCTION        : EdgeSubDiv2()
****************************************************************************/

void EdgeSubDiv2()
{
TetEdgeNum[0][0] = HexEdgeNum[HexNum][0];
TetEdgeNum[0][1] = HexEdgeNum[HexNum][1];
TetEdgeNum[0][2] = HexEdgeNum[HexNum][5];
TetEdgeNum[0][3] = HexEdgeNum[HexNum][2];
TetEdgeNum[0][4] =-HexEdgeNum[HexNum][6];
TetEdgeNum[0][5] = HexEdgeNum[HexNum][8];
 
TetEdgeNum[1][0] = HexEdgeNum[HexNum][9];
TetEdgeNum[1][1] = HexEdgeNum[HexNum][6];
TetEdgeNum[1][2] = HexEdgeNum[HexNum][7];
TetEdgeNum[1][3] =-HexEdgeNum[HexNum][15];
TetEdgeNum[1][4] = HexEdgeNum[HexNum][16];    
TetEdgeNum[1][5] = HexEdgeNum[HexNum][13];   
 
TetEdgeNum[2][0] = HexEdgeNum[HexNum][9];   
TetEdgeNum[2][1] = HexEdgeNum[HexNum][3];     
TetEdgeNum[2][2] = HexEdgeNum[HexNum][2];
TetEdgeNum[2][3] =-HexEdgeNum[HexNum][12];     
TetEdgeNum[2][4] = HexEdgeNum[HexNum][11];     
TetEdgeNum[2][5] =-HexEdgeNum[HexNum][4];  
 
TetEdgeNum[3][0] = HexEdgeNum[HexNum][10];     
TetEdgeNum[3][1] = HexEdgeNum[HexNum][8];
TetEdgeNum[3][2] = HexEdgeNum[HexNum][11];  
TetEdgeNum[3][3] =-HexEdgeNum[HexNum][14];   
TetEdgeNum[3][4] =-HexEdgeNum[HexNum][17];  
TetEdgeNum[3][5] = HexEdgeNum[HexNum][15];  
 
TetEdgeNum[4][0] = HexEdgeNum[HexNum][2];   
TetEdgeNum[4][1] = HexEdgeNum[HexNum][6];   
TetEdgeNum[4][2] = HexEdgeNum[HexNum][9];   
TetEdgeNum[4][3] = HexEdgeNum[HexNum][8];  
TetEdgeNum[4][4] =-HexEdgeNum[HexNum][11];    
TetEdgeNum[4][5] = HexEdgeNum[HexNum][15];  
}

/***************************************************************************
*   FUNCTION        : GlobalEdgeEndsDiv1()
****************************************************************************/

void GlobalEdgeEndsDiv1()
{
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][1]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][0]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][1]=HexNode[HexNum][7];
   }

/***************************************************************************
*   FUNCTION        : GlobalEdgeEndsDiv2()
****************************************************************************/

void GlobalEdgeEndsDiv2()
{
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][1]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][0]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][1]=HexNode[HexNum][7];
}


/***************************************************************************
*   FUNCTION        : Find_Volume_Determinent()
****************************************************************************/

void FindTetraHedronVolume()
{
   int Count_i, Count_j;
   double Buff[4][3],Buff1[3],Buff2[3],Buff3[3],Buff4[3];
   for(Count_i=0;Count_i<=3;Count_i++)
   for(Count_j=0;Count_j<=2;Count_j++)
   Buff[Count_i][Count_j]=NodeCord[TetGlobalNodeNum[t][Count_i]-1][Count_j];
   VTXsub1(1,0,Buff,Buff1);
   VTXsub1(2,0,Buff,Buff2);
   VTXsub1(3,0,Buff,Buff3);
   VTXcross(Buff1,Buff2,Buff4);
   TetVolume=(1.0/6.0)*(Buff3[0]*Buff4[0]+Buff3[1]*Buff4[1]+Buff3[2]*Buff4[2]);
}

/***************************************************************************
*   FUNCTION        : ComputeTetraCoFactor()
****************************************************************************/

void ComputeTetraCoFactor()
   {
   int i,j,RowNum,ColNum,Row,Col,Count=0;
   double CoF[4][4],Buff[4][4],CoFs[4][4]; 
     {
    for(i=0;i<=3;i++)  for(j=0;j<=3;j++) {
    if(j==3) Buff[i][j]= 1.0;
    else     Buff[i][j]= NodeCord[TetGlobalNodeNum[t][i]-1][j];
                                         }
     for(RowNum=0;RowNum<=3;RowNum++)
     for(ColNum=0;ColNum<=3;ColNum++) {
     Row=0;Col=0;Count=0;
        for(i=0;i<=3;i++)
        for(j=0;j<=3;j++) {
        if((i!=RowNum) && (j!=ColNum)) {
        Count++;
        if(Count<4)              Row=0;
        if((Count>3)&&(Count<7)) Row=1;
        else if(Count>6)         Row=2;
        CoF[Row][Col]=Buff[i][j];
        Col++; Col=Col%3;              }
                          }
     Row=0;Col=0;Count=0;
     CoFs[RowNum][ColNum]
     =CoF[0][0]*((CoF[1][1]*CoF[2][2])-(CoF[1][2]*CoF[2][1]))
     -CoF[0][1]*((CoF[1][0]*CoF[2][2])-(CoF[2][0]*CoF[1][2]))
     +CoF[0][2]*((CoF[1][0]*CoF[2][1])-(CoF[2][0]*CoF[1][1]));
     if((RowNum+ColNum)%2!=0)  CoFs[RowNum][ColNum]*=-1.0;
                                      }
     for(RowNum=0;RowNum<=3;RowNum++) {
     if     (RowNum==0) Col = 3;
     else if(RowNum==1) Col = 0;
     else if(RowNum==2) Col = 1;
     else if(RowNum==3) Col = 2;
     for(ColNum=0;ColNum<=3;ColNum++) {
     CoFactor[RowNum][ColNum] = -CoFs[ColNum][Col];
                                      }}
      }
    }

/***************************************************************************
*   FUNCTION        : S_Matrix()
****************************************************************************/

void S_Matrix()
  {
  int i,j;  double Length_i,Length_j,Mult1,Volume;
   {
   Volume = TetVolume; Mult1 = (1296.0*Volume*Volume*Volume);
   for(i=0;i<=5;i++) for(j=0;j<=5;j++)  {
  Length_i = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[i][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[i][1]]-1]);
  Length_j = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[j][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[j][1]]-1]);
   S_mat[i][j] = ( CoFactor[2][edge_end[i][0]] * CoFactor[3][edge_end[i][1]]
                 - CoFactor[3][edge_end[i][0]] * CoFactor[2][edge_end[i][1]] )
               * ( CoFactor[2][edge_end[j][0]] * CoFactor[3][edge_end[j][1]]
                 - CoFactor[3][edge_end[j][0]] * CoFactor[2][edge_end[j][1]] )
               + ( CoFactor[3][edge_end[i][0]] * CoFactor[1][edge_end[i][1]]
                 - CoFactor[1][edge_end[i][0]] * CoFactor[3][edge_end[i][1]] )
               * ( CoFactor[3][edge_end[j][0]] * CoFactor[1][edge_end[j][1]]
                 - CoFactor[1][edge_end[j][0]] * CoFactor[3][edge_end[j][1]] )
               + ( CoFactor[1][edge_end[i][0]] * CoFactor[2][edge_end[i][1]]
                 - CoFactor[2][edge_end[i][0]] * CoFactor[1][edge_end[i][1]] )
               * ( CoFactor[1][edge_end[j][0]] * CoFactor[2][edge_end[j][1]]
                 - CoFactor[2][edge_end[j][0]] * CoFactor[1][edge_end[j][1]] );
   S_mat[i][j] = (Length_i * Length_j * S_mat[i][j])/ Mult1 ;
                                         }
   }
   }

/***************************************************************************
*   FUNCTION        : T_Matrix()
****************************************************************************/

void T_Matrix()
  {
  int i,j;
  double Mult1,Length_i,Length_j;
 
   Mult1 = 1.0/(720.0 * TetVolume);
   T_mat[0][0] =               ff(1,1) + ff(2,2) - ff(1,2);
   T_mat[0][1] = 2 * ff(2,3) - ff(2,1) - ff(1,3) + ff(1,1);
   T_mat[0][2] = 2 * ff(2,4) - ff(2,1) - ff(1,4) + ff(1,1);
   T_mat[0][3] =-2 * ff(1,3) - ff(2,2) + ff(2,3) + ff(1,2);
   T_mat[0][4] = 2 * ff(1,4) - ff(2,4) - ff(1,2) + ff(2,2);
   T_mat[0][5] =     ff(2,4) - ff(2,3) - ff(1,4) + ff(1,3);
   T_mat[1][1] =               ff(1,1) + ff(3,3) - ff(1,3);
   T_mat[1][2] = 2 * ff(3,4) - ff(1,3) - ff(1,4) + ff(1,1);
   T_mat[1][3] = 2 * ff(1,2) - ff(2,3) - ff(1,3) + ff(3,3);
   T_mat[1][4] =     ff(2,3) - ff(3,4) - ff(1,2) + ff(1,4);
   T_mat[1][5] =-2 * ff(1,4) - ff(3,3) + ff(1,3) + ff(3,4);
   T_mat[2][2] =               ff(1,1) + ff(4,4) - ff(1,4);
   T_mat[2][3] =     ff(3,4) - ff(2,4) - ff(1,3) + ff(1,2);
   T_mat[2][4] =-2 * ff(1,2) - ff(4,4) + ff(1,4) + ff(2,4);
   T_mat[2][5] = 2 * ff(1,3) - ff(3,4) - ff(1,4) + ff(4,4);
   T_mat[3][3] =               ff(3,3) + ff(2,2) - ff(2,3);
   T_mat[3][4] =-2 * ff(3,4) - ff(2,2) + ff(2,3) + ff(2,4);
   T_mat[3][5] =-2 * ff(2,4) - ff(3,3) + ff(2,3) + ff(3,4);
   T_mat[4][4] =               ff(2,2) + ff(4,4) - ff(2,4);
   T_mat[4][5] =-2 * ff(2,3) - ff(4,4) + ff(2,4) + ff(3,4);
   T_mat[5][5] =               ff(3,3) + ff(4,4) - ff(3,4);
   for(i=0;i<=5;i++) for(j=0;j<=i-1;j++)
   {T_mat[i][j] = T_mat[j][i];}
 
   for(i=0;i<=5;i++) for(j=0;j<=5;j++) {
  Length_i = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[i][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[i][1]]-1]);
  Length_j = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[j][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[j][1]]-1]);
   if(i==j) T_mat[i][j] = 2 * Length_i * Length_j * Mult1 * T_mat[i][j];
   else     T_mat[i][j] =     Length_i * Length_j * Mult1 * T_mat[i][j];
   T_mat[i][j] = RELPerm[HexNum][0] * T_mat[i][j] ;
   T_mat[i][j] = WaveNumber*WaveNumber * T_mat[i][j] ;
                                       }
    }

/***************************************************************************
*   FUNCTION        : TetraSubMatrix()
****************************************************************************/

void TetraSubMatrix()
  {
  int i,j;
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j]=0.0;}
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j] =4.0*S_mat[i][j]-T_mat[i][j];}
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j] = A_mat[t][i][j] * Sign(TetEdgeNum[t][j]);}
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j] = A_mat[t][i][j] * Sign(TetEdgeNum[t][i]);}
  }

/***************************************************************************
*   FUNCTION        : FindGlobalMatrix()
****************************************************************************/

void FindGlobalMatrix()
{
 int m,n,k,l,rval;
  {
  for(m=0;m<=5;m++) for(n=0;n<=5;n++)   {
   k=absol(TetEdgeNum[t][m])-1; l=absol(TetEdgeNum[t][n])-1;
   if(fabs(A_mat[t][m][n])>1.0E-10) {
   rval=Srch_NZ_Element_Column(k,l);
   if(rval>=0)                {
   GBL_Matrix_Data[k][rval]+=A_mat[t][m][n];
                              }
   else                       {
   GBL_Matrix_ColNos[k][GBL_Matrix_Index[k]]=l;
   GBL_Matrix_Data[k][GBL_Matrix_Index[k]]=A_mat[t][m][n];
   GBL_Matrix_Index[k]++;
                              }
                                    }

                                        }
  }
}

/***************************************************************************
*   FUNCTION        : CountHalfBandWidth()
****************************************************************************/

void CountHalfBandWidth()
     {
      short int i,j,k,l,ColDiff,Val;
        {
        for(i=0;i<TotInnerEdgeNum;i++) {k=InnerEdgeStat[i];
        for(j=i;j<TotInnerEdgeNum;j++) {l=InnerEdgeStat[j];
        Val=Srch_NZ_Element_Column(k,l);
        if(Val>=0) { ColDiff=j-i;
        if(HalfBandWidth<=(ColDiff+1))
        {HalfBandWidth=(ColDiff+1);}
                   }
                                              }
                                              }
      }
    }

/***************************************************************************
*   FUNCTION      : PartitionGlobalMatrix()
****************************************************************************/

void PartitionGlobalMatrix()
 {
   int i,j,k,l,ColDiff,val,g=0;
   {
    for(i=0;i<TotInnerEdgeNum;i++) {k=InnerEdgeStat[i];
    for(j=i;j<TotInnerEdgeNum;j++) {l=InnerEdgeStat[j];
    val=Srch_NZ_Element_Column(k,l);
    if((val>=0)) { ColDiff=j-i;
    HalfBandMatrix[i+1][ColDiff+1]=GBL_Matrix_Data[k][val];
                 }
    else  {
    ColDiff=j-i; HalfBandMatrix[i+1][ColDiff+1]=0.0;}
                                             }
                                             }
    for(i=0;i<TotInnerEdgeNum;i++) {k=InnerEdgeStat[i];
    for(j=0;j<TotForcdEdgeNum;j++) {l=ForcdEdgeStat[j];
    val=Srch_NZ_Element_Column(k,l);
    if(val>=0)  {
    FORC_Matrix_ColNos[i][g]=j;
    FORC_Matrix_Data[i][g]=(GBL_Matrix_Data[k][val]);
    g++;        }
                                             }
    g=0;                                     } g=0;
   }
  }

/***************************************************************************
*   FUNCTION        : ComputeRHSVector()
****************************************************************************/

void ComputeRHSVector()
    {
     int Count_i,Count_j,Srch_Ptr=0;
     double Temp_Buff,Multpl;
     {
      for(Count_i=0;Count_i<TotInnerEdgeNum;Count_i++) {
      for(Count_j=0;Count_j<TotForcdEdgeNum;Count_j++)   {
      Srch_Ptr=Locate_NZ_Element_Column(Count_i,Count_j);
      if(Srch_Ptr>=0) { Temp_Buff=FORC_Matrix_Data[Count_i][Srch_Ptr];}
      else            { Temp_Buff=0.0;}
      if(Temp_Buff!=0.0) {
      Multpl=Temp_Buff*ForcdValue[Count_j][0]; RHSVector[Count_i]+=Multpl;
                         }
                                                             }
      RHSVector[Count_i]=RHSVector[Count_i]*-1.0;          }
      }
    }
 
/***************************************************************************
*   FUNCTION        : BandMatrixSolver()
****************************************************************************/

void BandMatrixSolver()
  {
  int Count_i,Count_j,Count_k,Count_l,Count_m,Count_n;
  int Row_Number,Column_Number;
  double Pivot_Emement,Multpl_1,Multpl_2;
  {
 
 Row_Number=TotInnerEdgeNum;
 Column_Number=HalfBandWidth;
 for(Count_i=1;Count_i<=Row_Number;Count_i++)  {
 for(Count_j=2;Count_j<=Column_Number;Count_j++) {
 if(HalfBandMatrix[Count_i][Count_j]!=0.0) {
 Pivot_Emement=HalfBandMatrix[Count_i][Count_j]/HalfBandMatrix[Count_i][1];
 Count_n=Count_i+Count_j-1;
 Count_l=0;
 for(Count_k=Count_j;Count_k<=Column_Number;Count_k++) {
 Count_l=Count_l+1;
 Multpl_1=Pivot_Emement*HalfBandMatrix[Count_i][Count_k];
 HalfBandMatrix[Count_n][Count_l] -= Multpl_1;
                                                       }
 Multpl_2=Pivot_Emement*RHSVector[Count_i-1];
 RHSVector[Count_n-1] -= Multpl_2;         }
                                                 }
                                               }
 RHSVector[Row_Number-1] /= HalfBandMatrix[Row_Number][1];
 for(Count_m=2;Count_m<=Row_Number;Count_m++)  {
 Count_i=Row_Number+1-Count_m;
 for(Count_j=2;Count_j<=Column_Number;Count_j++)  {
 if(HalfBandMatrix[Count_i][Count_j]!=0.0)  {
 Count_k=Count_i+Count_j-1;
 Multpl_1=HalfBandMatrix[Count_i][Count_j]*RHSVector[Count_k-1];
 RHSVector[Count_i-1] -= Multpl_1;          }
                                                  }
 RHSVector[Count_i-1] /= HalfBandMatrix[Count_i][1];
                                               }
 }
 }

/***************************************************************************
*   FUNCTION        : Read_Input_Pass_1()
****************************************************************************/

void Read_Input_Pass_1()
{
 int    x1, y1, z1, x2, y2, z2, OutF_Count=0, output_flag=0;
 char   type[20], att1[20], att2[18], att3[18], att4[18];
 float  tmp1, tmp2, cell_dim = 1;
 double freq = 0;
 char            buffer[80];

 Min_X = 10000; Min_Y=10000; Min_Z=10000;
 Max_X = 0;     Max_Y=0;     Max_Z=0;

   while (fgets(buffer, 80, InF))
   {
        if (buffer[0] == '#') continue;

        if(!strncmp(buffer, "default_output",14))
                   {
                   sscanf(buffer, "%s%s", type, Out_FileName0);
                   OutF_0=fopen(Out_FileName0, "w");
                   if (OutF_0 == NULL) {fprintf(stderr, " Error: Can't create default output file %s\n", Out_FileName0); exit(1);}
                   output_flag=1;
                   continue;
                   }
        if(!strncmp(buffer, "celldim", 7)) { 
                                              sscanf(buffer, "%s%g%s", type,
                                              &tmp1, att1);
                                              if (!strcmp(att1, "cm")) tmp2=0.01;
                                              else if(!strcmp(att1, "m"))
                                                      tmp2=1;
                                              else if(!strcmp(att1, "mm"))
                                                      tmp2=0.001;
                                              else if(!strcmp(att1, "in"))
                                                      tmp2=0.0254;
                                              else {fprintf(stderr,"ERROR: unrecognized cell dimension.\n");
                                                    exit(1);}
                                              cell_dim=tmp1*tmp2;
                                              continue;
                                            }

        if(sscanf(buffer, "%s%d%d%d%d%d%d%s%s%s%s", type, &x1, &y1, &z1, &x2, &y2, &z2, att1, att2, att3, att4)>2)
           {
              if(!strcmp(type,"boundary")) {fprintf(stderr, "\n\nERROR: boundary keyword is not allowed by EMAP3 \n"); exit(1);}
              if(!strcmp(type, "dielectric"))
              {
               if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
               if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
               if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
               if((x1==x2)||(y1==y2)||(z1==z2)) {fprintf(stderr, "ERROR:  diel thickness is not finite./n");
               exit(1); } continue;
              }
              if(!strcmp(type, "box"))
              {
                   if((x1==x2)||(y1==y2)||(z1==z2)){ fprintf(stderr, "ERROR: box dimension is invalid.");
                   exit(1);}
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   continue;
              }
              if(!strcmp(type,"conductor"))
              {
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   continue;
              }
              if(!strcmp(type, "aperture"))
              {
                   if((x1!=x2) & (y1!=y2) & (z1!=z2)){fprintf(stderr,"ERROR: aperture dimension is invalid.\n");
                   exit(1); }
                   if(((x1==x2)&&(y1==y2))||((x1==x2)&&(z1==z2))||((z1==z2)&&(y1==y2)))
                   {fprintf(stderr, "ERROR: aperture dimension is invalid.\n"); exit(1); }
                   continue;
              }
              if(!strcmp(type, "seam")) continue;
              if(!strcmp(type,"esource"))
              {
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   freq = atof(att1);
                   if((x1==x2)&&(y1==y2)&&(z1==z2)) {
                   fprintf(stderr, "ERROR: esource cannot be defined at a point.\n"); exit(1); }
                   if ((att2[0] != 'x')  && (att2[0] != 'y') && (att2[0] != 'z')){
                   fprintf(stderr, "ERROR: unrecognized polarization of esource statement  %c\n",att2[0]);
                   exit(1); }
                   continue;
              }
              if(!strcmp(type, "msource"))
              {
                   if (x1<Min_X) Min_X=x1; if (x2<Min_X) Min_X=x2; if (x1>Max_X) Max_X=x1; if (x2>Max_X) Max_X=x2;
                   if (y1<Min_Y) Min_Y=y1; if (y2<Min_Y) Min_Y=y2; if (y1>Max_Y) Max_Y=y1; if (y2>Max_Y) Max_Y=y2;
                   if (z1<Min_Z) Min_Z=z1; if (z2<Min_Z) Min_Z=z2; if (z1>Max_Z) Max_Z=z1; if (z2>Max_Z) Max_Z=z2;
                   freq = atof(att1);
                   continue;
              }
              if(!strcmp(type, "efield_output"))
              {    
                   switch (OutF_Count)
                        {
                        case 0:  {OutF_Count=1; strcpy(Out_FileName1, att1); OutF_1=fopen(Out_FileName1,"w"); break;}
                        case 1:  {if(strcmp(att1, Out_FileName1)){OutF_Count=2; strcpy(Out_FileName2, att1); OutF_2=fopen(Out_FileName2,"w");}
                                  break;
                                 }
                        case 2:  {if(strcmp(att1, Out_FileName1) && strcmp(att1, Out_FileName2))
                                      {OutF_Count=3; strcpy(Out_FileName3, att1); OutF_3=fopen(Out_FileName3,"w");}
                                  break;
                                 }
                        case 3:  {if((strcmp(att1, Out_FileName1) && strcmp(att1, Out_FileName2)) && strcmp(att1, Out_FileName3))
                                      {fprintf(stderr, "ERROR: too many different files specified by efield_output commands \n"); exit(1);}
                                  break;
                                 }
                        }
                   output_flag=1;
                   continue;
              }
              fprintf(stderr,"Warning: Unrecognized Keyword %s\n", type);
           }                 /*end of if*/
           OperateFreq = freq;
    }                    /* end of while */
    if(freq == 0) fprintf(stderr,"\nWARNING: The input file does not specify a source, code will not produce output!\n\n");
    if(output_flag == 0) fprintf(stderr,"\nWARNING: The input file does not contain an output statement, code will not produce output!\n\n");
    rewind(InF);
    XdiM = Max_X - Min_X; YdiM = Max_Y - Min_Y; ZdiM = Max_Z - Min_Z;
    OperateFreq = freq;
    CellDimension = cell_dim;
}

/***************************************************************************
*   FUNCTION        : Read_Input_Pass_2()
****************************************************************************/

void Read_Input_Pass_2()
{
 void swap();
 void edge();
 void deter_edge();
 void CountEdgeType();
 int    i, j, k, x1, y1, z1, x2, y2, z2, nn, tmp;
 char   type[20], att1[20], att2[18], att3[18], att4[18];
 char   buffer[80];

        /* Dynamically allocate memory for large variables */
        Epsilon= (double *) calloc( TotEdgeNum*sizeof(double), sizeof(double));
        Sigma= (double *) calloc( TotEdgeNum*sizeof(double), sizeof(double));
        EdgeStatus   = (double *) calloc( TotEdgeNum*sizeof(double), sizeof(double));
        ForceStat= (int *) calloc( TotEdgeNum*sizeof(int), sizeof(int));
        xn   = (int **) calloc(TotEdgeNum*sizeof(int),sizeof(int));
        for(i=0;i<TotEdgeNum;i++)
        {xn[i] = (int *) calloc(18*sizeof(int), sizeof(int));}

        /*Initialize dynamically allocated variables */
        for (i=0;i<=TotEdgeNum;i++) { ForceStat[i]=0.0;
                                      EdgeStatus[i]=0.0;
                                    }
        for(i=0;i<XdiM;i++) for(j=0;j<YdiM;j++) for(k=0;k<ZdiM;k++) {
        nn=k*YdiM*XdiM+j*XdiM+i; Epsilon[nn]=1.0; Sigma[nn]=0.0;
        edge(nn,i,j,k);  
        deter_edge(nn,i,j,k);                                       }
        
        /*Begin reading the input file for the second time */
        while (fgets(buffer, 80, InF))
        {
            if (buffer[0] == '#') continue;
            if(sscanf(buffer, "%s%d%d%d%d%d%d%s%s%s%s", 
            type, &x1, &y1, &z1, &x2, &y2, &z2, att1, att2, att3, att4)>2)
             {
             if(!strcmp(type, "boundary"))
             {
                   continue;
             }
             if(!strcmp(type, "efield_output"))
             {     
                   if((x1<Min_X || x1>Max_X)||(x2<Min_X || x2>Max_X)) {fprintf(stderr,"ERROR: efield_output parameter out of range.\n"); exit(1);}
                   if((y1<Min_Y || y1>Max_Y)||(y2<Min_Y || y2>Max_Y)) {fprintf(stderr,"ERROR: efield_output parameter out of range.\n"); exit(1);}
                   if((z1<Min_Z || z1>Max_Z)||(z2<Min_Z || z2>Max_Z)) {fprintf(stderr,"ERROR: efield_output parameter out of range.\n"); exit(1);}
                   continue;
             }
             if(!strcmp(type, "dielectric"))
             {
                   if (x2<x1) swap(&x1, &x2);
                   if (y2<y1) swap(&y1, &y2);
                   if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;
                   for(i=x1;i<x2;i++){
                      for(j=y1;j<y2;j++){
                         for(k=z1;k<z2;k++){
                         nn=k*XdiM*YdiM+j*XdiM+i;
                         Epsilon[nn]=atof(att1);
                         Sigma[nn]=atof(att2);
                                            }
                                         }
                                      } 
                   continue;
             }
             if(!strcmp(type, "box"))
             {
              if (x2<x1) swap(&x1, &x2);
              if (y2<y1) swap(&y1, &y2);
              if (z2<z1) swap(&z1, &z2);
              x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
              x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;
              /* left plate */
                           for(k=z1;k<z2;k++){
                             for(j=y1;j<y2;j++){
                                                nn=x2*y2*k+j*x2;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][8]]=1;
                                                EdgeStatus[xn[nn][8]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                               }
                                              }
              /* right plate */
                           for(k=z1;k<z2;k++){
                             for(j=y1;j<y2;j++){
                                                nn=x2*y2*k+(x2-1)+j*x2;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][9]]=1;
                                                EdgeStatus[xn[nn][9]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                               }
                                              }
              /* bottom plate */
                           for(k=z1;k<z2;k++){
                             for(i=x1;i<x2;i++){
                                                nn=y2*x2*k+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][6]]=1;
                                                EdgeStatus[xn[nn][6]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                               }
                                              }
              /* top plate */
                           for(k=z1;k<z2;k++){
                             for(i=x1;i<x2;i++){
                                                nn=y2*x2*k+x2*(y2-1)+i;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][11]]=1;
                                                EdgeStatus[xn[nn][11]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                               }
                                              }

              /* front plate */
                           for(j=y1;j<y2;j++){
                             for(i=x1;i<x2;i++){
                                                nn=j*x2+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][2]]=1;
                                                EdgeStatus[xn[nn][2]]=0.0;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                               }
                                              }
              /* back plate */
                           for(j=y1;j<y2;j++){
                             for(i=x1;i<x2;i++){
                                                nn=(z2-1)*y2*x2+j*x2+i;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                                ForceStat[xn[nn][15]]=1;
                                                EdgeStatus[xn[nn][15]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                               }
                                              }
             }
             if(!strcmp(type, "conductor"))
             {
                   /* fprintf(stderr, "conductor \n"); */
                   if (x2<x1) swap(&x1, &x2);
                   if (y2<y1) swap(&y1, &y2);
                   if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;

              /* three-dimensional conductor */
                   for(i=x1;i<x2;i++){
                     for(j=y1;j<y2;j++){
                        for(k=z1;k<z2;k++){
                                           nn=k*XdiM*YdiM+j*XdiM+i;
                                           Epsilon[nn]=0;
                                           Sigma[nn]=1.0E+09;
                                           for(tmp=1;tmp<19;tmp++){
                                                                  ForceStat[xn[nn][tmp]]=1;
                                                                  EdgeStatus[xn[nn][tmp]]=0.0;
                                                                  }
                                          }
                                       }
                                     }
               /* two-dimensional conductor */
                  if(x1==x2){     /*  y-z plane */
                           for(j=y1;j<y2;j++){
                             for(k=z1;k<z2;k++){
                                                if(x1!=XdiM)   {
                                                nn=k*XdiM*YdiM+j*XdiM+x1;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][8]]=1;
                                                EdgeStatus[xn[nn][8]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                                            }
                                                if(x1!=0)   {
                                                nn=k*XdiM*YdiM+j*XdiM+x1-1;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][9]]=1;
                                                EdgeStatus[xn[nn][9]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                                           }
                                               }
                                              }
                           }
                 if(y1==y2){     /*  x-z plane */
                           for(i=x1;i<x2;i++){
                             for(k=z1;k<z2;k++){
                                                if(y1!=YdiM)   {
                                                nn=k*XdiM*YdiM+y1*XdiM+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][6]]=1;
                                                EdgeStatus[xn[nn][6]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                                            }
                                                if(y1!=0)   {
                                                nn=k*XdiM*YdiM+(y1-1)*XdiM+i;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][11]]=1;
                                                EdgeStatus[xn[nn][11]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                                           }
                                               }
                                              }
                           }
                 if(z1==z2){     /*   x-y plane */
                           for(i=x1;i<x2;i++){
                             for(j=y1;j<y2;j++){
                                                if(z1!=ZdiM)   {
                                                nn=z1*XdiM*YdiM+j*XdiM+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][2]]=1;
                                                EdgeStatus[xn[nn][2]]=0.0;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                                            }
                                                if(z1!=0)   {
                                                nn=(z1-1)*XdiM*YdiM+j*XdiM+i;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                                ForceStat[xn[nn][15]]=1;
                                                EdgeStatus[xn[nn][15]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                                           }
                                               }
                                              }
                           }
 
               /* one-dimensional conductor */
 
                  if ((x1==x2)&&(y1==y2)){
                                          for(k=z1;k<z2;k++){
                                                            if((x1!=XdiM) & (y1!=YdiM)) {
                                                             nn=k*XdiM*YdiM+y1*XdiM+x1;
                                                             ForceStat[xn[nn][5]]=1;
                                                             EdgeStatus[xn[nn][5]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (y1!=YdiM)) {
                                                             nn=k*XdiM*YdiM+y1*XdiM+x1-1;
                                                             ForceStat[xn[nn][7]]=1;
                                                             EdgeStatus[xn[nn][7]]=0.0;
                                                                                  }
                                                            if((x1!=XdiM) & (y1!=0)) {
                                                             nn=k*XdiM*YdiM+(y1-1)*XdiM+x1;
                                                             ForceStat[xn[nn][10]]=1;
                                                             EdgeStatus[xn[nn][10]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (y1!=0)) {
                                                             nn=k*XdiM*YdiM+(y1-1)*XdiM+x1-1;
                                                             ForceStat[xn[nn][12]]=1;
                                                             EdgeStatus[xn[nn][12]]=0.0;
                                                                                  }
                                                            }
                                         }
                  if ((x1==x2)&&(z1==z2)){
                                          for(j=y1;j<y2;j++){
                                                            if((x1!=XdiM) & (z1!=ZdiM)) {
                                                             nn=z1*XdiM*YdiM+j*XdiM+x1;
                                                             ForceStat[xn[nn][1]]=1;
                                                             EdgeStatus[xn[nn][1]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (z1!=ZdiM)) {
                                                             nn=z1*XdiM*YdiM+j*XdiM+x1-1;
                                                             ForceStat[xn[nn][3]]=1;
                                                             EdgeStatus[xn[nn][3]]=0.0;
                                                                                  }
                                                            if((x1!=XdiM) & (z1!=0)) {
                                                             nn=(z1-1)*XdiM*YdiM+j*XdiM+x1;
                                                             ForceStat[xn[nn][14]]=1;
                                                             EdgeStatus[xn[nn][14]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (z1!=0)) {
                                                             nn=(z1-1)*XdiM*YdiM+j*XdiM+x1-1;
                                                             ForceStat[xn[nn][16]]=1;
                                                             EdgeStatus[xn[nn][16]]=0.0;
                                                                                  }
                                                            }
                                         }
                  if ((y1==y2)&&(z1==z2)){
                                          for(i=x1;i<x2;i++){
                                                            if((y1!=YdiM) & (z1!=ZdiM)) {
                                                             nn=z1*XdiM*YdiM+y1*XdiM+i;
                                                             ForceStat[xn[nn][0]]=1;
                                                             EdgeStatus[xn[nn][0]]=0.0;
                                                                                  }
                                                            if((y1!=0) & (z1!=ZdiM)) {
                                                             nn=z1*XdiM*YdiM+(y1-1)*XdiM+i;
                                                             ForceStat[xn[nn][4]]=1;
                                                             EdgeStatus[xn[nn][4]]=0.0;
                                                                                  }
                                                            if((y1!=YdiM) & (z1!=0)) {
                                                             nn=(z1-1)*XdiM*YdiM+y1*XdiM+i;
                                                             ForceStat[xn[nn][13]]=1;
                                                             EdgeStatus[xn[nn][13]]=0.0;
                                                                                  }
                                                            if((y1!=0) & (z1!=0)) {
                                                             nn=(z1-1)*XdiM*YdiM+(y1-1)*XdiM+i;
                                                             ForceStat[xn[nn][17]]=1;
                                                             EdgeStatus[xn[nn][17]]=0.0;
                                                                                  }
                                                            }
                                         }
 
                   continue;
             }

             if(!strcmp(type, "aperture"))
             {
                   if (x2<x1) swap(&x1, &x2);
                   if (y2<y1) swap(&y1, &y2);
                   if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;
 
                  if(x1==x2){                       /*  y-z plane */
                           for(j=y1;j<y2;j++){
                             for(k=z1;k<z2;k++){
                                                if(x1!=XdiM)   {
                                                nn=k*XdiM*YdiM+j*XdiM+x1;
                                                      if(k!=z1)     {
                                                      ForceStat[xn[nn][1]]=0;
                                                      EdgeStatus[xn[nn][1]]=-1;
                                                                    }
                                                      if(j!=y1)     {
                                                      ForceStat[xn[nn][5]]=0;
                                                      EdgeStatus[xn[nn][5]]=-1;
                                                                    }
                                                      ForceStat[xn[nn][8]]=0;
                                                      EdgeStatus[xn[nn][8]]=-1;
                                                      if(j!=(y2-1)) {
                                                      ForceStat[xn[nn][10]]=0;
                                                      EdgeStatus[xn[nn][10]]=-1;
                                                                    }
                                                      if(k!=(z2-1)) {
                                                      ForceStat[xn[nn][14]]=0;
                                                      EdgeStatus[xn[nn][14]]=-1;
                                                                    }
                                                             }
                                                if(x1!=0)   {
                                                nn=k*XdiM*YdiM+j*XdiM+x1-1;
                                                      if(k!=z1)     {
                                                      ForceStat[xn[nn][3]]=0;
                                                      EdgeStatus[xn[nn][3]]=-1;
                                                                    }
                                                      if(j!=y1)     {
                                                      ForceStat[xn[nn][7]]=0;
                                                      EdgeStatus[xn[nn][7]]=-1;
                                                                    }
                                                      ForceStat[xn[nn][9]]=0;
                                                      EdgeStatus[xn[nn][9]]=-1;
                                                      if(j!=(y2-1)) {
                                                      ForceStat[xn[nn][12]]=0;
                                                      EdgeStatus[xn[nn][12]]=-1;
                                                                    }
                                                      if(k!=(z2-1)) {
                                                      ForceStat[xn[nn][16]]=0;
                                                      EdgeStatus[xn[nn][16]]=-1;
                                                                    }
                                                           }
                                               }
                                              }
                           }
                 if(y1==y2){                             /*  x-z plane */
                           for(i=x1;i<x2;i++){
                             for(k=z1;k<z2;k++){
                                                if(y1!=YdiM)   {
                                                nn=k*XdiM*YdiM+y1*XdiM+i;
                                                      if(k!=z1)     {
                                                      ForceStat[xn[nn][0]]=0;
                                                      EdgeStatus[xn[nn][0]]=-1;
                                                                    }
                                                      if(i!=x1)     {
                                                      ForceStat[xn[nn][5]]=0;
                                                      EdgeStatus[xn[nn][5]]=-1;
                                                                    }
                                                      ForceStat[xn[nn][6]]=0;
                                                      EdgeStatus[xn[nn][6]]=-1;
                                                      if(i!=(x2-1)) {
                                                      ForceStat[xn[nn][7]]=0;
                                                      EdgeStatus[xn[nn][7]]=-1;
                                                                    }
                                                      if(k!=(z2-1)) {
                                                      ForceStat[xn[nn][13]]=0;
                                                      EdgeStatus[xn[nn][13]]=-1;
                                                                    }
                                                            }
                                                if(y1!=0)   {
                                                nn=k*XdiM*YdiM+(y1-1)*XdiM+i;
                                                      if(k!=z1)     {
                                                      ForceStat[xn[nn][4]]=0;
                                                      EdgeStatus[xn[nn][4]]=-1;
                                                                    }
                                                      if(i!=x1)     {
                                                      ForceStat[xn[nn][10]]=0;
                                                      EdgeStatus[xn[nn][10]]=-1;
                                                                    }
                                                      ForceStat[xn[nn][11]]=0;
                                                      EdgeStatus[xn[nn][11]]=-1;
                                                      if(i!=(x2-1)) {
                                                      ForceStat[xn[nn][12]]=0;
                                                      EdgeStatus[xn[nn][12]]=-1;
                                                                    }
                                                      if(k!=(z2-1)) {
                                                      ForceStat[xn[nn][17]]=0;
                                                      EdgeStatus[xn[nn][17]]=-1;
                                                                    }
                                                           }
                                               }
                                              }
                           }
                 if(z1==z2){                           /*   x-y plane */
                           for(i=x1;i<x2;i++){
                             for(j=y1;j<y2;j++){
                                                if(z1!=ZdiM)   {
                                                nn=z1*XdiM*YdiM+j*XdiM+i;
                                                     if(j!=y1)     {
                                                      ForceStat[xn[nn][0]]=0;
                                                      EdgeStatus[xn[nn][0]]=-1;
                                                                    }
                                                      if(i!=x1)     {
                                                      ForceStat[xn[nn][1]]=0;
                                                      EdgeStatus[xn[nn][1]]=-1;
                                                                    }
                                                      ForceStat[xn[nn][2]]=0;
                                                      EdgeStatus[xn[nn][2]]=-1;
                                                      if(i!=(x2-1)) {
                                                      ForceStat[xn[nn][3]]=0;
                                                      EdgeStatus[xn[nn][3]]=-1;
                                                                    }
                                                      if(j!=(y2-1)) {
                                                      ForceStat[xn[nn][4]]=0;
                                                      EdgeStatus[xn[nn][4]]=-1;
                                                                    }
                                                            }
                                                if(z1!=0)   {
                                                nn=(z1-1)*XdiM*YdiM+j*XdiM+i;
                                                     if(j!=y1)     {
                                                      ForceStat[xn[nn][13]]=0;
                                                      EdgeStatus[xn[nn][13]]=-1;
                                                                    }
                                                      if(i!=x1)     {
                                                      ForceStat[xn[nn][14]]=0;
                                                      EdgeStatus[xn[nn][14]]=-1;
                                                                    }
                                                      ForceStat[xn[nn][15]]=0;
                                                      EdgeStatus[xn[nn][15]]=-1;
                                                      if(i!=(x2-1)) {
                                                      ForceStat[xn[nn][16]]=0;
                                                      EdgeStatus[xn[nn][16]]=-1;
                                                                    }
                                                      if(j!=(y2-1)) {
                                                      ForceStat[xn[nn][17]]=0;
                                                      EdgeStatus[xn[nn][17]]=-1;
                                                                    }
                                                           }
                                               }
                                              }
                           }
 
                   continue;
             }
             if(!strcmp(type, "seam"))
             {
                   /* fprintf(stderr, "seam \n"); */
                   /*  need to define seam here  */
                   continue;
             }
             if(!strcmp(type, "esource"))
             {
                   /* fprintf(stderr, "esource \n"); */
                   if(OperateFreq != atof(att1)) { fprintf(stderr, "Error: all sources must have the same frequency\n");
                                            exit(1);
                                          }
                   if (x2<x1) swap(&x1, &x2);
                   if (y2<y1) swap(&y1, &y2);
                   if (z2<z1) swap(&z1, &z2);
                   x1=x1-Min_X;y1=y1-Min_Y;z1=z1-Min_Z;
                   x2=x2-Min_X;y2=y2-Min_Y;z2=z2-Min_Z;

                  /* three dimensional source blocks */       
                   if (att2[0] == 'y')
                     {
                      for(i=x1;i<x2;i++){
                       for(j=y1;j<y2;j++){
                        for(k=z1;k<z2;k++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][1]] = 1;
                       EdgeStatus[xn[nn][1]] = atof(att3);
                       ForceStat[xn[nn][3]] = 1;
                       EdgeStatus[xn[nn][3]] = atof(att3);
                       ForceStat[xn[nn][14]] = 1;
                       EdgeStatus[xn[nn][14]] = atof(att3);
                       ForceStat[xn[nn][16]] = 1;
                       EdgeStatus[xn[nn][16]] = atof(att3);
                                           }
                                         }
                                         }
                     }

                   if (att2[0] == 'z')
                     { 
                      for(i=x1;i<x2;i++){
                       for(j=y1;j<y2;j++){
                        for(k=z1;k<z2;k++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][5]] = 1;
                       EdgeStatus[xn[nn][5]] = atof(att3);
                       ForceStat[xn[nn][7]] = 1;
                       EdgeStatus[xn[nn][7]] = atof(att3);
                       ForceStat[xn[nn][10]] = 1;
                       EdgeStatus[xn[nn][10]] = atof(att3);
                       ForceStat[xn[nn][12]] = 1;
                       EdgeStatus[xn[nn][12]] = atof(att3);
                                          }
                                         }
                                        }
                     }

                  if (att2[0] == 'x')
                     {
                      for(i=x1;i<x2;i++){
                       for(j=y1;j<y2;j++){
                        for(k=z1;k<z2;k++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][0]] = 1;
                       EdgeStatus[xn[nn][0]] = atof(att3);
                       ForceStat[xn[nn][4]] = 1;
                       EdgeStatus[xn[nn][4]] = atof(att3);
                       ForceStat[xn[nn][13]] = 1;
                       EdgeStatus[xn[nn][13]] = atof(att3);
                       ForceStat[xn[nn][17]] = 1;
                       EdgeStatus[xn[nn][17]] = atof(att3);
                                          }
                                         }
                                        }
                     }

                  /* two dimensional source blocks */       
                   if ((att2[0] == 'y')&&(z1==z2))
                     {
                      for(i=x1;i<x2;i++){
                       for(j=y1;j<y2;j++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][1]] = 1;
                       EdgeStatus[xn[nn][1]] = atof(att3);
                       ForceStat[xn[nn][3]] = 1;
                       EdgeStatus[xn[nn][3]] = atof(att3);
                                         }
                                         }
                     }
                   if ((att2[0] == 'y')&&(x1==x2))
                     {
                      for(k=z1;k<z2;k++){
                       for(j=y1;j<y2;j++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][1]] = 1;
                       EdgeStatus[xn[nn][1]] = atof(att3);
                       ForceStat[xn[nn][14]] = 1;
                       EdgeStatus[xn[nn][14]] = atof(att3);
                                         }
                                         }
                     }
                   if ((att2[0] == 'z')&&(y1==y2))
                     { 
                      for(i=x1;i<x2;i++){
                        for(k=z1;k<z2;k++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][5]] = 1;
                       EdgeStatus[xn[nn][5]] = atof(att3);
                       ForceStat[xn[nn][7]] = 1;
                       EdgeStatus[xn[nn][7]] = atof(att3);
                                          }
                                        }
                     }
                   if ((att2[0] == 'z')&&(x1==x2))
                     { 
                      for(j=y1;j<y2;j++){
                        for(k=z1;k<z2;k++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][5]] = 1;
                       EdgeStatus[xn[nn][5]] = atof(att3);
                       ForceStat[xn[nn][10]] = 1;
                       EdgeStatus[xn[nn][10]] = atof(att3);
                                          }
                                        }
                     }
                  if ((att2[0] == 'x')&&(y1==y2))
                     {
                      for(i=x1;i<x2;i++){
                        for(k=z1;k<z2;k++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][0]] = 1;
                       EdgeStatus[xn[nn][0]] = atof(att3);
                       ForceStat[xn[nn][13]] = 1;
                       EdgeStatus[xn[nn][13]] = atof(att3);
                                          }
                                        }
                     }
                  if ((att2[0] == 'x')&&(z1==z2))
                     {
                      for(i=x1;i<x2;i++){
                        for(j=y1;j<y2;j++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][0]] = 1;
                       EdgeStatus[xn[nn][0]] = atof(att3);
                       ForceStat[xn[nn][4]] = 1;
                       EdgeStatus[xn[nn][4]] = atof(att3);
                                          }
                                        }
                     }

                  /* one dimensional source blocks */       
                   if ((att2[0] == 'x')&&(z1==z2)&&(y1==y2))
                     {
                       for(i=x1;i<x2;i++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][0]] = 1;
                       EdgeStatus[xn[nn][0]] = atof(att3);
                                         }
                     }
                   if ((att2[0] == 'y')&&(z1==z2)&&(x1==x2))
                     {
                       for(j=y1;j<y2;j++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][1]] = 1;
                       EdgeStatus[xn[nn][1]] = atof(att3);
                                         }
                     }
                   if ((att2[0] == 'z')&&(y1==y2)&&(x1==x2))
                     {
                       for(k=z1;k<z2;k++){
                       nn=k*XdiM*YdiM+j*XdiM+i;
                       ForceStat[xn[nn][5]] = 1;
                       EdgeStatus[xn[nn][5]] = atof(att3);
                                         }
                     }

                 continue;
              }
              if(!strcmp(type, "celldim"))
              {continue;
              }
           }              /*end of if*/ 
        }                 /* end of while */
        rewind(InF);
        remove("source_nodes.tmp");
        CountEdgeType(); 
}

/***************************************************************************
*   FUNCTION        : absol()
****************************************************************************/

int    absol(x)
        int x; { if (x<0) return(-x); else  return(x); }

/***************************************************************************
*
*   FUNCTION    : swap(x,y)
*   DESCRIPTION : swap the values of x and y
*
****************************************************************************/

void swap(x,y)
        int *x,*y; { int temp; temp = *x; *x=*y; *y=temp; }

/***************************************************************************
*   FUNCTION        : CountEdgeType()
****************************************************************************/

void CountEdgeType()
{

 TotInnerEdgeNum = 0;
 TotBoundEdgeNum = 0;
 TotForcdEdgeNum = 0;
 
 for(IEdge=0;IEdge<=TotEdgeNum-1;IEdge++)
 {
 if(ForceStat[IEdge] == 0)   {
    if(EdgeStatus[IEdge]==-1)       {TotInnerEdgeNum++;}
    else if(EdgeStatus[IEdge]==-2)  {TotBoundEdgeNum++;}
                         }
 else                    {           TotForcdEdgeNum++;}
  }

 InnerEdgeStat = INT_Vector(TotInnerEdgeNum);
 BoundEdgeStat = INT_Vector(TotBoundEdgeNum);
 ForcdEdgeStat = INT_Vector(TotForcdEdgeNum);
 ForcdValue    = FLOAT_Matrix(TotForcdEdgeNum,2);

 TotInnerEdgeNum = 0;
 TotBoundEdgeNum = 0;
 TotForcdEdgeNum = 0;

 for(IEdge=0;IEdge<=(TotEdgeNum)-1;IEdge++)
 {
     if(ForceStat[IEdge] == 0)   {
        if(EdgeStatus[IEdge]==-1)                   {
        InnerEdgeStat[TotInnerEdgeNum]= IEdge;
        TotInnerEdgeNum++;
                                        }
        else if(EdgeStatus[IEdge]==-2)              {
        BoundEdgeStat[TotBoundEdgeNum]= IEdge;
        TotBoundEdgeNum++;
                                        }
                       }
     else              {
        ForcdEdgeStat[TotForcdEdgeNum]= IEdge;
        ForcdValue[TotForcdEdgeNum][0] = EdgeStatus[IEdge];
        ForcdValue[TotForcdEdgeNum][1] = 0.0;
        TotForcdEdgeNum++;
                           }
       }
for(IEdge=0;IEdge<TotNumHexHedron;IEdge++)
{
RELPerm[IEdge][0] = Epsilon[IEdge];
RELPerm[IEdge][1] = Sigma[IEdge];
}
printf("TotInnerEdgeNum = %d\n",TotInnerEdgeNum);
printf("TotBoundEdgeNum = %d\n",TotBoundEdgeNum);
printf("TotForcdEdgeNum = %d\n",TotForcdEdgeNum);
}

/***************************************************************************
*   FUNCTION        : edge()
****************************************************************************/

void edge(nn,i,j,k)
int nn,i,j,k;
{
int kk,xx,dy1,dy2,divz,divz1;
   
    for (kk=0;kk<18;kk++) 
       xn[nn][kk] = 0;
    divz=(6*XdiM*YdiM+3*XdiM+3*YdiM+1)*k;
    divz1=XdiM*YdiM*6 + XdiM*3 +YdiM*3 + 1 ;
    dy1 = (3*XdiM+1)*j; dy2 = (3*XdiM+2)*j;
    xn[nn][0] = i + dy1 + divz;
    xn[nn][1] = XdiM + 2*i + dy1 + divz;
    xn[nn][2] = xn[nn][1] + 1; 
    xn[nn][3] = xn[nn][2] + 1;
    xn[nn][4] = 3*XdiM + 1 + i + dy1 + divz;
    xx = (YdiM+1)*XdiM + (2*XdiM+1)*YdiM;
    xn[nn][5] = xx + 2*i + dy2 + divz;
    xn[nn][6] = xn[nn][5] + 1; 
    xn[nn][7] = xn[nn][6] + 1;
    xx = xx + (2*XdiM+1); 
    xn[nn][8] = xx + i + dy2 + divz;
    xn[nn][9]= xn[nn][8] + 1;
    xx = xx + XdiM + 1; 
    xn[nn][10]= xx + 2*i + dy2 + divz;
    xn[nn][11]= xn[nn][10] + 1;
    xn[nn][12]= xn[nn][11] + 1; 
    xn[nn][13]= xn[nn][0] + divz1;
    xn[nn][14]= xn[nn][1] + divz1; 
    xn[nn][15]= xn[nn][2] + divz1;
    xn[nn][16]= xn[nn][3] + divz1; 
    xn[nn][17]= xn[nn][4] + divz1;
  return;
}

/***************************************************************************
*   FUNCTION        : deter_edge()
****************************************************************************/

void deter_edge(nn,i,j,k)
int nn,i,j,k;
{
int scratch;

     scratch=xn[nn][0];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 1\n",scratch);  fflush(stdout);
        }
     if ((j==0)||(k==0))   EdgeStatus[scratch] = -2;
     else                  EdgeStatus[scratch] = -1;

     scratch=xn[nn][1];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 2\n",scratch); fflush(stdout);
        }
     if ((i==0)||(k==0))   EdgeStatus[scratch] = -2;
     else                  EdgeStatus[scratch] = -1;

     scratch=xn[nn][2];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 3\n",scratch); fflush(stdout);
         fflush(stdout);
       }
     if (k==0)            EdgeStatus[scratch] = -2;
     else                 EdgeStatus[scratch] = -1;            
           
     scratch = xn[nn][3];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 4\n",scratch); fflush(stdout);
       }
     if ((i==(XdiM-1))||(k==0))    EdgeStatus[scratch] = -2; 
     else                          EdgeStatus[scratch] = -1;

     scratch = xn[nn][4];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 5\n",scratch); fflush(stdout);
       }
     if ((j==(YdiM-1))||(k==0))   EdgeStatus[scratch] = -2; 
     else                         EdgeStatus[scratch] = -1; 

     scratch = xn[nn][5];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 6\n",scratch); fflush(stdout);
       }
     if ((j==0)||(i==0))   EdgeStatus[scratch] = -2;
     else                  EdgeStatus[scratch] = -1;
 
     scratch = xn[nn][6];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 7\n",scratch); fflush(stdout);
       }
     if (j==0)   EdgeStatus[scratch] = -2;
     else        EdgeStatus[scratch] = -1;

     scratch = xn[nn][7];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 8\n",scratch); fflush(stdout);
       }
     if ((i==(XdiM-1))||(j==0))  EdgeStatus[scratch] = -2;
     else                        EdgeStatus[scratch] = -1;

     scratch = xn[nn][8];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 9\n",scratch); fflush(stdout);
       }
     if (i==0)    EdgeStatus[scratch] = -2;
     else         EdgeStatus[scratch] = -1;

     scratch = xn[nn][9];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 10\n",scratch); fflush(stdout);
       }
     if (i==(XdiM-1))   EdgeStatus[scratch] = -2;
     else               EdgeStatus[scratch] = -1;
 
     scratch = xn[nn][10];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 11\n",scratch); fflush(stdout);
       }
     if ((j==(YdiM-1))||(i==0))  EdgeStatus[scratch] = -2;
     else                        EdgeStatus[scratch] = -1;

     scratch = xn[nn][11];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 12\n",scratch); fflush(stdout);
       }
     if (j==(YdiM-1))    EdgeStatus[scratch] = -2;
     else                EdgeStatus[scratch] = -1;
 
     scratch = xn[nn][12];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 13\n",scratch); fflush(stdout);
       }
     if ((i==(XdiM-1))||(j==(YdiM-1))) EdgeStatus[scratch] = -2;
     else                              EdgeStatus[scratch] = -1;
 
     scratch = xn[nn][13];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
        {
          printf("scratch=%d out of bounds spot 14\n",scratch); fflush(stdout);
        }
     if ((k==(ZdiM-1))||(j==0))   EdgeStatus[scratch] = -2;
     else                         EdgeStatus[scratch] = -1;
 
     scratch = xn[nn][14];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
        {
          printf("scratch=%d out of bounds spot 15\n",scratch); fflush(stdout);
        }
     if ((k==(ZdiM-1))||(i==0))    EdgeStatus[scratch] = -2;
     else                          EdgeStatus[scratch] = -1;
 
     scratch = xn[nn][15];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
        {
          printf("scratch=%d out of bounds spot 16\n",scratch); fflush(stdout);
        }
     if (k==(ZdiM-1))  EdgeStatus[scratch] = -2;
     else              EdgeStatus[scratch] = -1;
 
     scratch = xn[nn][16];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
       {
         printf("scratch=%d out of bounds spot 17\n",scratch); fflush(stdout);
       }
     if ((i==(XdiM-1))||(k==(ZdiM-1)))  EdgeStatus[scratch] = -2;
     else                               EdgeStatus[scratch] = -1;

     scratch = xn[nn][17];
     if ((scratch >= TotEdgeNum) || (scratch < 0))
        {
          printf("scratch=%d out of bounds spot 18\n",scratch); fflush(stdout);
        }
     if ((j==(YdiM-1))||(k==(ZdiM-1)))  EdgeStatus[scratch] = -2;
     else                               EdgeStatus[scratch] = -1;
  return;
 }

/***************************************************************************
*   FUNCTION        : Produce_Output()
****************************************************************************/

void Produce_Output()
{
 int    i, j, k, x1, y1, z1, x2, y2, z2, Count_i, dflag=0;
 char   type[20], att1[20], att2[8], att3[8], att4[8];
 char   buffer[80];
 FILE   *OutF;

   EfieldData = (double *) calloc( TotEdgeNum*sizeof(double), sizeof(double));
   for(Count_i=0;Count_i<TotInnerEdgeNum;Count_i++)
   EfieldData[InnerEdgeStat[Count_i]] = RHSVector[Count_i];
   for(Count_i=0;Count_i<TotForcdEdgeNum;Count_i++)
   EfieldData[ForcdEdgeStat[Count_i]] = ForcdValue[Count_i][0];

        while (fgets(buffer, 80, InF))
        {
         if (buffer[0] == '#') {
                               if(OutF_0 != NULL) fprintf(OutF_0, "%s", buffer);
                               if(OutF_1 != NULL) fprintf(OutF_1, "%s", buffer);
                               if(OutF_2 != NULL) fprintf(OutF_2, "%s", buffer);
                               if(OutF_3 != NULL) fprintf(OutF_3, "%s", buffer);
                               continue;
                              }
        if(!strncmp(buffer, "default_output",14)) {dflag=1; continue;}
        if(!strncmp(buffer, "efield_output",13))
                   { 
                   sscanf(buffer, "%s%d%d%d%d%d%d%s%s%s%s", type, &x1, &y1, &z1, &x2, &y2, &z2, att1, att2, att3, att4);
                   if (x2<x1) swap(&x1, &x2);
                   if (y2<y1) swap(&y1, &y2);
                   if (z2<z1) swap(&z1, &z2);
                   if (!strcmp(att1, Out_FileName1)) OutF=OutF_1;
                   if (!strcmp(att1, Out_FileName2)) OutF=OutF_2;
                   if (!strcmp(att1, Out_FileName3)) OutF=OutF_3;
                   for(k=z1;k<=z2;k++){
                     for(j=y1;j<=y2;j++){
                        for(i=x1;i<=x2;i++){
                                            EfieldatNode(i, j, k, OutF); 
                                           }
                                        }
                                      } 
                   continue;
                   }

        }                 /* end of while */
        if(dflag == 1)    /* prints default output data */
              {
               for(Count_i=0;Count_i<TotEdgeNum;Count_i++){
               fprintf(OutF_0,"%d\t%f\n",Count_i,EfieldData[Count_i]);}
               fclose(OutF_0);
              }
        fclose(OutF_1);
        fclose(OutF_2);
        fclose(OutF_3);
        fclose(InF);
        fprintf(stdout,"\nProgram Successfully Completed: \n\n");
}


/***************************************************************************
*   FUNCTION        : EfieldatNode()
****************************************************************************/

void EfieldatNode(X,Y,Z,Out_F)
int X,Y,Z;
FILE *Out_F;
{
int INode,JNode,IEdge,NodeNum[3][2],CountX=0,CountY=0,CountZ=0;
double Xcord,Ycord,Zcord;
double ValX,ValY,ValZ;

Xcord = DivisorX*(double)X;
Ycord = DivisorY*(double)Y;
Zcord = DivisorZ*(double)Z;
 
for(JNode=0;JNode<TotNodeNum;JNode++)
 {
 if((NodeCord[JNode][0]==Xcord) && (NodeCord[JNode][1]==Ycord) && (NodeCord[JNode][2]==Zcord))
 {INode = JNode; break;} else {continue;}
 }
 
   for (IEdge=0;IEdge<TotEdgeNum;IEdge++)
    {
 
    if(((TetGlobalEdgeEnds[IEdge][0])==INode)
    || ((TetGlobalEdgeEnds[IEdge][1])==INode) )
    {
     if((NodeCord[TetGlobalEdgeEnds[IEdge][0]-1][1]==NodeCord[TetGlobalEdgeEnds[IEdge][1]-1][1]) &&
        (NodeCord[TetGlobalEdgeEnds[IEdge][0]-1][2]==NodeCord[TetGlobalEdgeEnds[IEdge][1]-1][2]))
        {NodeNum[0][CountX] = IEdge; CountX++;}
 
      else if((NodeCord[TetGlobalEdgeEnds[IEdge][0]-1][2]==NodeCord[TetGlobalEdgeEnds[IEdge][1]-1][2]) &&
              (NodeCord[TetGlobalEdgeEnds[IEdge][0]-1][0]==NodeCord[TetGlobalEdgeEnds[IEdge][1]-1][0]))
              {NodeNum[1][CountY] = IEdge; CountY++;}
 
      else if((NodeCord[TetGlobalEdgeEnds[IEdge][0]-1][0]==NodeCord[TetGlobalEdgeEnds[IEdge][1]-1][0]) &&
              (NodeCord[TetGlobalEdgeEnds[IEdge][0]-1][1]==NodeCord[TetGlobalEdgeEnds[IEdge][1]-1][1]))
              {NodeNum[2][CountZ] = IEdge; CountZ++;}
      else
              continue;
    }
   }
 
  ValX = (EfieldData[NodeNum[0][0]]+EfieldData[NodeNum[0][1]])/2.0;
  ValY = (EfieldData[NodeNum[1][0]]+EfieldData[NodeNum[1][1]])/2.0;
  ValZ = (EfieldData[NodeNum[2][0]]+EfieldData[NodeNum[2][1]])/2.0;

  fprintf(Out_F,"%6.3f\t%6.3f\t%6.3f\t%f\t%f\t%f\n",Xcord,Ycord,Zcord,ValX,ValY,ValZ);
}

/*************************** END OF PROGRAM ********************************/