emap3.c

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

/***************************************************************************
*   FUNCTION        : GlobalEdgeEndsDiv1()
****************************************************************************/

void GlobalEdgeEndsDiv1()
{
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][1]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][0]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][1]=HexNode[HexNum][7];
   }

/***************************************************************************
*   FUNCTION        : GlobalEdgeEndsDiv2()
****************************************************************************/

void GlobalEdgeEndsDiv2()
{
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][1]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][0]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][1]=HexNode[HexNum][7];
}


/***************************************************************************
*   FUNCTION        : Find_Volume_Determinent()
****************************************************************************/

void FindTetraHedronVolume()
{
   int Count_i, Count_j;
   double Buff[4][3],Buff1[3],Buff2[3],Buff3[3],Buff4[3];
   for(Count_i=0;Count_i<=3;Count_i++)
   for(Count_j=0;Count_j<=2;Count_j++)
   Buff[Count_i][Count_j]=NodeCord[TetGlobalNodeNum[t][Count_i]-1][Count_j];
   VTXsub1(1,0,Buff,Buff1);
   VTXsub1(2,0,Buff,Buff2);
   VTXsub1(3,0,Buff,Buff3);
   VTXcross(Buff1,Buff2,Buff4);
   TetVolume=(1.0/6.0)*(Buff3[0]*Buff4[0]+Buff3[1]*Buff4[1]+Buff3[2]*Buff4[2]);
}

/***************************************************************************
*   FUNCTION        : ComputeTetraCoFactor()
****************************************************************************/

void ComputeTetraCoFactor()
   {
   int i,j,RowNum,ColNum,Row,Col,Count=0;
   double CoF[4][4],Buff[4][4],CoFs[4][4]; 
     {
    for(i=0;i<=3;i++)  for(j=0;j<=3;j++) {
    if(j==3) Buff[i][j]= 1.0;
    else     Buff[i][j]= NodeCord[TetGlobalNodeNum[t][i]-1][j];
                                         }
     for(RowNum=0;RowNum<=3;RowNum++)
     for(ColNum=0;ColNum<=3;ColNum++) {
     Row=0;Col=0;Count=0;
        for(i=0;i<=3;i++)
        for(j=0;j<=3;j++) {
        if((i!=RowNum) && (j!=ColNum)) {
        Count++;
        if(Count<4)              Row=0;
        if((Count>3)&&(Count<7)) Row=1;
        else if(Count>6)         Row=2;
        CoF[Row][Col]=Buff[i][j];
        Col++; Col=Col%3;              }
                          }
     Row=0;Col=0;Count=0;
     CoFs[RowNum][ColNum]
     =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((RowNum+ColNum)%2!=0)  CoFs[RowNum][ColNum]*=-1.0;
                                      }
     for(RowNum=0;RowNum<=3;RowNum++) {
     if     (RowNum==0) Col = 3;
     else if(RowNum==1) Col = 0;
     else if(RowNum==2) Col = 1;
     else if(RowNum==3) Col = 2;
     for(ColNum=0;ColNum<=3;ColNum++) {
     CoFactor[RowNum][ColNum] = -CoFs[ColNum][Col];
                                      }}
      }
    }

/***************************************************************************
*   FUNCTION        : S_Matrix()
****************************************************************************/

void S_Matrix()
  {
  int i,j;  double Length_i,Length_j,Mult1,Volume;
   {
   Volume = TetVolume; Mult1 = (1296.0*Volume*Volume*Volume);
   for(i=0;i<=5;i++) for(j=0;j<=5;j++)  {
  Length_i = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[i][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[i][1]]-1]);
  Length_j = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[j][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[j][1]]-1]);
   S_mat[i][j] = ( CoFactor[2][edge_end[i][0]] * CoFactor[3][edge_end[i][1]]
                 - CoFactor[3][edge_end[i][0]] * CoFactor[2][edge_end[i][1]] )
               * ( CoFactor[2][edge_end[j][0]] * CoFactor[3][edge_end[j][1]]
                 - CoFactor[3][edge_end[j][0]] * CoFactor[2][edge_end[j][1]] )
               + ( CoFactor[3][edge_end[i][0]] * CoFactor[1][edge_end[i][1]]
                 - CoFactor[1][edge_end[i][0]] * CoFactor[3][edge_end[i][1]] )
               * ( CoFactor[3][edge_end[j][0]] * CoFactor[1][edge_end[j][1]]
                 - CoFactor[1][edge_end[j][0]] * CoFactor[3][edge_end[j][1]] )
               + ( CoFactor[1][edge_end[i][0]] * CoFactor[2][edge_end[i][1]]
                 - CoFactor[2][edge_end[i][0]] * CoFactor[1][edge_end[i][1]] )
               * ( CoFactor[1][edge_end[j][0]] * CoFactor[2][edge_end[j][1]]
                 - CoFactor[2][edge_end[j][0]] * CoFactor[1][edge_end[j][1]] );
   S_mat[i][j] = (Length_i * Length_j * S_mat[i][j])/ Mult1 ;
                                         }
   }
   }

/***************************************************************************
*   FUNCTION        : T_Matrix()
****************************************************************************/

void T_Matrix()
  {
  int i,j;
  double Mult1,Length_i,Length_j;
 
   Mult1 = 1.0/(720.0 * TetVolume);
   T_mat[0][0] =               ff(1,1) + ff(2,2) - ff(1,2);
   T_mat[0][1] = 2 * ff(2,3) - ff(2,1) - ff(1,3) + ff(1,1);
   T_mat[0][2] = 2 * ff(2,4) - ff(2,1) - ff(1,4) + ff(1,1);
   T_mat[0][3] =-2 * ff(1,3) - ff(2,2) + ff(2,3) + ff(1,2);
   T_mat[0][4] = 2 * ff(1,4) - ff(2,4) - ff(1,2) + ff(2,2);
   T_mat[0][5] =     ff(2,4) - ff(2,3) - ff(1,4) + ff(1,3);
   T_mat[1][1] =               ff(1,1) + ff(3,3) - ff(1,3);
   T_mat[1][2] = 2 * ff(3,4) - ff(1,3) - ff(1,4) + ff(1,1);
   T_mat[1][3] = 2 * ff(1,2) - ff(2,3) - ff(1,3) + ff(3,3);
   T_mat[1][4] =     ff(2,3) - ff(3,4) - ff(1,2) + ff(1,4);
   T_mat[1][5] =-2 * ff(1,4) - ff(3,3) + ff(1,3) + ff(3,4);
   T_mat[2][2] =               ff(1,1) + ff(4,4) - ff(1,4);
   T_mat[2][3] =     ff(3,4) - ff(2,4) - ff(1,3) + ff(1,2);
   T_mat[2][4] =-2 * ff(1,2) - ff(4,4) + ff(1,4) + ff(2,4);
   T_mat[2][5] = 2 * ff(1,3) - ff(3,4) - ff(1,4) + ff(4,4);
   T_mat[3][3] =               ff(3,3) + ff(2,2) - ff(2,3);
   T_mat[3][4] =-2 * ff(3,4) - ff(2,2) + ff(2,3) + ff(2,4);
   T_mat[3][5] =-2 * ff(2,4) - ff(3,3) + ff(2,3) + ff(3,4);
   T_mat[4][4] =               ff(2,2) + ff(4,4) - ff(2,4);
   T_mat[4][5] =-2 * ff(2,3) - ff(4,4) + ff(2,4) + ff(3,4);
   T_mat[5][5] =               ff(3,3) + ff(4,4) - ff(3,4);
   for(i=0;i<=5;i++) for(j=0;j<=i-1;j++)
   {T_mat[i][j] = T_mat[j][i];}
 
   for(i=0;i<=5;i++) for(j=0;j<=5;j++) {
  Length_i = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[i][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[i][1]]-1]);
  Length_j = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[j][0]]-1],
                    NodeCord[TetGlobalNodeNum[t][edge_end[j][1]]-1]);
   if(i==j) T_mat[i][j] = 2 * Length_i * Length_j * Mult1 * T_mat[i][j];
   else     T_mat[i][j] =     Length_i * Length_j * Mult1 * T_mat[i][j];
   T_mat[i][j] = RELPerm[HexNum][0] * T_mat[i][j] ;
   T_mat[i][j] = WaveNumber*WaveNumber * T_mat[i][j] ;
                                       }
    }

/***************************************************************************
*   FUNCTION        : TetraSubMatrix()
****************************************************************************/

void TetraSubMatrix()
  {
  int i,j;
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j]=0.0;}
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j] =4.0*S_mat[i][j]-T_mat[i][j];}
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j] = A_mat[t][i][j] * Sign(TetEdgeNum[t][j]);}
  for(i=0;i<=5;i++) for(j=0;j<=5;j++)
  {A_mat[t][i][j] = A_mat[t][i][j] * Sign(TetEdgeNum[t][i]);}
  }

/***************************************************************************
*   FUNCTION        : FindGlobalMatrix()
****************************************************************************/

void FindGlobalMatrix()
{
 int m,n,k,l,rval;
  {
  for(m=0;m<=5;m++) for(n=0;n<=5;n++)   {
   k=absol(TetEdgeNum[t][m])-1; l=absol(TetEdgeNum[t][n])-1;
   if(fabs(A_mat[t][m][n])>1.0E-10) {
   rval=Srch_NZ_Element_Column(k,l);
   if(rval>=0)                {
   GBL_Matrix_Data[k][rval]+=A_mat[t][m][n];
                              }
   else                       {
   GBL_Matrix_ColNos[k][GBL_Matrix_Index[k]]=l;
   GBL_Matrix_Data[k][GBL_Matrix_Index[k]]=A_mat[t][m][n];
   GBL_Matrix_Index[k]++;
                              }
                                    }

                                        }
  }
}

/***************************************************************************
*   FUNCTION        : CountHalfBandWidth()
****************************************************************************/

void CountHalfBandWidth()
     {
      short int i,j,k,l,ColDiff,Val;
        {
        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;
        if(HalfBandWidth<=(ColDiff+1))
        {HalfBandWidth=(ColDiff+1);}
                   }
                                              }
                                              }