sift5.c

三维矢量有限元-矩量法电磁场分析程序。 EMAP5 is a full-wave electromagnetic field solver that combines the method of m

/*****************************************************************************
 *  The following codes read a sif file as the input file, then generate     *
 *  the mesh data, source information and output requirements for emap5.     *
 *                        						                 *
 *  Authors:  Dr. Mohammad W. Ali,   Yun  Ji,  Dr. Todd. H. Hubing           *
 *  Version:  1.1							                       *
 *  Last updated:  Dec 16, 1998   				                       *
 *****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "complex.c"
#include "alloc.c"
#include "util.c"

typedef struct {
    int x1,y1,z1,x2,y2,z2;   /* two corner points of the conductor */
    double   DelX, DelY, DelZ;  /* step size */
} Cond;


typedef struct{
    double  start;  /* start point */
    double  end;    /* end point */
    double  step;   /* step */
} CellUnit;


typedef struct{
    int x;        
    int y;
} XPoint;



typedef struct{
   XPoint start;
   XPoint end;
}  RectAngle;



/******************   FUNCTION PROTOTYPE  ************************/
void AllocateMemory();
void AssignBoxCoordinates();
void AssignHexHedronEdgeNum();
void AssignHexHedronNodeNum();
void ComputeEdgeMultiplicity();
void ComputeParameters();
void ComputeTrngleCentroid(int, double **, int **, double **);
void DetectBoundEdgeStatus();
void DetectInnerEdgeStatus();
void DetectISourceEdges();
void DetectVSourceEdges();
void DiscretizeFEMVolume();
void DiscretizeExterPECPatch();
void FindInterMetalEdges();
void FindTetraHedronFaceNodes();
void GlobalEdgeEndsDiv1();
void GlobalEdgeEndsDiv2();
int  Intersection(RectAngle, RectAngle, RectAngle *, int *, RectAngle *);
int  Max(int, int);
int  Min(int, int);
void ModifyBoundFaces(int );
void PrintOutput();
void ReadSIFFile();
void TetraEdgeDiv1();


/******************** GLOBAL VARIABLES **********************/

int  DetectCondBound(Cond);
void DetectInnerCond(Cond);


CellUnit CellUnitX[10], CellUnitY[10], CellUnitZ[10];
int   CellUnitXNum=0;
int   CellUnitYNum=0;
int   CellUnitZNum=0;

Cond conductor[20];
int CondNum;


int flag_default_out=0; /* flag used to determine which kind of
                                           output users want */

int output_num;  /* how many times "output" are used, maxium number is 10 */ 
int output_x[2][10], 
    output_y[2][10], 
    output_z[2][10]; /* store information for  each "output" */

char output_file[10][20];  /* file to store the output data */
char output_axis[10];  /* only edges parallel to the designated axis 
                              wil be printed out */

char default_out_file[20]; /* file to store default output data */

double E_Theta, E_Phi, K_Theta, K_Phi, PlaneWaveE_Mag;


int JunctionNum=0;
/* maxium junction edge allowed is 100. 
   Junction[][0] ---- the hybrid edge number 
   Junction[][1] ---- the other external metal edge with same coordinates 
                      with Junction[][0]. these two form a juction      
   Junction[][2]---- the T++ or T-- triangle associate with these Junction */

int Junction[100][3];



/*  store the resistor edges */
int ResisterLocX[20][2];
int ResisterLocY[20][2];
int ResisterLocZ[20][2];
int ResisterNum;
double ResisterValue[20];


char  In_FileName[50];

int
    **HexNode,
    **HexEdgeNum,
    **TetTrngleNode,
    *TetBoundFaceTetHedronNum,
    *TetBoundFaceFaceNum,
    **TrngleEdge,
    **TrngleNode,
    **TetNode,
    **FaceNode,
    **FaceEdge,
    *EdgeFound,
    *MultFound,
    **PlusTrngleDetect,
    *PlusTrngleIndex,
    **MinusTrngleDetect,
    LastMultEdgeFoundNum,
    LastEdgeFoundNum,
    MetalInEdgeNum,
    MetalBoundEdgeNum,
    DipolWid[10],DipolLen[10],
    DielBoundEdgeNum,
    ExtEdge1,ExtEdge2,
    ExtEdge3,ExtEdge4,
    ExtEdge5,ExtEdge6,
    ImpEdge1,ImpEdge2,
    InnEdge1,InnEdge2,
    *MinusTrngleIndex,
    **TetEdge,
    *EdgeStat,
    *ISourceEdgeStat,
    *ForcdStat,
    TotISourceEdgeNum,
    *MetalEdgeStat,
    *InnerEdgeStat,
    *BoundEdgeStat,
    **GlobalEdgeEnds
;


int PatchNum,FrontWall,BackWall,BottomWall;
int TopWall,RightWall,LeftWall;

 
double
    VsourceMag,
    OperateFreq,
    CellDimension,
    **MediaPerm,
    **TrngleCenTroid,
    *TrngleArea,
    **TrngleNormal,
    ** NodeCord,
    *ISourceEdgeMag;



int CubeTrngleNum,CubeNumNode,CubeNumEdge,
    FrontEdge,FrontNode,BackEdge,BackNode;

int 
    IEdge,
    JEdge,
    VVal,II,KK,
    TotBoundEdgeNum,
    FExtDipolLen,
    BExtDipolLen,
    TotTrngleNum,
    TotDipolTrngleNum[10],
    TotCubeTrngleNum,
    HybrdBoundEdgeNum,
    TotMetalTrngleNum,
    ExtDipolNodeStart,
    ExtDipolEdgeStart,
    ExtDipolTrngleStart,
    TotInnerEdgeNum,
    Xcount,Ycount,Zcount,
    PLen,
    start, TotNumHexHedron,
    GlobalTetFaceNum,
    LocalTetFaceNum,
    i,j,k,ii,jj,kk,Count,
    TotTetFaceNum=0,
    XdiM,YdiM,ZdiM,
    TotEdgeNum,
    TotDipolNumEdge[10],
    TotNodeNum,
    TotDipolNumNode[10],
    TotCubeNumEdge,
    TotCubeNumNode,
    HexNum,con,
    PermConst,
    nnode,
    TotTetElement;

int 
    XBoundMinDim,YBoundMinDim,ZBoundMinDim,
    XBoundMaxDim,YBoundMaxDim,ZBoundMaxDim,
    XDielmat[2][10],YDielmat[2][10],ZDielmat[2][10],
    XintEsrc[2][10],YintEsrc[2][10],ZintEsrc[2][10],
    XextVsrc[2][10],YextVsrc[2][10],ZextVsrc[2][10] ,
    Xintcond[2][10],Yintcond[2][10],Zintcond[2][10],
    Xbndcond[2][10],Ybndcond[2][10],Zbndcond[2][10],
    Xextcond[2][10],Yextcond[2][10],Zextcond[2][10],
    XmaxWall,YmaxWall,ZmaxWall,
    XminWall,YminWall,ZminWall,
    DelXintcond[10],DelYintcond[10],DelZintcond[10],
    DelXbndcond[10],DelYbndcond[10],DelZbndcond[10],
    DelXextcond[10],DelYextcond[10],DelZextcond[10],
    XBoundDim[2],YBoundDim[2],ZBoundDim[2],
    VSourceEdge[20],
    VSourceNum=0,
    DielmatNum=0,
    IntEsrcNum=0,
    ExtVsrcNum=0,
    IntcondNum=0,
    BndcondNum=0,
    ExtcondNum=0,
    Max_X=0,
    Max_Y=0,
    Max_Z=0,
    Min_X=0,
    Min_Y=0,
    Min_Z=0;


double  VSourceMag[20],ExtVsrcMag[10] ;


char SourceType;
char   IntEsrcAxis[10];   /* store E source direction */
double IntEsrceMag[10];   /* store Esource magnetude*/


char PolInc[1];
double ThetaInc,PhiInc;
double Epsilon[10],Sigma[10];


int Dipolplane_XY[10],Dipolplane_YZ[10],Dipolplane_ZX[10];

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

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



int main(int argc, char **argv)
{
    char buff[10];

    if( argc != 3 ) 
    {
        fprintf(stderr,"Usage: sift5 inputfile outputfile\n");
        exit(1);
    }

    if (strlen(argv[1])<=4) 
    {
       fprintf(stderr, "Error: Use *.sif as the input file\n");
       exit(1);
    }
       
    /* find out the surfix of the input file */
    strcpy(buff, argv[1]+strlen(argv[1])-4 );
    if (strcmp(buff, ".sif")!=0) 
    {
       fprintf(stderr, "Error: Use *.sif as the input file\n");
       exit(1);
    }

    if (strlen(argv[2])<=3) 
    {
       fprintf(stderr, "Error: Use *.in as the output file\n");
       exit(1);
    }
       
    /* find out the surfix of the output file */
    strcpy(buff, argv[2]+strlen(argv[2])-3 );
    if (strcmp(buff, ".in")!=0) 
    {
       fprintf(stderr, "Error: Use *.in as the input file\n");
       exit(1);
    }
      
    if( (InF = fopen(argv[1], "r"))== NULL ) 
    {
        fprintf(stderr, " Error: Can't open the sif file.\n");
        exit(1);
    }
   

    if ( (OutFilePntr = fopen(argv[2],"w"))==NULL) 
    {
         fprintf(stderr, " Error: Can't open output file.\n");
	 exit(1);
    }

 
    ReadSIFFile();
    ComputeParameters();
    AllocateMemory();
    AssignBoxCoordinates();
    AssignHexHedronNodeNum();
    AssignHexHedronEdgeNum();

    DetectInnerEdgeStatus();
    DetectBoundEdgeStatus();

    DiscretizeFEMVolume();

    DiscretizeExterPECPatch();
 
    FindInterMetalEdges();  
    
    ComputeEdgeMultiplicity();


    /* initialize ForcdStat[] */
    ForcdStat = INT_Vector(TotEdgeNum);

    for(ii=0; ii<TotEdgeNum; ii++)
        ForcdStat[ii]=0;   /* all are set as "not forcd" */

    ISourceEdgeStat=INT_Vector(TotEdgeNum);

    /* initilize ISourceEdgeMag */
    ISourceEdgeMag=double_Vector(TotEdgeNum);

    for(ii=0;ii<TotEdgeNum; ii++)
        ISourceEdgeMag[ii]=0.0;

  
    if( SourceType=='I' ) DetectISourceEdges();
   
    /* to build the inner-edge-table */
    kk=0;
    for(ii=0;ii<TotEdgeNum;ii++)
        if(EdgeStat[ii]==1)  /* it is an internal edge */
        {  
           InnerEdgeStat[kk]=ii+1; 
           kk++;   
        }

  
    /* the inner edges should not include those forcd edges, 
       thus revise  TotInnerEdgeNum here */
    TotInnerEdgeNum=kk;  
   

    kk=0;
    HybrdBoundEdgeNum=0;
    DielBoundEdgeNum=0;
    BoundEdgeStat = INT_Vector(TotBoundEdgeNum);


    /* find out boundary edges, those metal edges forming the contour
       of external conductors, which are charactered by MetalEdgeStat[i] 
       equaling to 3, are excluded  since J on this edges
       are equal to zero. */

    /* the following are dielectric boundary edges 
       EdgeStat[ii]==-1 means boundary edge 
       MetalEdgeStat[ii]==0 means it is not a metal patch, 
       instead it is a dielectric edge  */
    for(ii=0;ii<TotEdgeNum;ii++)
        if( (EdgeStat[ii]==-1) && (MetalEdgeStat[ii]==0) )