emap1.c

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

      tetrnod[2][0]=HexHed[cn][1];   tetlink[2][0]=1;
      tetrnod[2][1]=HexHed[cn][0];   tetlink[2][1]=0;
      tetrnod[2][2]=HexHed[cn][5];   tetlink[2][2]=5;
      tetrnod[2][3]=HexHed[cn][3];   tetlink[2][3]=3;
 
      tetrnod[3][0]=HexHed[cn][4];   tetlink[3][0]=4;
      tetrnod[3][1]=HexHed[cn][5];   tetlink[3][1]=5;
      tetrnod[3][2]=HexHed[cn][0];   tetlink[3][2]=0;
      tetrnod[3][3]=HexHed[cn][6];   tetlink[3][3]=6;
 
      tetrnod[4][0]=HexHed[cn][3];   tetlink[4][0]=3;
      tetrnod[4][1]=HexHed[cn][0];   tetlink[4][1]=0;
      tetrnod[4][2]=HexHed[cn][5];   tetlink[4][2]=5;
      tetrnod[4][3]=HexHed[cn][6];   tetlink[4][3]=6;

   }

/***************************************************************************
*
*   FUNCTION        : Find_Volume_Determinant()
*   DESCRIPTION     : Find the Volume determinant of a Tetrahedron.
*
****************************************************************************/

void Find_Volume_Determinant()
     {
      int Count_i,Count_l,Count_j,Count_k,lu[4];
       {
       for(Count_l=0;Count_l<=3;Count_l++)
         {
         lu[Count_l]=tetrnod[t][Count_l];  
         }
       for(Count_i=0;Count_i<=3;Count_i++)
         for(Count_j=0;Count_j<=3;Count_j++)
           {
           if(Count_j==3)
            {
            Volume_Cord[Count_i][Count_j]=1.0;
            VVolume_Cord[Count_i][Count_j]=1.0;
            }
           else
            {
            Count_k=lu[Count_i];
            Volume_Cord[Count_i][Count_j]=MeshNodeCord[Count_k-1][Count_j];
            VVolume_Cord[Count_i][Count_j]=MeshNodeCord[Count_k-1][Count_j];
            } 
           }
         {
         det=Determinant_TetraMatrix();
         }
      }
    }

/***************************************************************************
*
*   FUNCTION        : Find_CoFactor() 
*   DESCRIPTION     : Find the  Cofactors of a Volume determinant.
*
****************************************************************************/

void Find_CoFactor()
   {
   short int i,j,row,col,ro,co,counter=0;
   double cof[4][4];
     {

     for(row=0;row<=3;row++)
     for(col=0;col<=3;col++) {
     co=0;ro=0;counter=0;  

        for(i=0;i<=3;i++)
        for(j=0;j<=3;j++) {
        if((i!=row) && (j!=col)) {
        counter++;
        if(counter<4) ro=0;
        if((counter>3)&&(counter<7)) ro=1;
        else if(counter>6) ro=2;
        cof[ro][co]=VVolume_Cord[i][j];
        co++;co=co%3;            }
                          }
        ro=0;co=0;counter=0;
        cofs[row][col]=cof[0][0]*((cof[1][1]*cof[2][2])-(cof[1][2]*cof[2][1]))
                      -cof[0][1]*((cof[1][0]*cof[2][2])-(cof[2][0]*cof[1][2]))
                      +cof[0][2]*((cof[1][0]*cof[2][1])-(cof[2][0]*cof[1][1]));
        
        if((row+col)%2!=0)
        cofs[row][col]*=-1.0;
         }
      }
    }

/**************************************************************************
*
*   FUNCTION        : TetraHedron_SubMatrix()
*   DESCRIPTION     : Form 12 x 12 Submatrix for a Tetrahedron.
*
**************************************************************************/

void TetraHedron_SubMatrix()
{
 short int mm0,mm1,mm2,mm3,cnt1,cnt2,p,q,m,n,i,k;
 double real_perm,real_buff,delta;
 { 

   mm0=tetrnod[t][0];
   mm1=tetrnod[t][1];
   mm2=tetrnod[t][2];
   mm3=tetrnod[t][3];
   if((REL_Permittivity[mm0]==1.0) && (REL_Permittivity[mm1]==1.0) &&
      (REL_Permittivity[mm2]==1.0) && (REL_Permittivity[mm3]==1.0))
    {
        for(m=1;m<=4;m++)
        for(n=1;n<=4;n++)
        for(i=1;i<=3;i++)                          {
        p=3*(m-1)+i; q=3*(n-1)+i;
        if(m==n) delta=1.0; else delta=0.0;
        real_buff=((1+delta)*det)/20.0;
        real_buff=real_buff*(WaveNumber*WaveNumber);
        Tet_Matrix_Element[t][p-1][q-1]=real_buff; }
    }
   else
    {
        for(m=1;m<=4;m++)
        for(n=1;n<=4;n++)
          {
           cnt1=tetrnod[t][m-1]; cnt2=tetrnod[t][n-1];
           if(REL_Permittivity[cnt1]==REL_Permittivity[cnt2])  
             {
                real_perm=REL_Permittivity[cnt1];
                for(i=1;i<=3;i++)                          {
                p=3*(m-1)+i; q=3*(n-1)+i;
                if(m==n) delta=1.0; else delta=0.0;
                real_buff=((1+delta)*det)/20.0;
                real_buff=real_buff*(WaveNumber*WaveNumber);
                if(REL_Permittivity[cnt1]==1.0)
                     { if(i<=2) {real_buff=real_buff*real_perm;}}
                else {           real_buff=real_buff*real_perm; }
                Tet_Matrix_Element[t][p-1][q-1]=real_buff; }
             }
           else
             {
                if(REL_Permittivity[cnt1]>=REL_Permittivity[cnt2]) 
                  {real_perm=REL_Permittivity[cnt1];}
                else
                  {real_perm=REL_Permittivity[cnt2];}

                for(i=1;i<=3;i++)                           {
                p=3*(m-1)+i; q=3*(n-1)+i;
                if(m==n) delta=1.0; else delta=0.0;
                real_buff=((1+delta)*det)/20.0;
                real_buff=real_buff*(WaveNumber*WaveNumber);
                if(i<=2) {real_buff=real_buff*real_perm;}
                Tet_Matrix_Element[t][p-1][q-1]=real_buff;  }
              }

          } 

     }
  for(m=1;m<=4;m++)          
    for(n=1;n<=4;n++)
      for(i=1;i<=3;i++)                          
        {
        p=3*(m-1)+i; q=3*(n-1)+i;
        real_buff=(cofs[m-1][0]*cofs[n-1][0])+
                  (cofs[m-1][1]*cofs[n-1][1])+
                  (cofs[m-1][2]*cofs[n-1][2]);
        real_buff=(real_buff)/(36.0*det);
        Tet_Matrix_Element[t][p-1][q-1]-=real_buff;
        }

  for(m=1;m<=4;m++)
    for(n=1;n<=4;n++)
      for(i=1;i<=3;i++)
        for(k=1;k<=3;k++)
        {
          p=3*(m-1)+i; q=3*(n-1)+k;
          real_buff=(cofs[m-1][k-1]*cofs[n-1][i-1]);
          real_buff=(real_buff)/(36.0*det);
          Tet_Matrix_Element[t][p-1][q-1]+=real_buff;
        }

 for(p=0;p<=11;p++) for(q=0;q<=11;q++)
 if(fabs(Tet_Matrix_Element[t][p][q]) < 1.0e-10) 
 {Tet_Matrix_Element[t][p][q] = 0.0;}

 }
}

/***************************************************************************
*
*   FUNCTION        : HexaHedron_Matrix()
*   DESCRIPTION     : Form 24 x 24 Submatrix for a Hexahedron.
*
****************************************************************************/

void HexaHedron_Matrix()
   {
    short int i,j,t,mm,nn,m,n,k,l;
    {

   for(i=0;i<=23;i++) for(j=0;j<=23;j++)
   { Hex_Matrix_Element[i][j]=0.0;}
 
   for(t=0;t<=4;t++) 
    for(m=0;m<=3;m++) 
      for(n=0;n<=3;n++)
        {
          k=tetlink[t][m]*3; 
          l=tetlink[t][n]*3;
          mm=m*3; nn=n*3;
          for(i=0;i<=2;i++) 
            for(j=0;j<=2;j++)
              {
               if(Tet_Matrix_Element[t][mm+i][nn+j]!=0.0)
               {
               Hex_Matrix_Element[k+i][l+j]+=
               Tet_Matrix_Element[t][mm+i][nn+j];
               }
              }
        }
  for(m=0;m<=23;m++) for(n=0;n<=23;n++)
  if(fabs(Hex_Matrix_Element[m][n]) < 1.0e-10)
   {Hex_Matrix_Element[m][n] = 0.0;}
    }
   }

/***************************************************************************
*
*   FUNCTION        : HexaHedron_link_1()
*   DESCRIPTION     : Provide links between the nodes. ( Method 1)
*
****************************************************************************/

void HexaHedron_link_1()
    {
     connect1[0][0] = 1;  connect1[0][1] = 2;
     connect1[0][2] = 4;  connect1[0][3] =-1;
     connect1[0][4] =-1;  connect1[0][5] =-1;

     connect1[1][0] = 0;  connect1[1][1] = 2;
     connect1[1][2] = 3;  connect1[1][3] = 4;
     connect1[1][4] = 5;  connect1[1][5] = 7;

     connect1[2][0] = 0;  connect1[2][1] = 1;
     connect1[2][2] = 3;  connect1[2][3] = 4;
     connect1[2][4] = 6;  connect1[2][5] = 7;

     connect1[3][0] = 1;  connect1[3][1] = 2;
     connect1[3][2] = 7;  connect1[3][3] =-1;
     connect1[3][4] =-1;  connect1[3][5] =-1;

     connect1[4][0] = 0;  connect1[4][1] = 1;
     connect1[4][2] = 2;  connect1[4][3] = 5;
     connect1[4][4] = 6;  connect1[4][5] = 7;

     connect1[5][0] = 1;  connect1[5][1] = 4;
     connect1[5][2] = 7;  connect1[5][3] =-1;
     connect1[5][4] =-1;  connect1[5][5] =-1;

     connect1[6][0] = 2;  connect1[6][1] = 4;
     connect1[6][2] = 7;  connect1[6][3] =-1;
     connect1[6][4] =-1;  connect1[6][5] =-1;

     connect1[7][0] = 1;  connect1[7][1] = 2;
     connect1[7][2] = 3;  connect1[7][3] = 4;
     connect1[7][4] = 5;  connect1[7][5] = 6;
     }

/***************************************************************************
*
*   FUNCTION        : HexaHedron_link_2()
*   DESCRIPTION     : Provide links between the nodes. ( Method 2)
*
****************************************************************************/
 
void HexaHedron_link_2()
    {
     connect1[0][0] = 1;  connect1[0][1] = 2;
     connect1[0][2] = 3;  connect1[0][3] = 4;
     connect1[0][4] = 5;  connect1[0][5] = 6;
 
     connect1[1][0] = 0;  connect1[1][1] = 3;
     connect1[1][2] = 5;  connect1[1][3] =-1;
     connect1[1][4] =-1;  connect1[1][5] =-1;
 
     connect1[2][0] = 0;  connect1[2][1] = 3;
     connect1[2][2] = 6;  connect1[2][3] =-1;
     connect1[2][4] =-1;  connect1[2][5] =-1;
 
     connect1[3][0] = 0;  connect1[3][1] = 1;
     connect1[3][2] = 2;  connect1[3][3] = 5;
     connect1[3][4] = 6;  connect1[3][5] = 7;
 
     connect1[4][0] = 0;  connect1[4][1] = 5;
     connect1[4][2] = 6;  connect1[4][3] =-1;
     connect1[4][4] =-1;  connect1[4][5] =-1;
 
     connect1[5][0] = 0;  connect1[5][1] = 1;
     connect1[5][2] = 3;  connect1[5][3] = 4;
     connect1[5][4] = 6;  connect1[5][5] = 7;
 
     connect1[6][0] = 0;  connect1[6][1] = 2;
     connect1[6][2] = 3;  connect1[6][3] = 4;
     connect1[6][4] = 5;  connect1[6][5] = 7;
 
     connect1[7][0] = 3;  connect1[7][1] = 5;
     connect1[7][2] = 6;  connect1[7][3] =-1;
     connect1[7][4] =-1;  connect1[7][5] =-1;
     }

/***************************************************************************
*
*   FUNCTION        : Find_Global_Matrix()
*   DESCRIPTION     : Construct 3N x 3N Global matrix. N =Total # of nodes.
*
****************************************************************************/

void Find_Global_Matrix()
{
 int l,s,x,y,z,i,j,ma,qa,k=0,c;
 int tem=0,rval,fval,bck;
{

for(k=0;k<=7;k++)
{

ma=HexHed[cn][k]; ma--;s=ma*3;c=k*3; 

for(i=0;i<=2;i++)
 for(j=0;j<=2;j++)   
  {
  if(fabs(Hex_Matrix_Element[c+i][c+j])>=1.0e-10)
  {
   tem=GBL_Matrix_Index[s+i];
   {rval=Srch_NZ_Element_Column(s+i,s+j);}
   if(rval>=0)
   {
    bck=tem;tem=rval;
    GBL_Matrix_ColNos[s+i][tem]=s+j;
    GBL_Matrix_Data[s+i][tem]+=Hex_Matrix_Element[c+i][c+j];
    tem=bck;
    if(tem==0) tem++; 
    GBL_Matrix_Index[s+i]=tem;
   }
   else
   {
   GBL_Matrix_ColNos[s+i][tem]=s+j;
   GBL_Matrix_Data[s+i][tem]=Hex_Matrix_Element[c+i][c+j];
   tem++;
   GBL_Matrix_Index[s+i]=tem;
   }
  }
  }

for(y=0;y<=5;y++)
{
 x=connect1[k][y];  if(x==-1) continue; 

 l=HexHed[cn][x]; l--; z=l*3;qa=x*3;

 for(i=0;i<=2;i++)
  for(j=0;j<=2;j++)
   {
   if(Hex_Matrix_Element[c+i][qa+j]!=0.0)
   {
   tem=GBL_Matrix_Index[s+i];
   {fval=Srch_NZ_Element_Column(s+i,z+j);}
   if(fval>=0)
   {
    bck=tem;tem=fval;
    GBL_Matrix_ColNos[s+i][tem]=z+j;
    GBL_Matrix_Data[s+i][tem]+=Hex_Matrix_Element[c+i][qa+j];
    tem=bck;
    GBL_Matrix_Index[s+i]=tem;
   }
   else
   {
   GBL_Matrix_ColNos[s+i][tem]=z+j;
   GBL_Matrix_Data[s+i][tem]=Hex_Matrix_Element[c+i][qa+j];
   tem++;
   GBL_Matrix_Index[s+i]=tem;
   }
  }
 }
}
ma=0;qa=0;c=0;s=0;z=0;x=0;
}
}
}

/***************************************************************************
*
*   FUNCTION        : Read_Input_Pass_1()
*   DESCRIPTION     : find the number of nodes and the number of hexahedron.
*
****************************************************************************/
 
void Read_Input_Pass_1()
{
     int    x1, y1, z1, x2, y2, z2, OutF_Count=0, output_flag=0;