emap3.c

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

                         nn=k*XdiM*YdiM+j*XdiM+i;
                         Epsilon[nn]=atof(att1);
                         Sigma[nn]=atof(att2);
                                            }
                                         }
                                      } 
                   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;
              /* left plate */
                           for(k=z1;k<z2;k++){
                             for(j=y1;j<y2;j++){
                                                nn=x2*y2*k+j*x2;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][8]]=1;
                                                EdgeStatus[xn[nn][8]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                               }
                                              }
              /* right plate */
                           for(k=z1;k<z2;k++){
                             for(j=y1;j<y2;j++){
                                                nn=x2*y2*k+(x2-1)+j*x2;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][9]]=1;
                                                EdgeStatus[xn[nn][9]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                               }
                                              }
              /* bottom plate */
                           for(k=z1;k<z2;k++){
                             for(i=x1;i<x2;i++){
                                                nn=y2*x2*k+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][6]]=1;
                                                EdgeStatus[xn[nn][6]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                               }
                                              }
              /* top plate */
                           for(k=z1;k<z2;k++){
                             for(i=x1;i<x2;i++){
                                                nn=y2*x2*k+x2*(y2-1)+i;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][11]]=1;
                                                EdgeStatus[xn[nn][11]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                               }
                                              }

              /* front plate */
                           for(j=y1;j<y2;j++){
                             for(i=x1;i<x2;i++){
                                                nn=j*x2+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][2]]=1;
                                                EdgeStatus[xn[nn][2]]=0.0;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                               }
                                              }
              /* back plate */
                           for(j=y1;j<y2;j++){
                             for(i=x1;i<x2;i++){
                                                nn=(z2-1)*y2*x2+j*x2+i;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                                ForceStat[xn[nn][15]]=1;
                                                EdgeStatus[xn[nn][15]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                               }
                                              }
             }
             if(!strcmp(type, "conductor"))
             {
                   /* fprintf(stderr, "conductor \n"); */
                   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*XdiM*YdiM+j*XdiM+i;
                                           Epsilon[nn]=0;
                                           Sigma[nn]=1.0E+09;
                                           for(tmp=1;tmp<19;tmp++){
                                                                  ForceStat[xn[nn][tmp]]=1;
                                                                  EdgeStatus[xn[nn][tmp]]=0.0;
                                                                  }
                                          }
                                       }
                                     }
               /* two-dimensional conductor */
                  if(x1==x2){     /*  y-z plane */
                           for(j=y1;j<y2;j++){
                             for(k=z1;k<z2;k++){
                                                if(x1!=XdiM)   {
                                                nn=k*XdiM*YdiM+j*XdiM+x1;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][8]]=1;
                                                EdgeStatus[xn[nn][8]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                                            }
                                                if(x1!=0)   {
                                                nn=k*XdiM*YdiM+j*XdiM+x1-1;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][9]]=1;
                                                EdgeStatus[xn[nn][9]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                                           }
                                               }
                                              }
                           }
                 if(y1==y2){     /*  x-z plane */
                           for(i=x1;i<x2;i++){
                             for(k=z1;k<z2;k++){
                                                if(y1!=YdiM)   {
                                                nn=k*XdiM*YdiM+y1*XdiM+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][5]]=1;
                                                EdgeStatus[xn[nn][5]]=0.0;
                                                ForceStat[xn[nn][6]]=1;
                                                EdgeStatus[xn[nn][6]]=0.0;
                                                ForceStat[xn[nn][7]]=1;
                                                EdgeStatus[xn[nn][7]]=0.0;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                                            }
                                                if(y1!=0)   {
                                                nn=k*XdiM*YdiM+(y1-1)*XdiM+i;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                                ForceStat[xn[nn][10]]=1;
                                                EdgeStatus[xn[nn][10]]=0.0;
                                                ForceStat[xn[nn][11]]=1;
                                                EdgeStatus[xn[nn][11]]=0.0;
                                                ForceStat[xn[nn][12]]=1;
                                                EdgeStatus[xn[nn][12]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                                           }
                                               }
                                              }
                           }
                 if(z1==z2){     /*   x-y plane */
                           for(i=x1;i<x2;i++){
                             for(j=y1;j<y2;j++){
                                                if(z1!=ZdiM)   {
                                                nn=z1*XdiM*YdiM+j*XdiM+i;
                                                ForceStat[xn[nn][0]]=1;
                                                EdgeStatus[xn[nn][0]]=0.0;
                                                ForceStat[xn[nn][1]]=1;
                                                EdgeStatus[xn[nn][1]]=0.0;
                                                ForceStat[xn[nn][2]]=1;
                                                EdgeStatus[xn[nn][2]]=0.0;
                                                ForceStat[xn[nn][3]]=1;
                                                EdgeStatus[xn[nn][3]]=0.0;
                                                ForceStat[xn[nn][4]]=1;
                                                EdgeStatus[xn[nn][4]]=0.0;
                                                            }
                                                if(z1!=0)   {
                                                nn=(z1-1)*XdiM*YdiM+j*XdiM+i;
                                                ForceStat[xn[nn][13]]=1;
                                                EdgeStatus[xn[nn][13]]=0.0;
                                                ForceStat[xn[nn][14]]=1;
                                                EdgeStatus[xn[nn][14]]=0.0;
                                                ForceStat[xn[nn][15]]=1;
                                                EdgeStatus[xn[nn][15]]=0.0;
                                                ForceStat[xn[nn][16]]=1;
                                                EdgeStatus[xn[nn][16]]=0.0;
                                                ForceStat[xn[nn][17]]=1;
                                                EdgeStatus[xn[nn][17]]=0.0;
                                                           }
                                               }
                                              }
                           }
 
               /* one-dimensional conductor */
 
                  if ((x1==x2)&&(y1==y2)){
                                          for(k=z1;k<z2;k++){
                                                            if((x1!=XdiM) & (y1!=YdiM)) {
                                                             nn=k*XdiM*YdiM+y1*XdiM+x1;
                                                             ForceStat[xn[nn][5]]=1;
                                                             EdgeStatus[xn[nn][5]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (y1!=YdiM)) {
                                                             nn=k*XdiM*YdiM+y1*XdiM+x1-1;
                                                             ForceStat[xn[nn][7]]=1;
                                                             EdgeStatus[xn[nn][7]]=0.0;
                                                                                  }
                                                            if((x1!=XdiM) & (y1!=0)) {
                                                             nn=k*XdiM*YdiM+(y1-1)*XdiM+x1;
                                                             ForceStat[xn[nn][10]]=1;
                                                             EdgeStatus[xn[nn][10]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (y1!=0)) {
                                                             nn=k*XdiM*YdiM+(y1-1)*XdiM+x1-1;
                                                             ForceStat[xn[nn][12]]=1;
                                                             EdgeStatus[xn[nn][12]]=0.0;
                                                                                  }
                                                            }
                                         }
                  if ((x1==x2)&&(z1==z2)){
                                          for(j=y1;j<y2;j++){
                                                            if((x1!=XdiM) & (z1!=ZdiM)) {
                                                             nn=z1*XdiM*YdiM+j*XdiM+x1;
                                                             ForceStat[xn[nn][1]]=1;
                                                             EdgeStatus[xn[nn][1]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (z1!=ZdiM)) {
                                                             nn=z1*XdiM*YdiM+j*XdiM+x1-1;
                                                             ForceStat[xn[nn][3]]=1;
                                                             EdgeStatus[xn[nn][3]]=0.0;
                                                                                  }
                                                            if((x1!=XdiM) & (z1!=0)) {
                                                             nn=(z1-1)*XdiM*YdiM+j*XdiM+x1;
                                                             ForceStat[xn[nn][14]]=1;
                                                             EdgeStatus[xn[nn][14]]=0.0;
                                                                                  }
                                                            if((x1!=0) & (z1!=0)) {
                                                             nn=(z1-1)*XdiM*YdiM+j*XdiM+x1-1;
                                                             ForceStat[xn[nn][16]]=1;
                                                             EdgeStatus[xn[nn][16]]=0.0;
                                                                                  }
                                                            }
                                         }
                  if ((y1==y2)&&(z1==z2)){
                                          for(i=x1;i<x2;i++){