emap1.c

一个求解三维时变场的有限元算法,C语言实现,根据MAILE提出的算法实现.

     char   type[20], att1[20], att2[18], att3[18], att4[18];
     double tmp1, tmp2, cell_dim=1.0;
     char   buffer[80];
 
     /* READ THE INPUT FILE:  FIRST PASS */
 
        while (fgets(buffer, 80, InF)) {
        if (buffer[0] == '#') {continue;}
        if (!strncmp(buffer, "celldim", 7)) { sscanf(buffer, "%s%lf%s", type,
                                              &tmp1, att1);
                                              if (!strcmp(att1, "cm")) tmp2=0.01;
                                              else if(!strcmp(att1, "m"))
                                                      tmp2=1.0;
                                              else if(!strcmp(att1, "mm"))
                                                      tmp2=0.001;
                                              else if(!strcmp(att1, "in"))
                                                      tmp2=2.54*0.01;
                                              else {fprintf(stderr, "unrecognized cell dimension.\n");
                                                    exit(1);}
                                              cell_dim=tmp1*tmp2;
                                              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;
                   }
        else 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"))
             {
/*
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   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: boundary must have finite thickness/n");
                                                     exit(1);
                                                    }
                   continue;
*/
fprintf(stderr, "\n\nERROR: boundary must be specified by the option 'box' \n"); exit(1);
             }
             if(!strcmp(type, "dielectric"))
             {
                   /* fprintf(stderr, "dielectric \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   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: dielectric must have finite thickness/n");
                                                     exit(1);
                                                    }
                   continue;
             }
             if(!strcmp(type, "box"))
             {
                   /* fprintf(stderr, "box \n"); */
                   if((x1==x2)||(y1==y2)||(z1==z2)){ fprintf(stderr, "Error: invalid box:  ");
                                                   fprintf(stderr, "box must be three dimensional\n");
                                                   exit(1);}
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   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"))
             {
                   /* fprintf(stderr, "conductor \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   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"))
             {
                   /* fprintf(stderr, "aperture \n"); */
                if((x1!=x2) & (y1!=y2) & (z1!=z2)) {fprintf(stderr, "ERROR: apertures must be two dimensional\n");
                                                       exit(1);
                                                      }
                if(((x1==x2)&&(y1==y2))||((x1==x2)&&(z1==z2))||((z1==z2)&&(y1==y2)))
                                                   {fprintf(stderr, "ERROR: apertures must be two dimensional\n");
                                                       exit(1);
                                                      }
                   continue;
             }
             if(!strcmp(type, "seam"))
             {
                   /* fprintf(stderr, "seam \n"); */
                   continue;
             }
             if(!strcmp(type,"esource"))
             {
                   /* fprintf(stderr, "esource \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   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"))
              {
                   /* fprintf(stderr, "msource \n"); */
                   if (x1>x2) swap(&x1,&x2); if(y1>y2) swap(&y1,&y2); if(z1>z2) swap(&z1,&z2);
                   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, "Unrecognized attribute %s\n\n", type);
            }                 /*end of else*/
        }                     /* end of while */
 
 
        if(freq ==0) fprintf(stderr, "\nWARNING: No sources were specified, 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");
        OperateFreq = freq;  CellDimension = cell_dim;
        Dimension_X=Max_X-Min_X; Dimension_Y=Max_Y-Min_Y; Dimension_Z=Max_Z-Min_Z;
        TotNum_HexaHedron=Dimension_X*Dimension_Y*Dimension_Z;
        TotNum_GBL_Nodes=(Dimension_X+1)*(Dimension_Y+1)*(Dimension_Z+1);
        rewind(InF);
 }

/***************************************************************************
*
*   FUNCTION        : Read_Input_Pass_2()
*   DESCRIPTION     : Reads information for each node from input datafile.
*                     Reads information about boundary conditions.
*
****************************************************************************/

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

  /*  Initialize dynamically allocated variables  */
    for(nn=0;nn<TotNum_GBL_Nodes;nn++)
       {for(i=0;i<=2;i++){Dat_Buffer[nn][i]=-999.;REL_Permittivity[nn]=1.0;}}

     I=Dimension_X+1; J=Dimension_Y+1;

   /*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, "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;

                   for(j=y1;j<=y2;j++){
                     for(k=z1;k<=z2;k++){
                                         nn=k*J*I+j*I+x1;
                                         Dat_Buffer[nn][1]=0.0;
                                         Dat_Buffer[nn][2]=0.0;
                                         nn=k*J*I+j*I+x2;
                                         Dat_Buffer[nn][1]=0.0;
                                         Dat_Buffer[nn][2]=0.0;
                                        }
                                      }                           
                   for(i=x1;i<=x2;i++){
                     for(k=z1;k<=z2;k++){
                                         nn=k*J*I+y1*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][2]=0.0;
                                         nn=k*J*I+y2*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][2]=0.0;                              
                                        }
                                      }                           
                   for(i=x1;i<=x2;i++){
                     for(j=y1;j<=y2;j++){
                                         nn=z1*J*I+j*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][1]=0.0;
                                         nn=z2*J*I+j*I+i;
                                         Dat_Buffer[nn][0]=0.0;
                                         Dat_Buffer[nn][1]=0.0;
                                        }
                                      }                           
                  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(k=z1;k<=z2;k++){
                      for(j=y1;j<=y2;j++){
                         for(i=x1;i<=x2;i++){
                         nn=k*I*J+j*I+i;
                         REL_Permittivity[nn]=atof(att1);
                                            }
                                         }
                                      } 
                   continue;
             }
             if(!strcmp(type, "conductor"))
             {
                   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*J*I+j*I+i;
                                           for(tmp=0;tmp<2;tmp++){
                                                                  Dat_Buffer[nn][tmp]=0.0;
                                                                 }
                                          }
                                       }
                                     }
               /* two-dimensional conductor */
                  if(x1==x2){
                           for(j=y1;j<=y2;j++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+j*I+x1;
                                                Dat_Buffer[nn][1]=0.0;
                                                Dat_Buffer[nn][2]=0.0;
                                               }
                                              }                           
                           }
                 if(y1==y2){
                           for(i=x1;i<=x2;i++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+y1*I+i;
                                                Dat_Buffer[nn][0]=0.0;
                                                Dat_Buffer[nn][2]=0.0;                                                
                                               }
                                              }                           
                           }
                 if(z1==z2){
                           for(i=x1;i<=x2;i++){
                             for(j=y1;j<=y2;j++){
                                                nn=z1*J*I+j*I+i;
                                                Dat_Buffer[nn][0]=0.0;
                                                Dat_Buffer[nn][1]=0.0;
                                               }
                                              }                           
                           }
               /* one-dimensional conductor */
                  if ((x1==x2)&&(y1==y2)){
                                          for(k=z1;k<z2;k++){
                                                             nn=k*J*I+y1*I+x1;
                                                             Dat_Buffer[nn][2]=0.0;
                                                            } 
                                         } 
                  if ((x1==x2)&&(z1==z2)){
                                          for(j=y1;j<y2;j++){
                                                             nn=z1*J*I+j*I+x1;
                                                             Dat_Buffer[nn][1]=0.0;
                                                            } 
                                         }
                  if ((y1==y2)&&(z1==z2)){
                                          for(i=x1;i<x2;i++){
                                                             nn=z1*J*I+y1*I+i;
                                                             Dat_Buffer[nn][0]=0.0;
                                                            } 
                                         }                                           
                   continue;
             }
             if(!strcmp(type, "aperture"))
             {
                   /* fprintf(stderr, "aperture \n"); */
                  if(x1==x2){
                           for(j=y1;j<=y2;j++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+j*I+x1;
                                                Dat_Buffer[nn][0]=-999.;
                                                Dat_Buffer[nn][1]=-999.;
                                                Dat_Buffer[nn][2]=-999.;
                                               }
                                              }                           
                           }
                 if(y1==y2){
                           for(i=x1;i<=x2;i++){
                             for(k=z1;k<=z2;k++){
                                                nn=k*J*I+y1*I+i;
                                                Dat_Buffer[nn][0]=-999.;
                                                Dat_Buffer[nn][1]=-999.;
                                                Dat_Buffer[nn][2]=-999.;                                                
                                               }
                                              }                           
                           }
                 if(z1==z2){
                           for(i=x1;i<=x2;i++){
                             for(j=y1;j<=y2;j++){
                                                nn=z1*J*I+j*I+i;
                                                Dat_Buffer[nn][0]=-999.;
                                                Dat_Buffer[nn][1]=-999.;
                                                Dat_Buffer[nn][2]=-999.;
                                               }
                                              }                           
                           }

                   continue;
             }
             if(!strcmp(type, "seam"))
             {
                   /* fprintf(stderr, "seam \n"); */
                   /*  need to define seam here  */
                   continue;
             }
             if(!strcmp(type, "esource"))
             {
                   /* fprintf(stderr, "esource \n"); */
                   if(freq != 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*I*J+j*I+i;
                                            Dat_Buffer[nn][1]=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*I*J+j*I+i;
                                           Dat_Buffer[nn][2]=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*I*J+j*I+i;
                                           Dat_Buffer[nn][0]=atof(att3);