emap5.c

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

/***************************************************************************
 *                                                                         *
 *  The EMAP5 hybrid FEM/MoM 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 H. Hubing,                                               *
 *       Department of Electrical and Computer Engineering                 *
 *       University of Missouri-Rolla, 					         *
 *       Rolla, MO  65409.                                                 *
 *       E-mail:  [e-mail]					         *
 *                                                                         *
 *  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.                                          *
 *                                                                         *
 *   PROGRAM     :  emap5.c						               *
 *   VERSION     :  1.1							               *
 *   AUTHORS     :  Dr. Mohammad W. Ali, Yun Ji, Dr. Todd T. Hubing        *
 *   LAST UPDATE :  Dec 16,  1998 					               *
 *   DESCRIPTION :  A 3-D hybrid FEM/MoM Code for Analyzing                *
 *                  Time Varying Complex Electromagnetic Fields.           *
 *									                     *
 ***************************************************************************/


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

/***************************  MARCRO DEFINITION ***************************/

/* TOLERANCE can be changed to 1E-3 or 1E-4  */
#define  TOLERANCE 1E-4

#ifndef M_PI
#define M_PI 3.1415926
#endif 

/****** Don't change the following four marcros ******/

/*  Use 7 point Gaussian Quadrature to evaluate elements of MOM matrices */
#define  TotQuadPoint 7

/* the maximum non-zero elements of B matrix */
#define  MaxBmatElementNum 5

/*  the maximum non-zero elements of A matrix */ 
#define  MaxAmatElementNum 30

/* the maximum forced edge number */
#define  MaxFORCDMatElementNum 15




/*************************** FUNCTION PROTOTYPE ****************************/

void ComputeBddMatrix(int, int, double *, double *, int *, int *,
                      int *, int *, double *, double [][3]);
void ComputeCMatrix(int, int,  double *, double *, int *, int *,
                    int *, int *, double **);
void ComputeCorrectionTerm();
void ComputeDMatrix(int, int, double *, double *, int *, int *,
                    int *, int *, double **, double [][3]);
void ComputeGaussQuadPoint(int, int *, double *);
void ComputeParameters();
void ComputeSingul(int *, double *, double, double*, double *, double *); 
void ComputeSourceVector( double, double *, complex *);
void ComputeNonSingul(int *, double, double *,complex *, 
                      complex *, complex *, complex *);
void ComputeNonSingul1(int *, double *, complex *, complex *, 
                       complex *, complex *);
void ComputeTetraHedronVolume(int,double **,int **,double *);
void ComputeTrngleCentroid(int, double **, int **, double **);
void ConjugateSolver(int,int *, int **, complex **, complex *,int, double);
void CreateRHSVector(int ,int, int, double *, complex *);
void FemMatrixCompute(int,int, int **, int **,double **, complex *, 
                      complex *, int **, complex **, int *, double *);
complex InnerProduct(complex *, complex *,int );
void InvertMatrix(complex **, int);
void matrixVectorproduct(char, int, complex *, complex *, 
                         int i, int **, complex **);
int  SearchNonZeroElement(int, int);
void SurfaceFieldCompute();
void TetFaceAreaNormal(int,double **,int **, double **,double *);
void PrintOutput();
void ReadInputFile();
double Sign(int);
double VectorNorm(complex *,int);



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



/****************************************************************************
       the following varialbes are used to represent resistor edges 
*****************************************************************************/
/* to represent */
int ResistorEdgeStat[100];      /* the resistor edge index */
double ResistorEdgeValue[100];  /* the resistor value */
int ResistorEdgeNum;            /* total number of resistor edges */


/****************************************************************************
       the following varialbes are used to represent junctions 
*****************************************************************************/
int JunctionNum=0;     /* total number of junctions, initialize to zero */

/* 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 the junction  */

int Junction[100][3];  /* maxium junction edge allowed is 100. */


char SourceType;   /* Source type. 
                      'V' means voltage source
                      'P' means incident plane wave source
                      'I' meams current source  */

/****************************************************************************
       the following varialbes are used to represent voltage sources 
 ****************************************************************************/
int VSourceNum,         /* voltage source edge number */
    VSourceEdge[20];    /* voltage source edge table. The maximum
                           voltage edges allowed is 20  */
double VSourceMag[20];  /* voltage edge magnitude */


/****************************************************************************
       the following varialbes are used to represent plane wave 
 ****************************************************************************/
double K_Theta, K_Phi,  /* the unit vector of K in the spherical  
                           coordiantes   */

       E_Theta, E_Phi,  /* the unit vecotr of E in the sperical  
                                    coordinates */
       PlaneWaveE_Mag;  /* the E field magnitude of the plane wave */


int 

    DielBoundEdgeNum,     /* total number of the diectric boundary edges */
    HybrdBoundEdgeNum,    /* total number of the hybrid boundary edges */
    TotEdgeNum,           /* total number of edges */
    TotTetElement,        /* total number of tetrahedrons */
    TotNodeNum,           /* total number of nodes */
    TotTrngleNum,         /* total number of triangles */
    TotBoundEdgeNum,      /* total number of boudary edgs */
    TotInnerEdgeNum,      /* total number of inner edgs */
    TotExtMetalEdgeNum,   /* total number of the exteral metal-boundary edges */

    TotISourceEdgeNum,
    TotMatrixEqnNum,      /* the size of final matrix equation after MoM is  
                             coupled to FEM matrix equation */

    **TrngleNode,      /* triangle node table */
    **TrngleEdge,      /* triangle edge table */

    *BoundEdgeStat,       /* boudary edge table */
    *InnerEdgeStat,       /* inner edge table */
    *ISourceEdgeStat,       /* the current source edge table */

    **PlusTrngleDetect,      /* plus triangle detected table */
    *PlusTrngleIndex,     /* plux triangle table */
    **MinusTrngleDetect,     /* minus triangle detected table */
    *MinusTrngleIndex,    /* minus triangle table */
    *VcRowNum ,
    **GBLmatColAddr,      /* the column table for each row in FEM matrix, which
                             is used by the row-indexed scheme */
    *GBLmatColIndex,      /* the number of non-zero element in each row in FEM
                             matrix, which is used by the row-indexed scheme */

    **TetNode,   /* tetrahedron node table */
    **TetEdge,   /* tetrahedron edge table */

    **GlobalEdgeEnds;  /* glabal edge table */


double

    OperateFreq,    /* the working frequency, the unit is Hz */
    FreeSpaceVel,   /* the free space light velocity: 3X10^8 m/s */
    AbsPermeable,   /* the absolute permiability: 4*Pi*10^-7  */
    AbsPermitt,     /* the absolute permittivity: 8.854*10^-12 */
    WaveNumber,     /* wavenumber, defined as wavelength/(2*pi) */
    WaveLength,     /* wavelength in meters */
    ImpeDance,      /* the impedance of free space, which equals to 377 ohms */
    *Wght,          /* the weighting coefficients of Gassian Qadrature */
    **Qpnt,         /* the position coefficients of Gassian Qadrature */
    *EdgeLength,    /* edge length table */
    *TetVolume,     /* tetrahedron volume table */
    **NodeCord,     /* global node table */
    **TrngleNormal, /* the normal vector of each triangle. TrngleNormal[i][0],
                       TrngleNormal[i][1] and TrngleNormal[i][2] are the x, y
                       and z component of the normal vector of triangle i+1. */

    **TrngleCenTroid,  /* the centroid point of each triangle.
                          TrngleCenTroid[i][0], TrngleCenTroid[i][1], 
                          and TrngleCenTroid[i][2] are the x, y and z 
                          coordinates of the centroid of triangle i+1 */


    *TrngleArea,    /* the triangle area table. TrngleArea[i] is the
                       area of triangle i+1 */
    **Bdd,          /* the surface term in the FEM equation. */
    *ISourceEdgeMag;   /* the current source value */


complex

    CFactor1,       /* one constant associated with MoM equation, see
                       ComputeParameters() */
    CFactor2,       /* one constant associated with MoM equation, see
                       ComputeParameters() */
    TParam[2],      /* two constants associated with FEM equation, see
                       ComputeParameters() */
    *Epsilon,       /* the complex permittivity table of each tetrahedron */
    *RHSVector,     /* the Right-Hand Side vector (RHS) of the 
                       final hybrid matrix equation */
    **Cdd, **Cdc, **Ccd, **Ccc,   /* four submatices of the MoM C matrix */
    *Einc,          /* the incident wave term generated by the testing
                       functions.  */  
    **Ddd, **Dcd,   /* the two submatrices of the MoM D matrix */
    **GBLmat,       /* the FEM matrix data stored in row-indexed scheme */


    /*********************************************************************
       the following six vectors are used in ComputerCorrectionTerm().
       They are used to couple the MoM equation to the FEM equation
     *********************************************************************/
    *GdVector,
    *VcVector,
    *VdVector, 
    *UdVector,  
    *MdVector,     
    *McVector,   

    *EdVector,   /* the magnetic current density on the dielectric surface */
    *JcVector,   /* the electric current density on the conductive surface */
    *JdVector;   /* the electric current density on the dielectric surface */


time_t start_time, mid_time, end_time;  /* to count CPU time used */
char  log_mesg[300];     /* store the log information generated in 
                            ConjugateSolver( ) */

FILE *InF,*OutF,*OutF1;  /* InF: file pointer of the input file 
                            OutF: file pointer of the output file
                            OutF1: file pointer of the log file  */



/**************************************************************************
 *		 	   MAIN PROGRAM 				  *
 **************************************************************************/
int main(int argc, char **argv)
{
    int IEdge, JEdge, ObservTrngle, SourceTrngle,     
        ObservEdge, QuadPointM, QuadPointN;
    char filename[40];
    double **SrcPointArray, **ObsPointArray, ***RhoCtrd;


    if( argc!=2 )
    {
        fprintf(stderr, "Usage:emap5 inputfile\n");
        exit(1);
    }
 
    if( strlen(argv[1])<=3 ) 
    {
        fprintf(stderr, "please use *.in  as the input file\n");
        exit(1);
    }

    /* get the surfix */
    strcpy( filename, argv[1]+strlen(argv[1])-3 );

    if( strcmp(filename, ".in")!=0 ) 
    {
        fprintf(stderr, "please use .in surfix for the input file\n");
        exit(1);
    }

    /* open the input file */
    if( (InF = fopen(argv[1],"r"))==NULL )
    {
        fprintf(stderr, "Sorry, can't open input file\n");
        exit(1);
    }
 

    /* open the log file */
    strncpy( filename, argv[1], strlen(argv[1])-3 );
    filename[ strlen(argv[1])-3 ] ='\0';    /* add numm character */

    strcat(filename, ".log");
 
    if( (OutF1=fopen(filename, "w"))==NULL )
    {
        fprintf(stderr, "Sorry, can't open log file\n");
        exit(1);
    } 

    ReadInputFile();

    if (SourceType=='V')
    {     
        int flag=0;

        for(IEdge=0;IEdge<VSourceNum;IEdge++)
        {
        
 	    for(JEdge=0;JEdge<TotBoundEdgeNum;JEdge++)  
            {
                if( BoundEdgeStat[JEdge]==VSourceEdge[IEdge] ) 
                {

                    /* vsource can not be defined on the hybrid boudary */
                   
                    if( JEdge<HybrdBoundEdgeNum ){
                        printf("\nvsource can not be defined on\
 the hybrid boudary\n Sorry, exit!!");
                        exit(1);
                    }

                    flag=1;
                    break;

                }  /*  if(BoundEdgeStat[JEdge]==  */
            }     
            
            /* Can't find the Vsource edge */
            if( flag==0 ) 
            {
                printf("Sorry: Check your Vsource edge\n");
                exit(1);
            } /* end of if(flag) */