sift5.c

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

                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;

            }  /* for(i=0 */
 
}




/**************************************************************************
Prototype:    void    DetectInnerEdgeStatus() 
Description:     To find edges located within the FEM region. 
Input value: none 
Return value: none 
Global value used: XdiM, YdiM, EdgeStat, HexEdgeNum, 
Global value modified: none 
Subroutines called: none ***************************************************************************/
void DetectInnerEdgeStatus()
{
 
     /* int EdgeNum,xx1,xx2,xx3,kk,dy1,dy2,divz,divz1; */
     int EdgeNum, i,j,k,m;
  
     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;
         
                 /* set default value to 1, which means it is an inner edge */
                 for(m=0; m<18; m++)  EdgeStat[HexEdgeNum[EdgeNum][m]-1] = 1;

                 /* detect whether it is a boundary edge */
                 if( (j==0)||(k==0) )  EdgeStat[HexEdgeNum[EdgeNum][0]-1] = -1;

                 if( (i==0)||(k==0) )  EdgeStat[HexEdgeNum[EdgeNum][1]-1] = -1;

                 if(k==0)  EdgeStat[HexEdgeNum[EdgeNum][2]-1] = -1;

                 if( (i==XdiM-1)||(k==0) )  
                     EdgeStat[HexEdgeNum[EdgeNum][3]-1] = -1;

                 if( (j==YdiM-1)||(k==0) )  
                     EdgeStat[HexEdgeNum[EdgeNum][4]-1] = -1;

                 if( (j==0)||(i==0) )   
                     EdgeStat[HexEdgeNum[EdgeNum][5]-1] = -1;

                 if( j==0 ) EdgeStat[HexEdgeNum[EdgeNum][6]-1] = -1;

                 if( (i==XdiM-1)||(j==0) )  
                     EdgeStat[HexEdgeNum[EdgeNum][7]-1] = -1;

                 if( i==0 ) EdgeStat[HexEdgeNum[EdgeNum][8]-1] = -1;

                 if( i==XdiM-1 )  EdgeStat[HexEdgeNum[EdgeNum][9]-1] = -1;

                 if( (j==YdiM-1)||(i==0) )  
                     EdgeStat[HexEdgeNum[EdgeNum][10]-1] = -1;

                 if( j==(YdiM-1) )   
                     EdgeStat[HexEdgeNum[EdgeNum][11]-1] = -1;

                 if( (i==XdiM-1)||(j==YdiM-1) )    
                       EdgeStat[HexEdgeNum[EdgeNum][12]-1] = -1;

                 if( (k==ZdiM-1)||(j==0) )   
                       EdgeStat[HexEdgeNum[EdgeNum][13]-1] = -1;

                 if( (k==ZdiM-1)||(i==0) )  
                       EdgeStat[HexEdgeNum[EdgeNum][14]-1] = -1;

                 if(k==ZdiM-1)  
                        EdgeStat[HexEdgeNum[EdgeNum][15]-1] = -1;

                 if((i==XdiM-1)||(k==ZdiM-1) )   
                        EdgeStat[HexEdgeNum[EdgeNum][16]-1] = -1;

                 if((j==YdiM-1)||(k==ZdiM-1) )  
                        EdgeStat[HexEdgeNum[EdgeNum][17]-1] = -1;
                             
            }  /* for() */
}
 
 

 
/**************************************************************************
Prototype:    void    DetectBoundEdgeStatus() 
Description:     To find edges located on the surface of the FEM region. 
Input value:     none 
Return value:     none 
Global value used:    XdiM, ZdiM, YdiM, EdgeStat, HexEdgeNum, BndcondNum, 
                      Xbndcond, Ybndcond, Zbndcond, MetalEdgeStat, EdgeStat,  
                      EdgeStat 
Global value modified: none
Subroutines called:     none
**************************************************************************/
void DetectBoundEdgeStatus()
{
    int EdgeNum;
    int  FrontWall, BackWall, BottomWall, TopWall, RightWall, LeftWall; 


    /*********************************************************
     assume the dielectric box is located in the orginal
     point the coordinate system is as follows:
  
                Y                    
                ^                     
                |                     
          X <--(x) Z                   
    *********************************************************/
  
    for(PatchNum=0;PatchNum<BndcondNum;PatchNum++)
    {
        FrontWall=0;BackWall=0;
        BottomWall=0;TopWall=0;
        RightWall=0;LeftWall=0;
 
        if( Zbndcond[0][PatchNum]==Zbndcond[1][PatchNum] )
        {
           if(Zbndcond[0][PatchNum]==ZBoundDim[0]) FrontWall=1;
           else   BackWall=1;
        }

        else if(Ybndcond[0][PatchNum]==Ybndcond[1][PatchNum])
        {
            if( Ybndcond[0][PatchNum]==YBoundDim[0] )  BottomWall=1;
            else    TopWall=1;
        }

        else if(Xbndcond[0][PatchNum]==Xbndcond[1][PatchNum]) 
        {
            if( Xbndcond[0][PatchNum]==XBoundDim[0] )  RightWall=1;
            else    LeftWall=1;
        }
         
        else 
        {
            fprintf(stderr,"Check conductor dimension\n");
            exit(1);    
        }
    
        /* revise here on Apr. 10 due to introducing non-uniform mesh */

        XmaxWall=YmaxWall= ZmaxWall=0;
        XminWall=YminWall= ZminWall=0;
   
        /* calculate XminWall */
        if( RightWall==1 )  XminWall=XmaxWall=0;
        else if( LeftWall==1 ) XminWall=XmaxWall=XdiM;
        else 
        {
            for(i=0; i<CellUnitXNum; i++)
                if( Xbndcond[0][PatchNum]>(CellUnitX[i].start-1e-5 ) &&
                    Xbndcond[0][PatchNum]<(CellUnitX[i].end +1e-5) )  break;

            if( i==CellUnitXNum )
            {
                fprintf(stderr,"Please check conductor's X dimension\n");
                exit(1);
            }
   
            for(j=0; j<i; j++)
                XminWall+= (int ) ( fabs((double)CellUnitX[j].end-CellUnitX[j].start) 
                           /CellUnitX[j].step+1e-5 );
            XminWall+= (int) ( (Xbndcond[0][PatchNum]-CellUnitX[i].start) 
                       /CellUnitX[i].step+1e-5 );
     
            /* calculate XmaxWall */
            for(i=0; i<CellUnitXNum; i++)
                if( Xbndcond[1][PatchNum]>(CellUnitX[i].start-1e-5) &&
                    Xbndcond[1][PatchNum]<(CellUnitX[i].end+1e-5 )) break;
            if( i==CellUnitXNum )
            {
                fprintf(stderr,"Check conductor's X dimension\n");
                exit(1);
            }
   
            for(j=0; j<i; j++)
                XmaxWall+= (int) ((CellUnitX[j].end-CellUnitX[j].start) 
                            /CellUnitX[j].step+1e-3 );
            XmaxWall += (int) ( (Xbndcond[1][PatchNum]-CellUnitX[i].start)
                        /CellUnitX[i].step+1e-3 );
        }/* else */
       
        if( BottomWall==1 )  YminWall=YmaxWall=0;
        else if(TopWall==1 ) YminWall=YmaxWall=YdiM;
        else 
        {    
            for(i=0; i<CellUnitYNum; i++)
                if( Ybndcond[0][PatchNum]>(CellUnitY[i].start-1e-5 ) &&
                    Ybndcond[0][PatchNum]<(CellUnitY[i].end +1e-5) )  break;
            if( i==CellUnitYNum )
            {
                fprintf(stderr,"Check conductor's Y dimension\n");
                exit(1);
            }
   
            for(j=0; j<i; j++)
                YminWall+= (int) ( fabs((double)CellUnitY[j].end-CellUnitY[j].start) 
                           /CellUnitY[j].step+1e-5 );

            YminWall += (int) ( (Ybndcond[0][PatchNum]-CellUnitY[i].start) 
                        /CellUnitY[i].step+1e-5 );
     
 
            /* calculate YmaxWall */
     
            for(i=0;i<CellUnitYNum; i++)
                if( Ybndcond[1][PatchNum]>(CellUnitY[i].start-1e-5) &&
                    Ybndcond[1][PatchNum]<(CellUnitY[i].end+1e-5 ) ) break;

            if( i==CellUnitYNum )
            {
                fprintf(stderr,"Check conductor's Y dimension\n");
                exit(1);
            }
   
            for(j=0; j<i; j++)
                YmaxWall+= (int ) ( fabs ((double)CellUnitY[j].end-CellUnitY[j].start) 
                        /CellUnitY[j].step+1e-5 );
            YmaxWall += (int) ( (Ybndcond[1][PatchNum]-CellUnitY[i].start) 
                        /CellUnitY[i].step+1e-5 );
     
        } /* else */
   
   
        /* determine ZminWall and ZmaxWall  */
        if( FrontWall==1 ) ZminWall=ZmaxWall=0;
        else if( BackWall==1 ) ZminWall=ZmaxWall=ZdiM;
        else
        {  
            for(i=0; i<CellUnitZNum; i++)
                if( Zbndcond[0][PatchNum]>(CellUnitZ[i].start-1e-5) &&
                    Zbndcond[0][PatchNum]<(CellUnitZ[i].end +1e-5) ) break;
            if( i==CellUnitZNum )
            {
                fprintf(stderr,"Check conductor's Z dimension\n");
                exit(1);
            }
   
            for(j=0; j<i; j++)
                ZminWall+= (int ) ( fabs((double)CellUnitZ[j].end-CellUnitZ[j].start) 
                          /CellUnitZ[j].step+1e-5 );
            ZminWall += (int ) ( (Zbndcond[0][PatchNum]-CellUnitZ[i].start) 
                          /CellUnitZ[i].step+1e-5 );
     
 
            /* calculate ZmaxWall */
            for(i=0; i<CellUnitZNum; i++)
                if( Zbndcond[1][PatchNum]>(CellUnitZ[i].start-1e-5) &&
                    Zbndcond[1][PatchNum]<(CellUnitZ[i].end+1e-5) ) break;

            if( i==CellUnitZNum )
            {
                fprintf(stderr,"Please check conductor's Z dimension\n");
                exit(1);
            }
   
            for(j=0; j<i; j++)
                ZmaxWall+= (int ) ( fabs((double)CellUnitZ[j].end-CellUnitZ[j].start)/
                         CellUnitZ[j].step+1e-3 );

            ZmaxWall += (int ) ( (Zbndcond[1][PatchNum]-CellUnitZ[i].start)/ 
                         CellUnitZ[i].step+1e-3 );

        } /* end of else */    


        if( BottomWall )
        {
 
            j=YminWall;
            for(k=0;k<ZdiM;k++)
                for(i=0;i<XdiM;i++) 
                {

                    EdgeNum = k*YdiM*XdiM + j*XdiM + i;
                    if( (i>=XminWall) && (i<XmaxWall) && 
                        (k>=ZminWall) && (k<ZmaxWall) )
                    {
                        MetalEdgeStat[HexEdgeNum[EdgeNum][0]-1]++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][5]-1]++;
                        /* diagonal edge is always external metal edge */
                        MetalEdgeStat[HexEdgeNum[EdgeNum][6]-1] =2;          
                        MetalEdgeStat[HexEdgeNum[EdgeNum][7]-1]++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][13]-1]++;
                    }
                }

        }

        if( TopWall )
        {
            j=YmaxWall-1;
            for(k=0;k<ZdiM;k++)
                for(i=0;i<XdiM;i++) 
                {
                    EdgeNum = k*YdiM*XdiM + j*XdiM + i;
                    if( (i>=XminWall) && (i<XmaxWall) &&  
                        (k>=ZminWall) && (k<ZmaxWall) )
                    {

                        MetalEdgeStat[HexEdgeNum[EdgeNum][4]-1] ++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][10]-1]++;

                        /* diagonal edge */
                        MetalEdgeStat[HexEdgeNum[EdgeNum][11]-1]=2;            
                        MetalEdgeStat[HexEdgeNum[EdgeNum][12]-1]++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][17]-1]++;
                    }  /* if( (!>=XminWall */
                }      /* for(i=0; */
 
        }  /* if(TotWall */

        if( FrontWall )
        {
            k=ZminWall;
            for(j=0;j<YdiM;j++)
                for(i=0;i<XdiM;i++) 
                {
                    EdgeNum = k*YdiM*XdiM + j*XdiM + i;
                    if( (i>=XminWall) && (i<XmaxWall) &&
                        (j>=YminWall) && (j<YmaxWall) )
                    {

                        MetalEdgeStat[HexEdgeNum[EdgeNum][0]-1] ++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][1]-1] ++;
                        /* diagonal edge */
                        MetalEdgeStat[HexEdgeNum[EdgeNum][2]-1] =2;     
                        MetalEdgeStat[HexEdgeNum[EdgeNum][3]-1] ++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][4]-1] ++;
                    }  /* if( (i>=XminWall */
                }      /* for(i=0 */
 
        }

        if( BackWall )
        {
 
            k=ZmaxWall-1;
            for(j=0;j<YdiM;j++)
                for(i=0;i<XdiM;i++) 
                {
                    EdgeNum = k*YdiM*XdiM + j*XdiM + i;
                    if( (i>=XminWall) && (i<XmaxWall) && 
                        (j>=YminWall) && (j<YmaxWall))
                    {

                        MetalEdgeStat[HexEdgeNum[EdgeNum][13]-1] ++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][14]-1] ++; 
                        /* diagonal edge */
                        MetalEdgeStat[HexEdgeNum[EdgeNum][15]-1] =2;          
                        MetalEdgeStat[HexEdgeNum[EdgeNum][16]-1] ++;
                        MetalEdgeStat[HexEdgeNum[EdgeNum][17]-1] ++;
                    }
                }

        }

        if( RightWall )
        {
            i=XminWall;
            for(k=0;k<ZdiM;k++)
                for(j=0;j<YdiM;j++)
                {
                    EdgeNum = k*YdiM*XdiM + j*XdiM + i;
                    if( (j>=YminWall) && (j<YmaxWall) &&  
                        (k>=ZminWall ) && (k<ZmaxWall) )
                    {

 	                MetalEdgeStat[HexEdgeNum[EdgeNum][1]-1]++;