emap3.c

计算电磁学的３维矢量有限元程序 EMAP-3 is a vector (edge element) code. Vector codes are generally not affected

      }
    }

/***************************************************************************
*   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++){