sift5.c

三维矢量有限元-矩量法电磁场分析程序。 EMAP5 is a full-wave electromagnetic field solver that combines the method of m

        { 
            DielBoundEdgeNum++;  
            HybrdBoundEdgeNum++;  
            BoundEdgeStat[kk]=ii+1; 
            kk++;
        }  /* if() */


    /* the following are half-metal/half-dielectric edges */
    for(ii=0;ii<TotEdgeNum;ii++)
        if( (EdgeStat[ii]==-1) && (MetalEdgeStat[ii]==1) )
        {   
            HybrdBoundEdgeNum++; 
            BoundEdgeStat[kk]=ii+1; 
            kk++; 
        }  /* if() */


    /* following are external metal edges */
    for(ii=0;ii<TotEdgeNum;ii++)
        if( (EdgeStat[ii]==-1) && (MetalEdgeStat[ii]==2) )
        {
            BoundEdgeStat[kk]=ii+1; 
            kk++;
        }  /* if() */

    /* this function can't be called before the boudary edge table
       is generated  */
    if (SourceType=='V')  DetectVSourceEdges();


    PrintOutput();
    return 0;

}   /* end of main() */



/**************************************************************************
Prototype:    void    ComputeParameters() 
Description:     To compute parameters for the mesh. For example, the 
                 total node  number, the total triangular element number, the 
                 total tetrahedron number. 
Input value:     none 
Return value:     none 
Global value used:     TotCubeNumEdge, TotCubeTrngleNum, TotCubeNumNode,
    TotNumHexHedron, TotTetElement, TotEdgeNum, TotNodeNum. TotTrngleNum,
    TotDipolNumEdge, TotDipolTrngleNum, Xdim, Ydim, Zdim, ExdcondNum, 
    Xextcond, Yextcond, DipolWid, DipolLen, Dipolplane_YZ, Dipolplane_XY,     
    TotInnerEdgeNum, InnerEdgeStat 
Global value modified:    TotCubeNumEdge, TotCubeTrngleNum, TotCubeNumNode,
    TotNumHexHedron, TotTetElement, TotEdgeNum, TotNodeNum. TotTrngleNum,
    TotInnerEdgeNum, InnerEdgeStat 
Subroutines called:     INT_Vector() 
***************************************************************************/
void ComputeParameters()
{

    int i;
    double error=1e-5;


    XdiM=YdiM=ZdiM=0;

    for(i=0; i<CellUnitXNum; i++)
        XdiM+=(int) ( (CellUnitX[i].end -CellUnitX[i].start)/CellUnitX[i].step +error );

    for(i=0; i<CellUnitYNum; i++)
        YdiM+=(int) ( (CellUnitY[i].end-CellUnitY[i].start)/CellUnitY[i].step +error );

    for(i=0; i<CellUnitZNum; i++)
        ZdiM+= (int) ( (CellUnitZ[i].end-CellUnitZ[i].start)/CellUnitZ[i].step +error );



    /* number of cubes */
    TotNumHexHedron = XdiM*YdiM*ZdiM;
    /* each cube can be devided into 5 tetrahedrons*/
    TotTetElement = 5 * TotNumHexHedron;

    /* number of edge in the dielectirc box */
    TotCubeNumEdge = 6*XdiM*YdiM*ZdiM + (XdiM+YdiM+ZdiM)
                    + 3*(XdiM*YdiM+YdiM*ZdiM+ZdiM*XdiM);

    /* number of triangles contained in the dielectirc box */
    TotCubeTrngleNum= 4*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM);
    TotCubeNumNode  = (XdiM+1)*(YdiM+1)*(ZdiM+1);
 
    /* TotEdgeNum, TotNodeNum and TotTrngleNum will also includes 
       those elements of external conductors */ 
    TotEdgeNum  = TotCubeNumEdge;
    TotNodeNum  = TotCubeNumNode;
    TotTrngleNum= TotCubeTrngleNum;


    /**************External Conductor information*************************/
    for(IEdge=0;IEdge<ExtcondNum;IEdge++)
    {
        /* This external conductor is a plane parallel to Y-Z plane */
        if( Xextcond[0][IEdge]==Xextcond[1][IEdge] )
        {
            DipolWid[IEdge]= (int) fabs((double)Yextcond[1][IEdge]-Yextcond[0][IEdge]) 
                            /DelYextcond[IEdge];
            DipolLen[IEdge]= (int)  fabs((double)Zextcond[1][IEdge]-Zextcond[0][IEdge])
                            /DelZextcond[IEdge];
            Dipolplane_YZ[IEdge]=1;

        }



        /* This external conductor is a plane parallel to X-Z plane */
        else if (Yextcond[0][IEdge]==Yextcond[1][IEdge])
        {
            DipolWid[IEdge]= (int) fabs((double)Xextcond[1][IEdge]-Xextcond[0][IEdge]) 
                             /DelXextcond[IEdge];
            DipolLen[IEdge]= (int) fabs((double)Zextcond[1][IEdge]-Zextcond[0][IEdge])
                             /DelZextcond[IEdge];
            Dipolplane_ZX[IEdge]=1;
        }


        /* This external conductor is a plane parallel to X-Y plane */
        else if (Zextcond[0][IEdge]==Zextcond[1][IEdge])
        {
            DipolWid[IEdge]= (int) (fabs((double)Xextcond[1][IEdge]-Xextcond[0][IEdge]) 
                             /DelXextcond[IEdge] );
            DipolLen[IEdge]= (int) (fabs((double)Yextcond[1][IEdge]-Yextcond[0][IEdge])
                             /DelYextcond[IEdge]);
            Dipolplane_XY[IEdge]=1;
        }

        else
        {
            fprintf(stderr, "Plate Orientation is not parallel to Axis.\n");  
            exit(1);
        }


        /* edges on the dipole */
        TotDipolNumEdge[IEdge] = 3*DipolLen[IEdge]*DipolWid[IEdge]
                       + DipolWid[IEdge] + DipolLen[IEdge];
                       
        /* Total edges equal to edges in the dielectrical box and
           edges on dipoles  */                      
        TotEdgeNum  = TotEdgeNum + TotDipolNumEdge[IEdge];

        /* total nodes on the dipole */
        TotDipolNumNode[IEdge] = (DipolLen[IEdge]+1)*(DipolWid[IEdge]+1);

        /* total nodes equal to nodes in the diectrical box and 
           nodes on the dipole */
        TotNodeNum  = TotNodeNum + TotDipolNumNode[IEdge];

        /* calculate total triangle number */
        TotDipolTrngleNum[IEdge]= 2*DipolLen[IEdge]*DipolWid[IEdge];
        TotTrngleNum= TotTrngleNum + TotDipolTrngleNum[IEdge]; 
 
    }  /* for(IEdge=0;IEdge<ExtcondNum;  */


    TotInnerEdgeNum = TotCubeNumEdge - 6*(XdiM*YdiM+YdiM*ZdiM+ZdiM*XdiM);
    InnerEdgeStat = INT_Vector(TotInnerEdgeNum);
 
}   /* end of ComputeParameter() */




/**************************************************************************
Prototype:    void    AllocateMemory() 
Description:     To dynamically allocate required memory. 
Input value:     none 
Return value:     none 
Global value used:     TotNumHexHedron, TotTetElement, TotTrngleNum,
    TotEdgeNum, TotNodeNum, MediaPerm, EdgeStat, MetalEdgeStat, EdgeFound,
    MultFound, PlusTrngleIndex, MinusTrngleIndex, GlobalEdgeEnds HexNode,
    HexEdgeNum, TetNode, TetEdge, TetBoundFaceTetHedronNum,
    TetBoundFaceFaceNum, TrngleNode, TrngleEdge, FaceNode, FaceEdge, EdgeStat,
    MetalEdgeStat, PlusTrngleDetect, PlusTrngleIndex, MinusTrngleDetect,
    MinusTrngleIndex, EdgeFound, MultFound ,GlobalEdgeEnds, NodeCord 
Global value modified:     HexNode, HexEdgeNum, TetNode, TetEdge,
    TetBoundFaceTetHedronNum, TetBoundFaceFaceNum, TrngleNode, TrngleEdge,
    FaceNode, FaceEdge, EdgeStat, MetalEdgeStat, PlusTrngleDetect,    
    PlusTrngleIndex, MinusTrngleDetect, MinusTrngleIndex, EdgeFound, 
    MultFound ,GlobalEdgeEnds, NodeCord 
Subroutines called:     INT_Matrix(), FLOAT_Matrix(), INT_Vector() 
**************************************************************************/
void AllocateMemory()
{

    int ii;

    HexNode = INT_Matrix(TotNumHexHedron,8);
    HexEdgeNum = INT_Matrix(TotNumHexHedron,19);
    MediaPerm   = double_Matrix(TotTetElement,2);
    TetNode  = INT_Matrix(TotTetElement,4);
    TetEdge  = INT_Matrix(TotTetElement,6);

    TetBoundFaceTetHedronNum = INT_Vector(TotTrngleNum);
    TetBoundFaceFaceNum = INT_Vector(TotTrngleNum);

    TrngleNode = INT_Matrix(TotTrngleNum,3);
    TrngleEdge = INT_Matrix(TotTrngleNum,3);

    FaceNode = INT_Matrix(TotTrngleNum,3);
    FaceEdge = INT_Matrix(TotTrngleNum,3);
    EdgeStat = INT_Vector(TotEdgeNum);
    MetalEdgeStat = INT_Vector(TotEdgeNum);
    PlusTrngleDetect = INT_Matrix(TotEdgeNum,3);
    PlusTrngleIndex  = INT_Vector(TotEdgeNum);
    MinusTrngleDetect = INT_Matrix(TotEdgeNum,3);
    MinusTrngleIndex  = INT_Vector(TotEdgeNum);
    EdgeFound  = INT_Vector(TotEdgeNum);
    MultFound  = INT_Vector(TotEdgeNum);
    GlobalEdgeEnds = INT_Matrix(TotEdgeNum,6);
    NodeCord = double_Matrix(TotNodeNum,3);

    /* the default permittivity is 1.0 + j0, which means in free space */
    for(ii=0;ii<TotTetElement;ii++) 
    {  
        MediaPerm[ii][0]=1.0;
        MediaPerm[ii][1]=0.0;
    }

    for(ii=0;ii<TotEdgeNum;ii++) 
    {
        EdgeStat[ii]=1; 
        MetalEdgeStat[ii]=0;
    }


    Null_INTVector( EdgeFound, TotEdgeNum );
    Null_INTVector( MultFound, TotEdgeNum );
    
    Null_INTVector( PlusTrngleIndex, TotEdgeNum );
    Null_INTVector( MinusTrngleIndex, TotEdgeNum );
    Null_INTMatrix( PlusTrngleDetect, TotEdgeNum, 3 );
    Null_INTMatrix( MinusTrngleDetect, TotEdgeNum, 3 );
    Null_INTMatrix( GlobalEdgeEnds, TotEdgeNum, 6 );

}



/**************************************************************************
Prototype:    void    AssignBoxCoordinates() 
Description:     To assign Cartesian coordinates for each node of 
                 each hexahedron 
Input value:     none 
Return value:     none 
Global value used:     NodeCord, Divisor 
Global value modified:     NodeCord 
Subroutines called: none 
**************************************************************************/
void AssignBoxCoordinates()
{
    int i, j, k, n, NodeNum;
    double Xtemp, Ytemp, Ztemp, DivisorX, DivisorY, DivisorZ;;

    
    NodeNum=0; 
     
    /* steps along three axises, revise here on Apr 10 */ 
    DivisorX = CellDimension*CellUnitX[0].step;
    DivisorY = CellDimension*CellUnitY[0].step;
    DivisorZ = CellDimension*CellUnitZ[0].step;
     
    /* start point */
    Xtemp    = CellDimension*XBoundDim[0]; 
    Ytemp    = CellDimension*YBoundDim[0]; 
    Ztemp    = CellDimension*ZBoundDim[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++;   
          
         
                /* to test whether needs to change step */
                for(n=0; n<CellUnitXNum; n++)
                    if(fabs((double)Xtemp-CellDimension*CellUnitX[n].start)<1e-5) break;
            
                /* change to new X step */  
                if(n!=CellUnitXNum) DivisorX=CellDimension*CellUnitX[n].step;
                    Xtemp = Xtemp + DivisorX; 
            }
        
            /* reset along x axis */
            Xtemp = CellDimension*XBoundDim[0];
            DivisorX = CellDimension*CellUnitX[0].step;
         
            /* increase 1 step along Y axis */
            /* revise here on Apr 10 */
            for(n=0; n<CellUnitYNum; n++)
                if(fabs((double)Ytemp-CellDimension*CellUnitY[n].start)<1e-5) break;
            
            /* change to new Y step */  
            if(n!=CellUnitYNum) DivisorY=CellDimension*CellUnitY[n].step;
                Ytemp = Ytemp + DivisorY; 
        
        }
       
        /* reset along y axis */
        Ytemp = CellDimension*XBoundDim[0]; 
        DivisorY = CellDimension*CellUnitY[0].step;
       
        /* increase 1 step along Z axis */
        for(n=0; n<CellUnitZNum; n++)
            if(fabs((double)Ztemp-CellDimension*CellUnitZ[n].start)<1e-5) break;
            
        /* change to new Z step */  
        if(n!=CellUnitZNum) DivisorZ=CellDimension*CellUnitZ[n].step;
          
        Ztemp = Ztemp + DivisorZ; 
    }
}




/**************************************************************************
Prototype:    void    AssignHexHedronNodeNum() 
Description:     To map the  local nodes to the  global nodes for each    
                 hexahedron. 
Input value:     none 
Return value:     none 
Global value used:     TotNumHexHedron, HexNode, XdiM, YdiM 
Global value modified:     none 
Subroutines called:     none **************************************************************************/
void AssignHexHedronNodeNum()
{
   int i, start;

   start=1;
  
   for(i=0;i<TotNumHexHedron;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++;

   }  /* for(i=0  */

}





/**************************************************************************
Prototype:    void    AssignHexHedronEdgeNum() 
Description:     To map the local edges to global edges for each hexahedron 
Input value:  none 
Return value:  none 
Global value used:  XdiM, YdiM, HexEdgeNum 
Global value modified:  none 
Subroutines called:  none 
**************************************************************************/
void AssignHexHedronEdgeNum()
{
    int EdgeNum,xx1,xx2,xx3,kk,dy1,dy2,divz,divz1;
  
    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);