emap3.c

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

 HexNum = k*YdiM*XdiM + j*XdiM + i;
 con = pow(-1.0,(double)i)*pow(-1.0,(double)j)*pow(-1.0,(double)k);

 if(con==1) {HexHedraSubDiv1(); EdgeSubDiv1(); GlobalEdgeEndsDiv1();}
 else       {HexHedraSubDiv2(); EdgeSubDiv2(); GlobalEdgeEndsDiv2();}

 for(t=0;t<=4;t++)
 {
  FindTetraHedronVolume();
  ComputeTetraCoFactor();
  S_Matrix();
  T_Matrix();
  TetraSubMatrix();
  FindGlobalMatrix();
 }

                          }
CountHalfBandWidth();
HalfBandMatrix = FLOAT_Matrix(TotInnerEdgeNum+1,HalfBandWidth+1);
FORC_Matrix_Data=FLOAT_Matrix(TotInnerEdgeNum,FRCD_Matrix_MaxColNos);
FORC_Matrix_ColNos=INT_Matrix(TotInnerEdgeNum,FRCD_Matrix_MaxColNos);
RHSVector = FLOAT_Vector(TotInnerEdgeNum);
for(i=0;i<=TotInnerEdgeNum;i++)
for(j=0;j<=HalfBandWidth;j++) {HalfBandMatrix[i][j] =0.0;}
PartitionGlobalMatrix();
ComputeRHSVector();
BandMatrixSolver();
Produce_Output();
return(0);
}

/***************************************************************************/

void ParameterInfo()
{
double FreeSpaceVel,AbsPermeable,AbsPermitt,WaveLength; 

 AbsPermeable  = 1.25663706144E-06;
 AbsPermitt    = 8.8542E-12;
 FreeSpaceVel  = 1.0/sqrt(AbsPermeable*AbsPermitt);
 WaveLength = FreeSpaceVel/(OperateFreq*1.0E+06);
 WaveNumber = 2.0*M_PI/WaveLength;
 printf("OperateFrequency=%g MHz\n",OperateFreq);
 printf("WaveLength=%f m\n",WaveLength );
 printf("FreeSpaceVelocity=%E m/s\n",FreeSpaceVel );
 printf("WaveNumber=%g\n",WaveNumber);

 TotEdgeNum      = 6*XdiM*YdiM*ZdiM + (XdiM+YdiM+ZdiM)
                 + 3*(XdiM*YdiM+YdiM*ZdiM+ZdiM*XdiM);
 TotNodeNum      = (XdiM+1)*(YdiM+1)*(ZdiM+1);
 TotNumHexHedron = XdiM*YdiM*ZdiM;
 TotTetElement = 5 * TotNumHexHedron;
 TotTrngleNum    = 4*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM);
 TotBoundEdgeNum = 6*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM);
 TotBoundNodeNum = 2*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM+1);
 DivisorX = CellDimension; 
 DivisorY = CellDimension; 
 DivisorZ = CellDimension; 

 printf("CellDimension = %f m \n",CellDimension);
 printf("Total Edge Number = %d \n",TotEdgeNum);
 printf("Total Node Number = %d \n",TotNodeNum);
 printf("Total Cube Number = %d \n",TotNumHexHedron);
 printf("Total Tetrahedron Number = %d \n",TotTetElement);
 printf("Total Boundary Edge Number = %d \n",TotBoundEdgeNum);
 printf("Total Boundary Node Number = %d \n",TotBoundNodeNum);

  /* Local edge orientation of tetrahedron edges */
  edge_end[0][0] = 0; edge_end[0][1] = 1;
  edge_end[1][0] = 0; edge_end[1][1] = 2;
  edge_end[2][0] = 0; edge_end[2][1] = 3;
  edge_end[3][0] = 1; edge_end[3][1] = 2;
  edge_end[4][0] = 3; edge_end[4][1] = 1;
  edge_end[5][0] = 2; edge_end[5][1] = 3;

 }

/***************************************************************************
*   FUNCTION        : AssignGlobalCoorDinates()
****************************************************************************/

void AssignGlobalCoorDinates()
 {
  int i,j,k,NodeNum; double Xtemp,Ytemp,Ztemp;
   {
     NodeNum=0; Xtemp=0.0; Ytemp=0.0; Ztemp=0.0;
     for(k=0;k<=ZdiM;k++) {
       for(j=0;j<=YdiM;j++) {
        for(i=0;i<=XdiM;i++) {
          NodeCord[NodeNum][0] = Xtemp;
          NodeCord[NodeNum][1] = Ytemp;
          NodeCord[NodeNum][2] = Ztemp;
          NodeNum++;   Xtemp = Xtemp + DivisorX; }
          Xtemp = 0.0; Ytemp = Ytemp + DivisorY; }
          Ytemp = 0.0; Ztemp = Ztemp + DivisorZ; }
    }
 }

/***************************************************************************
*   FUNCTION        : AssignHexHedronNodeNumbering()
****************************************************************************/

void AssignHexHedronNodeNumbering()
{
 int i;
  {
   start=1;
   for(i=0;i<=TotNumHexHedron-1;i++)
     {
      HexNode[i][0]=start;
      HexNode[i][1]=HexNode[i][0]+1;
      HexNode[i][2]=HexNode[i][0]+XdiM+1;
      HexNode[i][3]=HexNode[i][2]+1;
      HexNode[i][4]=HexNode[i][0]+(XdiM+1)*(YdiM+1);
      HexNode[i][5]=HexNode[i][4]+1;
      HexNode[i][6]=HexNode[i][4]+(XdiM+1);
      HexNode[i][7]=HexNode[i][6]+1;
      if(((i+1)%((YdiM*XdiM))==0)) {start=start+(XdiM+3);}
      else                     {start=start+1;}
      if((start%(XdiM+1))==0)  { start++; }
     }
  }
}

/***************************************************************************
*   FUNCTION        : AssignHexHedronEdgeNumbering()
****************************************************************************/

void AssignHexHedronEdgeNumbering()
  {
  int EdgeNum,xx1,xx2,xx3,kk,dy1,dy2,divz,divz1,i,j,k;
  {
        for(k=0;k<=ZdiM-1;k++) 
        for(j=0;j<=YdiM-1;j++)
        for(i=0;i<=XdiM-1;i++) {  
        EdgeNum = k*YdiM*XdiM + j*XdiM + i;
        for(kk=0;kk<=17;kk++)   {HexEdgeNum[EdgeNum][kk] = 0;}
        divz=(6*XdiM*YdiM+3*XdiM+3*YdiM+1)*k;
        divz1=  6*XdiM*YdiM + 3*XdiM + 3*YdiM + 1 ;
        dy1 = (3*XdiM+1)*j; dy2 = (3*XdiM+2)*j;
        xx1 = 1 + (YdiM+1)*XdiM + (2*XdiM+1)*YdiM;
        xx2 = xx1 + (2*XdiM+1);
        xx3 = xx2 + (XdiM + 1);
        HexEdgeNum[EdgeNum][0] = 1 + i + dy1 + divz;
        HexEdgeNum[EdgeNum][1] = 1 + XdiM + 2*i + dy1 + divz;
        HexEdgeNum[EdgeNum][2] = HexEdgeNum[EdgeNum][1] + 1;  
        HexEdgeNum[EdgeNum][3] = HexEdgeNum[EdgeNum][2] + 1;
        HexEdgeNum[EdgeNum][4] = 1 + 3*XdiM + 1 + i + dy1 + divz;
        HexEdgeNum[EdgeNum][5] = xx1 + 2*i + dy2 + divz;
        HexEdgeNum[EdgeNum][6] = HexEdgeNum[EdgeNum][5] + 1; 
        HexEdgeNum[EdgeNum][7] = HexEdgeNum[EdgeNum][6] + 1;
        HexEdgeNum[EdgeNum][8] = xx2 + i + dy2 + divz;
        HexEdgeNum[EdgeNum][9] = HexEdgeNum[EdgeNum][8] + 1;
        HexEdgeNum[EdgeNum][10]= xx3 + 2*i + dy2 + divz;
        HexEdgeNum[EdgeNum][11]= HexEdgeNum[EdgeNum][10] + 1;
        HexEdgeNum[EdgeNum][12]= HexEdgeNum[EdgeNum][11] + 1; 
        HexEdgeNum[EdgeNum][13]= HexEdgeNum[EdgeNum][0] + divz1;
        HexEdgeNum[EdgeNum][14]= HexEdgeNum[EdgeNum][1] + divz1; 
        HexEdgeNum[EdgeNum][15]= HexEdgeNum[EdgeNum][2] + divz1;
        HexEdgeNum[EdgeNum][16]= HexEdgeNum[EdgeNum][3] + divz1; 
        HexEdgeNum[EdgeNum][17]= HexEdgeNum[EdgeNum][4] + divz1;
                               }
        
 }
}

/***************************************************************************
*   FUNCTION        : HexHedraSubDiv1()
****************************************************************************/

void HexHedraSubDiv1()
{
TetGlobalNodeNum[0][0]=HexNode[HexNum][0];TetGlobalNodeNum[1][0]=HexNode[HexNum][3];
TetGlobalNodeNum[0][1]=HexNode[HexNum][3];TetGlobalNodeNum[1][1]=HexNode[HexNum][5];
TetGlobalNodeNum[0][2]=HexNode[HexNum][2];TetGlobalNodeNum[1][2]=HexNode[HexNum][7];
TetGlobalNodeNum[0][3]=HexNode[HexNum][6];TetGlobalNodeNum[1][3]=HexNode[HexNum][6];
 
TetGlobalNodeNum[2][0]=HexNode[HexNum][0];TetGlobalNodeNum[3][0]=HexNode[HexNum][0];
TetGlobalNodeNum[2][1]=HexNode[HexNum][1];TetGlobalNodeNum[3][1]=HexNode[HexNum][6];
TetGlobalNodeNum[2][2]=HexNode[HexNum][3];TetGlobalNodeNum[3][2]=HexNode[HexNum][4];
TetGlobalNodeNum[2][3]=HexNode[HexNum][5];TetGlobalNodeNum[3][3]=HexNode[HexNum][5];
 
TetGlobalNodeNum[4][0]=HexNode[HexNum][0];TetGlobalNodeNum[4][2]=HexNode[HexNum][3];
TetGlobalNodeNum[4][1]=HexNode[HexNum][5];TetGlobalNodeNum[4][3]=HexNode[HexNum][6];
}

/***************************************************************************
*   FUNCTION        : HexHedraSubDiv2()
****************************************************************************/

void HexHedraSubDiv2()
{
TetGlobalNodeNum[0][0]=HexNode[HexNum][0];TetGlobalNodeNum[1][0]=HexNode[HexNum][1];
TetGlobalNodeNum[0][1]=HexNode[HexNum][1];TetGlobalNodeNum[1][1]=HexNode[HexNum][7];
TetGlobalNodeNum[0][2]=HexNode[HexNum][2];TetGlobalNodeNum[1][2]=HexNode[HexNum][4];
TetGlobalNodeNum[0][3]=HexNode[HexNum][4];TetGlobalNodeNum[1][3]=HexNode[HexNum][5];
 
TetGlobalNodeNum[2][0]=HexNode[HexNum][1];TetGlobalNodeNum[3][0]=HexNode[HexNum][2];
TetGlobalNodeNum[2][1]=HexNode[HexNum][7];TetGlobalNodeNum[3][1]=HexNode[HexNum][6];
TetGlobalNodeNum[2][2]=HexNode[HexNum][3];TetGlobalNodeNum[3][2]=HexNode[HexNum][4];
TetGlobalNodeNum[2][3]=HexNode[HexNum][2];TetGlobalNodeNum[3][3]=HexNode[HexNum][7];

TetGlobalNodeNum[4][0]=HexNode[HexNum][1];TetGlobalNodeNum[4][2]=HexNode[HexNum][4];
TetGlobalNodeNum[4][1]=HexNode[HexNum][2];TetGlobalNodeNum[4][3]=HexNode[HexNum][7];
}

/***************************************************************************
*   FUNCTION        : EdgeSubDiv1()
****************************************************************************/

void EdgeSubDiv1()
{
TetEdgeNum[0][0] = HexEdgeNum[HexNum][2];  
TetEdgeNum[0][1] = HexEdgeNum[HexNum][1]; 
TetEdgeNum[0][2] = HexEdgeNum[HexNum][8];    
TetEdgeNum[0][3] =-HexEdgeNum[HexNum][4];  
TetEdgeNum[0][4] =-HexEdgeNum[HexNum][11];  
TetEdgeNum[0][5] = HexEdgeNum[HexNum][10];   

TetEdgeNum[1][0] = HexEdgeNum[HexNum][9];   
TetEdgeNum[1][1] = HexEdgeNum[HexNum][12];  
TetEdgeNum[1][2] = HexEdgeNum[HexNum][11]; 
TetEdgeNum[1][3] = HexEdgeNum[HexNum][16];  
TetEdgeNum[1][4] =-HexEdgeNum[HexNum][15];  
TetEdgeNum[1][5] =-HexEdgeNum[HexNum][17];   
 
TetEdgeNum[2][0] = HexEdgeNum[HexNum][0];   
TetEdgeNum[2][1] = HexEdgeNum[HexNum][2];   
TetEdgeNum[2][2] = HexEdgeNum[HexNum][6]; 
TetEdgeNum[2][3] = HexEdgeNum[HexNum][3];  
TetEdgeNum[2][4] =-HexEdgeNum[HexNum][7];  
TetEdgeNum[2][5] = HexEdgeNum[HexNum][9]; 
 
TetEdgeNum[3][0] = HexEdgeNum[HexNum][8];   
TetEdgeNum[3][1] = HexEdgeNum[HexNum][5];  
TetEdgeNum[3][2] = HexEdgeNum[HexNum][6]; 
TetEdgeNum[3][3] =-HexEdgeNum[HexNum][14];  
TetEdgeNum[3][4] = HexEdgeNum[HexNum][15];  
TetEdgeNum[3][5] = HexEdgeNum[HexNum][13]; 
 
TetEdgeNum[4][0] = HexEdgeNum[HexNum][6];   
TetEdgeNum[4][1] = HexEdgeNum[HexNum][2];   
TetEdgeNum[4][2] = HexEdgeNum[HexNum][8];  
TetEdgeNum[4][3] =-HexEdgeNum[HexNum][9]; 
TetEdgeNum[4][4] =-HexEdgeNum[HexNum][15];  
TetEdgeNum[4][5] = HexEdgeNum[HexNum][11];  
}

/***************************************************************************
*   FUNCTION        : EdgeSubDiv2()
****************************************************************************/

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