emap1.c

计算电磁学的３维变分标量有限元方法程序。EMAP1 employs a variational formulation。

                                          }
                                         }
                                        }
                     }

                 continue;
              }
              if(!strcmp(type, "ignore")) continue;
              fprintf(stderr, "Warning: Unrecognized attribute - %s\n", type);
           }              /*end of if*/ 
        }                 /* end of while */
       rewind(InF);
       /* Count the number of known and unknown vectors */
       TotNum_Known_Var=0;
       TotNum_Unknown_Var=0;
       for(nn=0;nn<TotNum_GBL_Nodes;nn++){
           for(i=0;i<=2;i++){
            if(Dat_Buffer[nn][i]==-999.) TotNum_Unknown_Var++;
            else                         TotNum_Known_Var++;
                            }                                                                                          
                                         } 
  }

/***************************************************************************
*
*   FUNCTION        : Fill_Vectors()
*   DESCRIPTION     : Defines forced and unforced vectors
*
****************************************************************************/

void Fill_Vectors()

  {
       int    i, nn;
 
       /* Define forced and unforced vectors */
       TotNum_Known_Var=0;
       TotNum_Unknown_Var=0;
       for(nn=0;nn<TotNum_GBL_Nodes;nn++){
           for(i=0;i<=2;i++){
                             if(Dat_Buffer[nn][i]==-999.){
                                                         UN_Vector[TotNum_Unknown_Var]=nn*3+i;
                                                         Loc_UN_Node[TotNum_Unknown_Var]=nn;
                                                         Loc_UN_NodeCord[TotNum_Unknown_Var]=i;
                                                         TotNum_Unknown_Var++;
                                                        }
                             else if(Dat_Buffer[nn][i]!=0.0){
                                                             FORC_Vector[TotNum_Known_Var]=nn*3+i;
                                                             FORC_Val[TotNum_Known_Var]=Dat_Buffer[nn][i];
                                                             TotNum_Known_Var++;
                                                            }                                                                                          
                            }
                                         } 
  }

/***************************************************************************
*
*   FUNCTION        : Count_Half_BandWidth()
*   DESCRIPTION     : Counts Bandwidth of the matrix A of Ax = B.
*
****************************************************************************/

 void Count_Half_BandWidth()
     {
      short int i,j,k,l,ColDiff,Val;
        {
        for(i=0;i<=TotNum_Unknown_Var-1;i++)
          {
          k=UN_Vector[i];
          for(j=i;j<=TotNum_Unknown_Var-1;j++)
            {
            l=UN_Vector[j];
            Val=Srch_NZ_Element_Column(k,l);
            if((Val>=0))
             {
              ColDiff=j-i;
              if(Half_BandWidth<=(ColDiff+1))
               {
                Half_BandWidth=(ColDiff+1);
               }
             }
           }
         }
         printf("Half_BandWidth = %d \n",Half_BandWidth);
      }
    }

/***************************************************************************
*
*   FUNCTION      : Partition_Global_Matrix()
*   DESCRIPTION   : Partitions Global matrix and applies boundary conditions.
*                   Construct A matrix of Ax = B matrix equation.
*
****************************************************************************/

void Partition_Global_Matrix()
     {
      short int i,j,k,l,ColDiff,g=0,val;
        {

        for(i=0;i<=TotNum_Unknown_Var-1;i++)  { k=UN_Vector[i];
        for(j=i;j<=TotNum_Unknown_Var-1;j++)  { l=UN_Vector[j];
        val=Srch_NZ_Element_Column(k,l);
        if((val>=0))
         {
          ColDiff=j-i;
          HalfBand_Matrix[i+1][ColDiff+1]+=GBL_Matrix_Data[k][val];
         }
                                              }

         g=0;
         for(j=0;j<=TotNum_Known_Var-1;j++)  { l=FORC_Vector[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        : Find_RHS_Vector()
*   DESCRIPTION     : Calculates RHS vectors B of Ax = B matrix equation.
*
****************************************************************************/

void Find_RHS_Vector()
    {
     int Count_i,Count_j,Srch_Ptr=0; double Temp_Buff,Multpl;
     {
      for(Count_i=0;Count_i<=TotNum_Unknown_Var-1;Count_i++)
        {
        for(Count_j=0;Count_j<=TotNum_Known_Var-1;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*FORC_Val[Count_j];
             RHS_Vector[Count_i]+=Multpl;
            }
          }
          RHS_Vector[Count_i]=RHS_Vector[Count_i]*-1.0;
        }
      }
    }

/***************************************************************************
*
*   FUNCTION        : Produce_Output()
*   DESCRIPTION     : Prints the efield data to the output files.
*
****************************************************************************/

void Produce_Output()
{
 int       i, j, k, x1, y1, z1, x2, y2, z2, Count_i, dflag=0;
 short int ColNos, RowNos;
 char      type[20], att1[20], att2[8], att3[8], att4[8];
 char      buffer[80];
 FILE      *OutF;

  for(Count_i=0;Count_i<=TotNum_Unknown_Var-1;Count_i++)   /*put data in Dat_Buffer*/
   {
    RowNos=Loc_UN_Node[Count_i]; ColNos=Loc_UN_NodeCord[Count_i];
    Dat_Buffer[RowNos][ColNos]=RHS_Vector[Count_i];
   }

        while (fgets(buffer, 80, InF))
        {
         if (buffer[0] == '#') {
                               if(OutF_0 != NULL) fprintf(OutF_0, "%s", buffer);
                               if(OutF_1 != NULL) fprintf(OutF_1, "%s", buffer);
                               if(OutF_2 != NULL) fprintf(OutF_2, "%s", buffer);
                               if(OutF_3 != NULL) fprintf(OutF_3, "%s", buffer);
                               continue;
                              }
        if(!strncmp(buffer, "default_output",14)) {dflag=1; continue;}
        if(!strncmp(buffer, "efield_output",13))
                   { 
                   sscanf(buffer, "%s%d%d%d%d%d%d%s%s%s%s", type, &x1, &y1, &z1, &x2, &y2, &z2, att1, att2, att3, att4);
                   if (x2<x1) swap(&x1, &x2);
                   if (y2<y1) swap(&y1, &y2);
                   if (z2<z1) swap(&z1, &z2);
                   if (!strcmp(att1, Out_FileName1)) OutF=OutF_1;
                   if (!strcmp(att1, Out_FileName2)) OutF=OutF_2;
                   if (!strcmp(att1, Out_FileName3)) OutF=OutF_3;
                   for(k=z1;k<=z2;k++){ 
                      for(j=y1;j<=y2;j++){
                         for(i=x1;i<=x2;i++){
                                             EfieldatNode(i, j, k, OutF); 
                                            }
                                         }
                                      } 
                   continue;
                   }

        }                 /* end of while */
        if(dflag == 1)    /* prints default output data */
              {

                for(Count_i=0;Count_i<=TotNum_GBL_Nodes-1;Count_i++)
                 {
                 fprintf(OutF_0," %6.3f\t%6.3f\t%6.3f\t",MeshNodeCord[Count_i][0],
                 MeshNodeCord[Count_i][1],MeshNodeCord[Count_i][2]);
                 fprintf(OutF_0," %f\t%f\t%f\t",
                 Dat_Buffer[Count_i][0],Dat_Buffer[Count_i][1],Dat_Buffer[Count_i][2]);
                 fprintf(OutF_0,"\n");
                 }
              }
        fclose(OutF_1);
        fclose(OutF_2);
        fclose(OutF_3);
        fclose(InF);
        fprintf(stdout,"\nProgram Successfully Completed: \n\n");
}

/***************************************************************************
*
*   FUNCTION        : Band_Matrix_Solver()
*   DESCRIPTION     : Solves for x variables of Ax = B matrix equation.
*
****************************************************************************/

void Band_Matrix_Solver()
  {
  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=TotNum_Unknown_Var;
 Column_Number=Half_BandWidth;
 for(Count_i=1;Count_i<=Row_Number;Count_i++)
  {
  for(Count_j=2;Count_j<=Column_Number;Count_j++)
   {
   if(HalfBand_Matrix[Count_i][Count_j]!=0.0)
   {
   Pivot_Emement=HalfBand_Matrix[Count_i][Count_j]/HalfBand_Matrix[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*HalfBand_Matrix[Count_i][Count_k];
    HalfBand_Matrix[Count_n][Count_l] -= Multpl_1;
    }
    Multpl_2=Pivot_Emement*RHS_Vector[Count_i-1];
    RHS_Vector[Count_n-1] -= Multpl_2;
   }
  }
  }


 RHS_Vector[Row_Number-1] /= HalfBand_Matrix[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(HalfBand_Matrix[Count_i][Count_j]!=0.0)
   {
   Count_k=Count_i+Count_j-1;
   Multpl_1=HalfBand_Matrix[Count_i][Count_j]*RHS_Vector[Count_k-1];
   RHS_Vector[Count_i-1] -= Multpl_1;
   }
  }
  RHS_Vector[Count_i-1] /= HalfBand_Matrix[Count_i][1];
  }
 }
 }

/***************************************************************************
*
*   FUNCTION    : Srch_NZ_Element_Column()
*   DESCRIPTION : Searches for column number of non-zero element
*                 within each row of Global matrix. 
*
****************************************************************************/

int Srch_NZ_Element_Column(RowNum,ColNum)
   { 
   int Count_i;
    {
     for(Count_i=0;Count_i<=GBL_Matrix_Index[RowNum]-1;Count_i++)
       {
        if(ColNum==(GBL_Matrix_ColNos[RowNum][Count_i]))
             {return Count_i;
             }
       }
     return -1;
    }
   }

/***************************************************************************
*
*   FUNCTION    : Locate_NZ_Element_Column()
*   DESCRIPTION : Searches for column number of non-zero element
*                 within each row of A matrix (?).
*
****************************************************************************/

int Locate_NZ_Element_Column(RowNum,ColNum)
     {
      int Count_i;
       {
        for(Count_i=0;Count_i<=FRCD_Matrix_MaxColNos-1;Count_i++)
          {
           if(ColNum==(FORC_Matrix_ColNos[RowNum][Count_i]))
             {return Count_i;
             }
          } 
        return -1;
       }
     }

/***************************************************************************
*
*   FUNCTION    : Determinant_TetraMatrix()
*   DESCRIPTION : Calculates determiment of a 4 x 4 matrix. 
*
****************************************************************************/

double Determinant_TetraMatrix()
{
int a,b,i,j,k,jj,kk,ll,m=0;
double Buff[4][4],Buff1[4],Buff2[4];
double TempDet=0.0,ppp=0.0,pp;
{

for(a=0;a<=3;a++)
 for(b=0;b<=3;b++)
  {Buff[a][b]=Volume_Cord[a][b];}

for(i=0;i<=3;i++) {ll=0;       
 for(j=0;j<=3;j++)  {
  if(j!=i) {

   for(k=0;k<=3;k++)  
    { if((k!=i)&&(k!=j)) {Buff1[m]=k; m++;} }

m=0; jj=Buff1[0];kk=Buff1[1];
pp=Buff[1][j]*(Buff[2][jj]*Buff[3][kk]-Buff[3][jj]*Buff[2][kk]);
ppp+=pp*pow(-1.0,(double)ll); pp=0.0;ll++;
           }
                    }
Buff2[i]=Buff[0][i]*ppp; ppp=0.0;
TempDet+=(Buff2[i])*pow(-1.0,(double)i);
                 }
TempDet=TempDet/6.0; return (TempDet*-1.0);
}
}

/***************************************************************************
*
*   FUNCTION    : swap(x,y)
*   DESCRIPTION : swap the values of x and y
*
****************************************************************************/

void swap(x,y)
        int *x,*y;
        {
        int temp;
        temp = *x; *x=*y; *y=temp;
        }

/***************************************************************************
*   FUNCTION        : EfieldatNode()
****************************************************************************/

void EfieldatNode(X,Y,Z,Out_F)

int X,Y,Z;
FILE *Out_F;
{
  int Count_i,Xcount,Ycount;
  int Xnode, Ynode, Znode;
  Xcount=(int)Dimension_X;Ycount=(int)Dimension_Y;
  Xnode=(X-Min_X);Ynode=(Y-Min_Y);Znode=(Z-Min_Z);
  Count_i=Xnode + ((Xcount+1)*Ynode) + ((Xcount+1)*(Ycount+1)*Znode);
  fprintf(Out_F," %6.3f\t%6.3f\t%6.3f\t",MeshNodeCord[Count_i][0],
  MeshNodeCord[Count_i][1],MeshNodeCord[Count_i][2]);
  fprintf(Out_F," %f\t%f\t%f\t",
  Dat_Buffer[Count_i][0],Dat_Buffer[Count_i][1],Dat_Buffer[Count_i][2]);
  fprintf(Out_F,"\n");

}
/*************************** END OF PROGRAM ********************************/