emap5.c

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

              JetaPQ and JzetaPQ. 
Input value: 
    int *TrngleNode --- nodes of the source triangle 
    double *MidPoint --- middle point of the observing triangle 
Return value: none 
    Global value used: TotQuadPoint, Qpnt,  WaveNumber 
Global value modified: none
Subroutines called: ComputeGaussQuadPoint(), VTXmag(),
    COMplex_Cmplx(), Real_Mul(), COMplex_Mul(), COMplex_Add(),
    COMplex_Sub(), COMplex_Expon(), Aimag() 
*****************************************************************************/
void ComputeNonSingul1(int *TrngleNode, double *MidPoint, complex *JnoPQ,
                       complex *JxiPQ, complex *JetaPQ, complex *JzetaPQ )
{
    double SrcPoint[3], Rcp;
    int QuadPointN;

   
    for(QuadPointN=0;QuadPointN<TotQuadPoint;QuadPointN++)
    {
        ComputeGaussQuadPoint(QuadPointN,TrngleNode,SrcPoint);
        Rcp = VTXmag(MidPoint,SrcPoint);

        if( fabs(Rcp) <= 1.0e-07 )
        {
            fprintf(OutF1,"Check observation point \n"); 
            exit(1);
        }  /* if( fabs(Rcp) */
        else
        {
            complex Eprj, CFactor;

            CFactor = COMplex_Cmplx(0.0,(Rcp*WaveNumber));
            Eprj  = Real_Mul((Wght[QuadPointN]/(Rcp*Rcp*Rcp)),
                      COMplex_Cmplx((1.0+Real(CFactor)),Aimag(CFactor)));
            Eprj  = COMplex_Mul(Eprj,COMplex_Expon(-1.0,CFactor));
            *JnoPQ = COMplex_Add(*JnoPQ, Eprj);
            *JxiPQ = COMplex_Add(*JxiPQ,  Real_Mul(Qpnt[QuadPointN][0],Eprj) );
            *JetaPQ = COMplex_Add(*JetaPQ, Real_Mul(Qpnt[QuadPointN][1],Eprj) );

        }  /* else */

    }   /* for(QuadPoint */

    *JzetaPQ = COMplex_Sub(*JnoPQ,COMplex_Add(*JxiPQ, *JetaPQ) );

}
 


 
/****************************************************************************
Prototype:    void    ComputeParameters() 
Description:    To computes parameters(TetVolume, TrngleArea, TrngleCentroid
                and etc) which are necessary to build the FEM matrix 
                and MOM matrix, 
Input value:    none 
Return value:    none 
Global value used:     AbsPermeable, AbsPermitt, CFactor1, CFactor2,
    EdgeLength, FreeSpaceVel, ImpeDance, NodeCord, OperateFreq,
    GlobalEdgeEnds, TetVolume, TotEdgeNum, TotTetElement, TotTrnlgeNum,
    TParam, TrngleArea, TrngleCentroid, WaveLength, WaveNumber, 
Global value modified:     CFactor1, CFactor2, FreeSpaceVel,
    TetVolume, TrngleArea, TrngleCentroid 
Subroutines called:     COMplex_Cmplx(), ComputeTetraHeronVolume(),
                        ComputeTrngleCentroid(), Real_Mul(),  
                        TetfaceAreaNormal(), VTXmag() 
****************************************************************************/
void ComputeParameters()
{
    int II;

    OperateFreq   = OperateFreq * 1.0E+06;

    FreeSpaceVel  = 3.0E+08;
    AbsPermeable  = 1.25663706144E-06;
    AbsPermitt    = 8.8542E-12;

    WaveLength = FreeSpaceVel/OperateFreq;
    WaveNumber = 2.0*M_PI /WaveLength;
    ImpeDance  = WaveNumber/(2.0*M_PI*OperateFreq*AbsPermitt);

    CFactor1 = Real_Mul((-1.0/(8.0*M_PI)),
                        COMplex_Cmplx(0.0,(WaveNumber*ImpeDance)));
    CFactor2 = Real_Mul((-1.0/(2.0*M_PI)),
                        COMplex_Cmplx(0.0,-(ImpeDance/WaveNumber)));


    TParam[0]  = COMplex_Cmplx(0.0,-(1.0/(ImpeDance*WaveNumber)));
    TParam[1]  = COMplex_Cmplx(0.0,(WaveNumber/ImpeDance));
 
    TetVolume = double_Vector(TotTetElement);
    for(II=0;II<TotTetElement;II++) 
        ComputeTetraHedronVolume(II,NodeCord,TetNode,TetVolume);
 
    EdgeLength = double_Vector(TotEdgeNum);
    for(II=0;II<TotEdgeNum;II++)  
        EdgeLength[II] = VTXmag(NodeCord[GlobalEdgeEnds[II][1]-1],
                         NodeCord[GlobalEdgeEnds[II][0]-1]);

    TrngleCenTroid = double_Matrix(TotTrngleNum,3);
    for(II=0;II<TotTrngleNum;II++)
        ComputeTrngleCentroid(II,NodeCord,TrngleNode,TrngleCenTroid);

    TrngleArea     = double_Vector(TotTrngleNum);
    TrngleNormal= double_Matrix(TotTrngleNum,3);
 
 
    for(II=0;II<TotTrngleNum;II++) 
        TetFaceAreaNormal(II,NodeCord,TrngleNode,TrngleNormal,TrngleArea);
 
}




/****************************************************************************
Prototype:  void  ComputeTetraHedronVolume(int TetHedNum, double**Cord, 
                                 int **TGNodeNum,double *TVolume ) 
Description:  To calculate the volume of a tetrahedron 
Input value: 
    int TetHedNum --- the number of  the tetrahedron 
    double **cord  --- the global node table 
    int **TGNodeNum --- the global node numbers of the tetrahedron 
    double *TVolume, -- where to store the results. 
Return value:     none 
Global value used:     none 
Global value modified:     none 
Subroutines called:     VTXsub1(), VTXcross() 
****************************************************************************/
void ComputeTetraHedronVolume(int TetHedNum, double **Cord, int **TGNodeNum,  
                              double *TVolume)
{
    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]=
                Cord[TGNodeNum[TetHedNum][Count_i]-1][Count_j];

    VTXsub1(1,0,Buff,Buff1);
    VTXsub1(2,0,Buff,Buff2);
    VTXsub1(3,0,Buff,Buff3);
    VTXcross(Buff1,Buff2,Buff4);

    TVolume[TetHedNum] =Buff3[0]*Buff4[0]+Buff3[1]*Buff4[1]+Buff3[2]*Buff4[2];

}
 



/****************************************************************************
Prototype: void ComputeTrngleCentroid( int TrngleNum, double **Cord,
                                 int **BFaceNode, double **CenTroid) 
Description: To compute the centroid of a triangular patch. 
Input value: 
    int TrngleNum  --- the number of the triangle 
    double **Cord  --- the gloabl node table 
    int BFaceNode  --- the nodes of the triangle 
    double **Centroid -- where to store the results 
Return value:    none 
Global value used:    none 
Global value modified:    none 
Subroutines called:     VTXadd2() 
*****************************************************************************/
void ComputeTrngleCentroid(int TrngleNum, double **Cord, int **BFaceNode, 
                           double **CenTroid )
{
    int Count_k; double Buff[3],Buff1[3],Buff2[3],Buff3[3];
 
    for(Count_k=0;Count_k<=2;Count_k++)
    {
        Buff1[Count_k]=Cord[BFaceNode[TrngleNum][0]-1][Count_k];
        Buff2[Count_k]=Cord[BFaceNode[TrngleNum][1]-1][Count_k];
        Buff3[Count_k]=Cord[BFaceNode[TrngleNum][2]-1][Count_k];
    }

    VTXadd2(Buff1,Buff2,Buff3,Buff);
    CenTroid[TrngleNum][0] = (1.0/3.0) * Buff[0];
    CenTroid[TrngleNum][1] = (1.0/3.0) * Buff[1];
    CenTroid[TrngleNum][2] = (1.0/3.0) * Buff[2];

}
 






/****************************************************************************
Prototype: void TetFaceAreaNormal(int TrngleNum, double **Cord, int **BFaceNode,
                                  double **Normal, double *Area) 
Description: To compute the unit normal and the area of each face of the 
             tetrahedron 
Input value: 
    int TrngleNum --- index of the tetrahedron
    double **Cord --- the global node table
    int **BFaceNode  --- the nodes of each face
    double **Normal -- where to store the centroid of each face
    double *Area  --- where to store the area of each face
Return value:    none 
Global value used:     none 
Global value modified:     none 
Subroutines called:     VTXcross1(), VTXsub1(), VTXcross()
****************************************************************************/
void TetFaceAreaNormal(int TrngleNum, double **Cord, int **BFaceNode, 
                        double **Normal, double *Area)
{
    int j,k;
    double Buff[3][3],Buff1[3],Buff2[3],Buff3[3],MagAreaFace;
    double FaceNorm[3];

    for(j=0;j<=2;j++)
        for(k=0;k<=2;k++)
            Buff[j][k]= Cord[BFaceNode[TrngleNum][j]-1][k];
 
    VTXcross1(1,2,Buff,Buff1);
    VTXcross1(2,0,Buff,Buff2);
    VTXcross1(0,1,Buff,Buff3);

    FaceNorm[0] = (Buff1[0] + Buff2[0] + Buff3[0]);
    FaceNorm[1] = (Buff1[1] + Buff2[1] + Buff3[1]);
    FaceNorm[2] = (Buff1[2] + Buff2[2] + Buff3[2]);
    MagAreaFace = sqrt( FaceNorm[0]*FaceNorm[0] +
               FaceNorm[1]*FaceNorm[1] + FaceNorm[2]*FaceNorm[2]);
    Normal[TrngleNum][0] = FaceNorm[0] / MagAreaFace ;
    Normal[TrngleNum][1] = FaceNorm[1] / MagAreaFace ;
    Normal[TrngleNum][2] = FaceNorm[2] / MagAreaFace ;

    VTXsub1(0,1,Buff,Buff1);
    VTXsub1(2,1,Buff,Buff2);
    VTXcross(Buff1,Buff2,Buff3);
    Area[TrngleNum] = 0.5*sqrt(VTXdot(Buff3,Buff3));

}
  

  
/****************************************************************************
Prototype: void ComputeSourceVector(double WaveNum, 
                                    double*MidPoint, complex *EiC) 
Description:  To compute the E fields of an incident plane wave at the 
              boundary surface. 
Input value: 
    double WaveNum --- wavenumber of the plane wave. 
    double MidPoint --- midpoint of the observing triangle. 
    complex *EiC  --- where to store the results 
Return value:  none 
Global value used:  none 
Global value modified:  none 
Subroutines called:   none *****************************************************************************/
void ComputeSourceVector(double WaveNum, double *MidPoint, complex *EiC)
{
    double R,PolK[3],PolE[3];
    complex RE;

    /* the Progagation vector PolK[] in coordinate system */

    PolK[0] = sin(K_Theta*M_PI/180.0)*cos(K_Phi*M_PI/180.0);
    PolK[1] = sin(K_Theta*M_PI/180.0)*sin(K_Phi*M_PI/180.0);
    PolK[2] = cos(K_Theta*M_PI/180.0);

    /* the E polarization vector PolE[] in coordiante system */
    PolE[0] = sin(E_Theta*M_PI/180.0)*cos(E_Phi*M_PI/180.0);
    PolE[1] = sin(E_Theta*M_PI/180.0)*sin(E_Phi*M_PI/180.0);
    PolE[2] = cos(E_Theta*M_PI/180.0);

    /* RE=e(-jk*R)  vector productor, the phase term for every point */
    R = PolK[0]*MidPoint[0] + PolK[1]*MidPoint[1] + PolK[2]*MidPoint[2];
    RE= COMplex_Expon(-1.0,COMplex_Cmplx(0.0,(R*WaveNum)));
 
    /* E field including the phase term */
    EiC[0] = Real_Mul(PolE[0],RE);
    EiC[1] = Real_Mul(PolE[1],RE);
    EiC[2] = Real_Mul(PolE[2],RE);

}  /* end of ComputeSourceVector( ) */
 



/***************************************************************************
Prototype: void InvertMatrix(complex **Qmat, int MatrixSize) 
Description:    To invert a complex matrix. The results are stored in the 
                input matrix, thus no additional memory needed. 
Input value: 
    complex **Qmat --- a complex matrix 
    int matrixSize --- the size of the matrix 
Return value:     none 
Global value used:     none 
Global value modified:     none 
Subroutines called:    COMplex_Abs(), COMplex_Div(), COMplex_Cmplx(),
                       COMplex_Sub(), COMplex_Mul(), INT_Matrix(), 
                       free_INT_Matrix() 
****************************************************************************/
void InvertMatrix(complex **Qmat, int MatrixSize)
{
    int **Inter,i,j,k,l,KK,JJ,Row,Col;
    double Val1,Val2; complex tmp;

 
    Inter  = INT_Matrix(MatrixSize,2);

    for(i=0;i<MatrixSize;i++) 
        for(j=0;j<2;j++)  Inter[i][j]=0;
 
    for(k=0;k<MatrixSize;k++)  
    {
        JJ=k;
        if( k!=MatrixSize-1 )
        {
            KK= k+1; 
            Val1 = COMplex_Abs(Qmat[k][k]);
            for(i=KK;i<MatrixSize;i++)
            {
                Val2 = COMplex_Abs(Qmat[i][k]);
                if( (Val1 - Val2)<0.0 ) 
                {
                    Val1 = Val2; 
                    JJ=i;
                }     /* if( (Val1-Val2) */
            }         /* for(i=KK */
  
        }   /* if( k!=  */

        Inter[k][0] = k; Inter[k][1] = JJ;

        if( JJ!=k )
        {
            for(j=0;j<MatrixSize;j++)
            {
                tmp = Qmat[JJ][j]; 
                Qmat[JJ][j] = Qmat[k][j]; 
                Qmat[k][j] = tmp;

            }   /* for(j=0 */

            for(j=0;j<MatrixSize;j++)
                if(j!=k) Qmat[k][j]=COMplex_Div(Qmat[k][j],Qmat[k][k] );
            
        }    /* if( JJ!=K ) */
        else
        {
            for(j=0;j<MatrixSize;j++)
                if(j!=k)Qmat[k][j] = COMplex_Div( Qmat[k][j],Qmat[k][k] );
            
        }   /* else */


        Qmat[k][k] = COMplex_Div(COMplex_Cmplx(1.0,0.0),Qmat[k][k]);