emap5.c

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

        for(i=0;i<MatrixSize;i++)  
        {       
            if(i!=k)
            {
                for(j=0;j<MatrixSize;j++)
                    if( j!=k ) 
                    {      
                        Qmat[i][j]=COMplex_Sub(Qmat[i][j],
                        COMplex_Mul(Qmat[i][k],Qmat[k][j]));
                    }  /* if(j!=k) */              
            }  /* if( i!=k ) */
        }  /* for(i=0; */

        for(i=0;i<MatrixSize;i++)  
            if(i!=k) Qmat[i][k] = Real_Mul(-1.0, 
                          COMplex_Mul(Qmat[i][k],Qmat[k][k]));
           
   }   /* for(k=0;  */

   for(l=0;l<MatrixSize;l++)
   {
       k = MatrixSize - 1 - l;
       Row = Inter[k][0];
       Col = Inter[k][1];
       if( Row != Col ) 
           for(i=0;i<MatrixSize;i++)  
           {
               tmp = Qmat[i][Col]; 
               Qmat[i][Col]=Qmat[i][Row]; 
               Qmat[i][Row]=tmp;

           }   /* for(i=0 */
        
   }   /* for(l=0 */

   free_INT_Matrix(Inter,MatrixSize,2);

}    /*  end of InvertMatrix( )  */
 
 
 
/*****************************************************************************
Prototype: void ComputeCorrectionTerm() 
Description:  To couple the MoM matrix equation to the FEM matrix equation.
              The main idea is to obtain Jd  from C and D. 
              Then couple to the FEM matrix. 
Input value:     none 
Return value:     none 
Global value used:     Ccc, Ccd, Cdc, Cdd, Bdd, TotExtMetalEdgeNum,
                       HybrdBoundEdgeNum, Gd 
Global value modified:   none.  
Subroutines called:     InvertMatrix(), COMplex_Add() 
******************************************************************************/
void ComputeCorrectionTerm()
{
    int II,JJ,KK;
    complex CmpBuff, *TempBuff;


    TempBuff =  CMPLX_Vector(TotBoundEdgeNum);

    Null_CMPLXVector(UdVector,HybrdBoundEdgeNum);

    /*   Ccc <= Ccc^[-1]   */
    InvertMatrix(Ccc,TotExtMetalEdgeNum);
 
   
    for(II=0;II<TotExtMetalEdgeNum;II++) 
        McVector[II] = COMplex_Null();

    /*    Mc <=  Ccc * Vc   */
    for(JJ=0;JJ<TotExtMetalEdgeNum;JJ++)
    {
        CmpBuff = VcVector[JJ];
        for(II=0;II<TotExtMetalEdgeNum;II++)
            McVector[II] =COMplex_Add(McVector[II],  
                                COMplex_Mul(Ccc[II][JJ],CmpBuff));

    }   /* for(JJ=0 */

 
    for(KK=0;KK<HybrdBoundEdgeNum;KK++)
    {
        for(II=0;II<TotExtMetalEdgeNum;II++) 
            TempBuff[II] = COMplex_Null();

        for(JJ=0;JJ<TotExtMetalEdgeNum;JJ++)
        {
            CmpBuff = Ccd[JJ][KK];
            for(II=0;II<TotExtMetalEdgeNum;II++)
                TempBuff[II] =COMplex_Add(TempBuff[II],  
                                    COMplex_Mul(Ccc[II][JJ],CmpBuff));
        }  /* for(JJ=0 */

        for(II=0;II<TotExtMetalEdgeNum;II++) Ccd[II][KK] = TempBuff[II];

    }   /* for(KK=0 */

    for(KK=0;KK<HybrdBoundEdgeNum;KK++)
    {
        for(II=0;II<TotExtMetalEdgeNum;II++) 
            TempBuff[II] = COMplex_Null();

        for(JJ=0;JJ<TotExtMetalEdgeNum;JJ++)
        {
            CmpBuff = Dcd[JJ][KK];
            for(II=0;II<TotExtMetalEdgeNum;II++)
                TempBuff[II] = COMplex_Add(TempBuff[II], 
                                    COMplex_Mul(Ccc[II][JJ],CmpBuff));
        }  /* for(JJ=0 */

        for(II=0;II<TotExtMetalEdgeNum;II++) 
            Dcd[II][KK] = TempBuff[II];

    }   /* for(KK=0 */

    for(II=0;II<HybrdBoundEdgeNum;II++) 
        UdVector[II] = COMplex_Null();

    for(JJ=0;JJ<TotExtMetalEdgeNum;JJ++)
    {
        CmpBuff = McVector[JJ];

        for(II=0;II<HybrdBoundEdgeNum;II++)
            UdVector[II] = COMplex_Add(UdVector[II], 
                                COMplex_Mul(Cdc[II][JJ],CmpBuff));

    }   /* for(JJ=0 */

    /* Ud <=  Cdc * Mc - Vd   */
    for(II=0;II<HybrdBoundEdgeNum;II++) 
        UdVector[II] = COMplex_Sub(UdVector[II],VdVector[II]);

    for(KK=0;KK<HybrdBoundEdgeNum;KK++)
    {
        for(II=0;II<HybrdBoundEdgeNum;II++) 
            TempBuff[II] = COMplex_Null();

        for(JJ=0;JJ<TotExtMetalEdgeNum;JJ++)
        {
            CmpBuff = Ccd[JJ][KK];
            for(II=0;II<HybrdBoundEdgeNum;II++)
                TempBuff[II] =COMplex_Add(TempBuff[II], 
                    COMplex_Mul(Cdc[II][JJ],CmpBuff));

        }  /* for(JJ=0 */

        for(II=0;II<HybrdBoundEdgeNum;II++)
            Cdd[II][KK] = COMplex_Sub(Cdd[II][KK],TempBuff[II]);

     }   /* for(KK=0 */

     for(KK=0;KK<HybrdBoundEdgeNum;KK++)
     {
         for(II=0;II<HybrdBoundEdgeNum;II++) 
             TempBuff[II] = COMplex_Null();

         for(JJ=0;JJ<TotExtMetalEdgeNum;JJ++)
         {
             CmpBuff = Dcd[JJ][KK];

             for(II=0;II<HybrdBoundEdgeNum;II++)
                 TempBuff[II] = COMplex_Add(TempBuff[II], 
                                    COMplex_Mul(Cdc[II][JJ],CmpBuff));

          }   /* for(JJ=0 */

          for(II=0;II<HybrdBoundEdgeNum;II++)
              Ddd[II][KK] = COMplex_Sub(Ddd[II][KK],TempBuff[II]);

     }  /* for(KK=0  */


     /*  Cdd <= Cdd^[-1]   */
     InvertMatrix(Cdd,HybrdBoundEdgeNum);

     for(KK=0;KK<HybrdBoundEdgeNum;KK++)
     {
         for(II=0;II<HybrdBoundEdgeNum;II++) 
             TempBuff[II] = COMplex_Null();
 
         for(JJ=0;JJ<HybrdBoundEdgeNum;JJ++)
         {
             CmpBuff = Ddd[JJ][KK];
             for(II=0;II<HybrdBoundEdgeNum;II++)
                 TempBuff[II] = COMplex_Add( TempBuff[II], 
                       COMplex_Mul(Cdd[II][JJ],CmpBuff));
         }  /* for(JJ=0 */
         
         for(II=0;II<HybrdBoundEdgeNum;II++)
             Ddd[II][KK] = TempBuff[II];

     }  /* for(KK=0 */

     for(II=0;II<HybrdBoundEdgeNum;II++) 
         MdVector[II] = COMplex_Null();

     /*  Md <= Cdd * Ud   */
     for(JJ=0;JJ<HybrdBoundEdgeNum;JJ++)
     {
         CmpBuff = UdVector[JJ];

         for(II=0;II<HybrdBoundEdgeNum;II++)
             MdVector[II] = COMplex_Add(MdVector[II], 
                                   COMplex_Mul(Cdd[II][JJ],CmpBuff));
     }   /* for(JJ=0 */

     for(KK=0;KK<HybrdBoundEdgeNum;KK++)
     {
         for(II=0;II<HybrdBoundEdgeNum;II++) 
             TempBuff[II] = COMplex_Null();

         for(JJ=0;JJ<HybrdBoundEdgeNum;JJ++)
         {
             CmpBuff = Ddd[JJ][KK];

             for(II=0;II<HybrdBoundEdgeNum;II++)
                 if( fabs(Bdd[II][JJ])>=1.0e-30 )
                     TempBuff[II] = COMplex_Add(TempBuff[II], 
                               Real_Mul(Bdd[II][JJ],CmpBuff));
                
         }   /* for(JJ=0  */

         for(II=0;II<HybrdBoundEdgeNum;II++)
             Cdd[II][KK] = Real_Mul(-1.0,TempBuff[II]);

     }  /* for(KK=0  */


     for(II=0;II<HybrdBoundEdgeNum;II++) 
         GdVector[II] = COMplex_Null();

     /*  Gd <= Bdd * Md  */
     for(JJ=0;JJ<HybrdBoundEdgeNum;JJ++)
     {
         CmpBuff = MdVector[JJ];
         for(II=0;II<HybrdBoundEdgeNum;II++)
             GdVector[II] = COMplex_Add(GdVector[II], 
                                 Real_Mul(Bdd[II][JJ],CmpBuff));
     }  /* for(JJ=0 */


     free_CMPLX_Vector(TempBuff, TotBoundEdgeNum);
      
 
}   /* end of ComputeCorrectionTerm */




/****************************************************************************
Prototype:  void  SurfaceFieldCompute() 
Description:  To compute the equivalent surface current density (Jc and Jd ) 
Input value:  none 
Return value:  none 
Global value used:  TotInnerEdgeNum EdVector HybrdBoundEdgeNum, 
                    JdVector, JcVector,TotExtMetalEdgeNum, Dcd, 
                    EdVector, 
Global value modified:  none
Subroutines called:   COMplex_Add(),COMplex_Mul()
*****************************************************************************/
void SurfaceFieldCompute()
{
    int II,JJ;
    complex CmpBuff, *TempBuff;

    /*************************/
    for(II=TotInnerEdgeNum;II<TotMatrixEqnNum;II++)
        EdVector[II-TotInnerEdgeNum] = RHSVector[II];


    TempBuff =  CMPLX_Vector(TotBoundEdgeNum);
    for(II=0;II<HybrdBoundEdgeNum;II++) 
        TempBuff[II] = COMplex_Null();
    for(JJ=0;JJ<HybrdBoundEdgeNum;JJ++)
    {
        CmpBuff = EdVector[JJ];
        for(II=0;II<HybrdBoundEdgeNum;II++)
            TempBuff[II] = COMplex_Add( TempBuff[II], 
                     COMplex_Mul(Ddd[II][JJ],CmpBuff));
    }

    /*  Jd = Ddd *Ed +Md */
    for(II=0;II<HybrdBoundEdgeNum;II++)
        JdVector[II] = COMplex_Add(TempBuff[II],MdVector[II]);

    for(II=0;II<TotExtMetalEdgeNum;II++) 
        TempBuff[II] = COMplex_Null();
    for(JJ=0;JJ<HybrdBoundEdgeNum;JJ++)
    {
        CmpBuff = JdVector[JJ];
        for(II=0;II<TotExtMetalEdgeNum;II++)
        {
            TempBuff[II] = COMplex_Add(TempBuff[II], 
                 COMplex_Mul(Ccd[II][JJ],CmpBuff));
        }
    }
 
    /*  Jc <= - Ccd * Jd - Mc  */
    for(II=0;II<TotExtMetalEdgeNum;II++)
        JcVector[II] = COMplex_Add(TempBuff[II],McVector[II]);
    for(II=0;II<TotExtMetalEdgeNum;II++)
        JcVector[II] = Real_Mul(-1.0,JcVector[II]);

    for(II=0;II<TotExtMetalEdgeNum;II++) 
        TempBuff[II] = COMplex_Null();
    for(JJ=0;JJ<HybrdBoundEdgeNum;JJ++)
    {
        CmpBuff = EdVector[JJ];
        for(II=0;II<TotExtMetalEdgeNum;II++)
        {
            TempBuff[II] = COMplex_Add( TempBuff[II], 
                               COMplex_Mul(Dcd[II][JJ], CmpBuff));
        }
    }

    /*  Jc <=  Dcd * Ed + Jc  */
    for(II=0;II<TotExtMetalEdgeNum;II++)
        JcVector[II] = COMplex_Add(TempBuff[II],JcVector[II]);

    free_CMPLX_Vector(TempBuff, TotBoundEdgeNum);


 }
 




/****************************************************************************
Prototype: void FemMatrixCompute(int NumOfTetElement, int MaxElementPerRow, 
                   int **TetEdge, int **GlobalEdgeEnds, 
                   double **NodeCord, complex *TParam, complex *Epsilon, 
                   int **GBLmatColAddr, complex **GBLmat, 
                   int *GBLmatColIndex, double *TetVolume) 
Description: To build the FEM matrix and store them using row-indexed scheme. 
Input value: 
    int NumOfTetElement --- total number of tetrahedron elements 
    int MaxElementPerRow --- the FEM matrix is a sparse matrix. Every row 
                             has a maximum number (MaxElementPerRow) of 
                             non-zero elemnet. Only non-zero elements are 
                             stored. 
    int **TetEdge --- the edge table of tetrahedron elements 
    int **GlobalEdgeEnds --- global edge table 
    double **NodeCord --- global node table 
    complex *TParam --- store constants used by FEM 
    complex *Epsilon --- the permittivity associated with each 
                          tetrahedron element. 
    int **GBLmatColAddr, complex **GBLmat, int *GBLmatColIndex --- the FEM