sift5.c

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

    int Count_k; double Buff[3],Buff1[3],Buff2[3],Buff3[3];
 
    for(Count_k=0;Count_k<=2;Count_k++)
    {
        Buff1[Count_k]=Cord[BFaceNode[TrngleNum][0]-1][Count_k];
        Buff2[Count_k]=Cord[BFaceNode[TrngleNum][1]-1][Count_k];
        Buff3[Count_k]=Cord[BFaceNode[TrngleNum][2]-1][Count_k];

    }  /* for(Count_k  */

    VTXadd2(Buff1,Buff2,Buff3,Buff);
    CenTroid[TrngleNum][0] = (1.0/3.0) * Buff[0];
    CenTroid[TrngleNum][1] = (1.0/3.0) * Buff[1];
    CenTroid[TrngleNum][2] = (1.0/3.0) * Buff[2];

}   /* end of ComputeTrngleCentroid() */
 
 
 


/*************************************************************************
Prototype:    void    Read_SIF_File() 
Description:    Reads input file 
Input value:     none 
Return value:     none 
Global value used: CellDim, XdiM, YdiM, ZdiM, XDielmat, YDielmat, 
                   ZDielmat,  Xintcond, Yintcond, Zintcond, BndcondNum,  
                   Xextcond, Yextcond, Zextcond, XextVsrc, YextVsrc, 
                   ZextVsrc, XintEsrc, YintEsrc, ZintEsrc, ExtcondNum, 
                   BndcondNum, ExtVsrcNum, freq, SourceType, PolInc, 
                   ThetaInc, PhiInc, OperateFreq, CellDimension , Min_X, 
                   Min_Y, Min_Z, Epsilon, Sigma, DielmatNum, IntcondNum, 
                   VsourceMag, IntEsrcNum. 
Global value modified: Min_X, Min_Y, Min_Z, XdiM, YdiM, ZdiM, CellDim,  
                       DielmatNum, IntcondNum, BndcondNum, ExtcondNum,  
                       ExtVsrcNum, VsourceMag, IntEsrcNum, ThetaInc,
                       PhiInc, OperateFreq, CellDimension. 
Subroutines called: none 
*************************************************************************/
void ReadSIFFile()
{
    int    x1, y1, z1, x2, y2, z2, delx, dely, delz, i;
    char   axis, type[20], buffer[80];
    double  tmp1, tmp2, freq = 0, CellDim = 1;

    Min_X = 10000; Min_Y=10000; Min_Z=10000;
    Max_X = 0;     Max_Y=0;     Max_Z=0;
    CondNum = 0;
 
    printf("\nProcessing Key Words:\n\n");
 
    while ( fscanf(InF,"%s", type)!=EOF )
    {
        /* # line means comments */
        if( type[0] == '#' ) 
        {
            /* get rid of those comment lines */
            fgets(buffer, 80, InF);
            continue;
        }

        if( !strncmp(type, "unit", 4) ) 
        {
             printf("\tunit\n");
 
            if( fscanf(InF, "%lf %s", &tmp1, buffer)!=2 )
            {
                fprintf(stderr, "unit: missing some parameters\n");
                exit(1);
            }

            if( !strcmp(buffer, "cm") ) tmp2=0.01;
            else if(!strcmp(buffer, "m"))  tmp2=1;
            else if(!strcmp(buffer, "mm")) tmp2=0.001;
            else 
            {
                fprintf(stderr,"ERROR: unrecognized cell dimension.\n");
                exit(1);
            }
            CellDim=tmp1*tmp2;  /* get the unit in meter */
            continue; 
        }   

        if( !strcmp(type,"boundary") ) 
        {
            printf("\tboundary\n");
            if( fscanf(InF, "%d %d %d %d %d %d", &x1, &y1, &z1,
                &x2, &y2, &z2)!=6 )
            {
                fprintf(stderr,"Boundary: missing some parameters\n");
                exit(1);
            }      
  
            XBoundDim[0] = x1; 
            XBoundDim[1] = x2;
            YBoundDim[0] = y1; 
            YBoundDim[1] = y2;
            ZBoundDim[0] = z1; 
            ZBoundDim[1] = z2;
        
	    /* to check the boundary parameters */
            if( (x1==x2)||(y1==y2)||(z1==z2) ) 
            {
                fprintf(stderr, "boundary: dimension is invalid.\n");
                exit(1);
            }
            continue;
        
        } /* end of "boundary" */
        
        if( !strcmp(type,"conductor") )
        {
            printf("\tconductor\n");
            if( fscanf(InF, "%d %d %d %d %d %d %d %d %d", &x1, &y1, &z1,
                          &x2, &y2, &z2, &delx, &dely, &delz)!=9 )
            {
                fprintf(stderr, "conductor: missing some parameters\n");
                exit(1);
            }

            /* check the parameters */
            if( (x1!=x2) && (y1!=y2) &&(z1!=z2) )
            {
           	fprintf(stderr, "Conductor should be two demension and \
parallel one axis!\n");
           	exit(1);
            }
           
            conductor[CondNum].x1=x1;
            conductor[CondNum].y1=y1;
            conductor[CondNum].z1=z1;
            conductor[CondNum].x2=x2;
            conductor[CondNum].y2=y2;
            conductor[CondNum].z2=z2;
            conductor[CondNum].DelX= delx;
            conductor[CondNum].DelY= dely;
            conductor[CondNum].DelZ= delz;
            
            CondNum++;
            continue;
        }   /* end of "conductor */    
              
       
        if( !strcmp(type, "dielectric") ) 
        {
            printf("\tdielectric\n");
            if( fscanf(InF, "%d %d %d %d %d %d %lf %lf", &x1, &y1, &z1, 
                    &x2, &y2, &z2, &tmp1, &tmp2)!=8 )
            {
                fprintf(stderr,"dielectric: missing some parameters\n");
                exit(1);
            }
            if( (x1==x2)||(y1==y2)||(z1==z2) )
            {
                fprintf(stderr,"ERROR:  diel thickness is not finite.\n");
                exit(1);
            }
            XDielmat[0][DielmatNum]=x1;
            XDielmat[1][DielmatNum]=x2;
            YDielmat[0][DielmatNum]=y1;
            YDielmat[1][DielmatNum]=y2;
            ZDielmat[0][DielmatNum]=z1;
            ZDielmat[1][DielmatNum]=z2;
            Epsilon[DielmatNum]=tmp1; 
            Sigma[DielmatNum]=tmp2;
            DielmatNum++;
            continue;
        }  /* end "dielectric" */

        if( !strcmp(type,"vsource") )
        {
            printf("\tvsource\n");
            SourceType='V';
            if( fscanf(InF, "%d %d %d %d %d %d %lf %c %lf",
                &x1, &y1, &z1, &x2, &y2, &z2, &freq, &axis, &tmp1)!=9 )
            {
                fprintf(stderr, "vsource: missing some parameters\n");
                exit(1);
            }

            if( (x1==x2)&&(y1==y2)&&(z1==z2) ) 
            {
               fprintf(stderr, "ERROR: vsource cannot be defined at \
a point.\n"); 
               exit(1);
            }
            if( (axis != 'x')  && (axis!= 'y') && (axis != 'z') )
            {
                fprintf(stderr, "vsource: unrecognized polarization\n");
                exit(1);
            }
 
            XextVsrc[0][ExtVsrcNum] = x1;
            XextVsrc[1][ExtVsrcNum] = x2;
            YextVsrc[0][ExtVsrcNum] = y1;
            YextVsrc[1][ExtVsrcNum] = y2;
            ZextVsrc[0][ExtVsrcNum] = z1;
            ZextVsrc[1][ExtVsrcNum] = z2;
            ExtVsrcMag[ExtVsrcNum] = tmp1;
            ExtVsrcNum++;
            
            continue;
        }   /* end "vsource" */
   
        if( !strcmp(type,"isource") ) 
        {
            SourceType='I';
            printf("\tisource\n");
            if( fscanf(InF, "%d %d %d %d %d %d  %lf %c %lf",
                &x1, &y1, &z1, &x2, &y2, &z2, &freq, &axis,&tmp1)!=9)
            {
                fprintf(stderr, "isource: missing some parameters\n");
                exit(1);
            }

            XintEsrc[0][IntEsrcNum] = x1;
            XintEsrc[1][IntEsrcNum] = x2;
            YintEsrc[0][IntEsrcNum] = y1;
            YintEsrc[1][IntEsrcNum] = y2;
            ZintEsrc[0][IntEsrcNum] = z1; 
            ZintEsrc[1][IntEsrcNum] = z2;
            IntEsrcAxis[IntEsrcNum]=axis;
            IntEsrceMag[IntEsrcNum]=tmp1;        
            IntEsrcNum++;
        
            if( (x1==x2)&&(y1==y2)&&(z1==z2) ) 
            {
        	fprintf(stderr, "ERROR: isource cannot be defined \
at a point.\n"); 
                exit(1);
            }

            if( (axis != 'x') && (axis!= 'y') && ( axis != 'z') )
            {
        	fprintf(stderr, "isource: unrecognized polarization\n");
                exit(1); 
            }
            continue;

        }  /* end "esource" */

        if( !strcmp(type,"eplane") )
        {
            printf("\teplane\n");
            if( fscanf(InF, "%lf %lf %lf  %lf %lf %lf", &freq, &E_Theta, \
&E_Phi,  &K_Theta, &K_Phi, &PlaneWaveE_Mag)!=6 )
            {
                fprintf(stderr,"%lf  %lf %lf  %lf %lf %lf", freq, E_Theta, \
E_Phi, K_Theta, K_Phi, PlaneWaveE_Mag );
                fprintf(stderr, "eplane: missing some parameters\n");
                exit(1);
            }

            SourceType='P';
            continue;
        } /* end "eplane" */
    
    
        if( !strcmp(type,"celldim") )
        {
            char axis;
            double start, end, step;
      
            printf("\tcelldim\n");
      
            if( fscanf(InF, "\t %lf  %lf %lf %c", &start, &end, &step, \
&axis)!=4 )
            {
                fprintf(stderr, "celldim: missing some parameters\n");
                exit(1);
            }
      
            if( axis=='x' )
            {
               CellUnitX[CellUnitXNum].start=start;
               CellUnitX[CellUnitXNum].end=end;
               CellUnitX[CellUnitXNum].step=step;
               CellUnitXNum++;
            }
            else if(axis=='y')
            {
                CellUnitY[CellUnitYNum].start=start;
                CellUnitY[CellUnitYNum].end=end;
                CellUnitY[CellUnitYNum].step=step;
                CellUnitYNum++;
            }
      
            else if( axis=='z' )
            {
                CellUnitZ[CellUnitZNum].start=start;
                CellUnitZ[CellUnitZNum].end=end;
                CellUnitZ[CellUnitZNum].step=step;
                CellUnitZNum++;
            }
        
            else 
            {
                fprintf(stderr,"celldim: unrecognized axix\n");
                exit(1);
            }
            continue;
      
        }  

        if( !strcmp(type, "resistor") ) 
        {
            double tmp;

            printf("\tresistor\n");
            fscanf(InF,"%d %d %d %d %d %d %lf", &x1, &y1, &z1, &x2, &y2, \
&z2, &tmp);
            ResisterLocX[ResisterNum][0]=x1;
            ResisterLocX[ResisterNum][1]=x1;
            ResisterLocY[ResisterNum][0]=y1;
            ResisterLocY[ResisterNum][1]=y2;
            ResisterLocZ[ResisterNum][0]=z1;
            ResisterLocZ[ResisterNum][1]=z2;
            ResisterValue[ResisterNum]=tmp;
            ResisterNum++;
            continue;
        }

        if( !strcmp(type,"output") ) 
        {
   	    char filename[20];
   	    char axis;
   	
   	    printf("\toutput\n");
            if( fscanf(InF, "%d %d %d %d %d %d  %c %s", &x1, &y1, &z1, 
                      &x2, &y2, &z2, &axis, filename)!=8 )
            {
                fprintf(stderr, "output: misisng some parameters\n");
                exit(1);
            }
        
            if( axis!='x' && axis!='y'&& axis!='z' ) 
            {
                fprintf(stderr, "output: mising axis parameter\n");    
                exit(1);
            } 
            output_x[0][output_num]=x1;            
            output_x[1][output_num]=x2;            
            output_y[0][output_num]=y1;            
            output_y[1][output_num]=y2;            
            output_z[0][output_num]=z1;            
            output_z[1][output_num]=z2;  
            output_axis[output_num]=axis;
            strcpy(output_file[output_num], filename);
            output_num++;
            continue;
        }    
   
        if( !strcmp(type, "default_out") ) 
        {
            char filename[20];

            printf("\tdefault_out\n");       
            flag_default_out = 1;    /* set the flag to true */
            if( fscanf(InF, "%s", filename)!=1 ) 
            {
                fprintf(stderr, "default_out: mising output file name\n");
                exit(1);
            }
            strcpy(default_out_file, filename);
            continue;
        }
       
  
        fprintf(stderr,"error: Unrecognized keyword %s\n", type);
        exit(1);
    } /* end of sscanf(***, type ) */
 
    if( freq==0 ) 
       fprintf(stderr, "WARNING: input file does not have a source!\n\n");


    /* detect conductative boundary and inner conductors */ 
    for(i=0; i<CondNum; i++) 
    {
        x1=conductor[i].x1;
        y1=conductor[i].y1;
        z1=conductor[i].z1;
        x2=conductor[i].x2;
        y2=conductor[i].y2;
        z2=conductor[i].z2;
   
        /* conductative boundary maybe exist */
        if( ((x1==x2) && ((x1==XBoundDim[0]) || (x1==XBoundDim[1])))  ||
            ((y1==y2) && ((y1==YBoundDim[0]) ||