📄 emap5.c
字号:
} /* for(i=0 */
} /* resistor */
/**************CREATE RHS VECTOR ********************/
RHSVector = CMPLX_Vector(TotMatrixEqnNum);
Null_CMPLXVector(RHSVector,TotMatrixEqnNum);
end_time=time(0);
fprintf(OutF1, "\n\tcomputing A matrix\n");
fprintf(OutF1, "\t\ttime used: %d sec\n", (end_time-mid_time));
mid_time=end_time;
/*revise here */
CreateRHSVector(TotInnerEdgeNum,HybrdBoundEdgeNum,TotISourceEdgeNum,
ISourceEdgeMag,RHSVector);
/**************SOLVING HYBRID MATRIX ********************/
ConjugateSolver(TotMatrixEqnNum, GBLmatColIndex,GBLmatColAddr, GBLmat,
RHSVector,TotMatrixEqnNum, TOLERANCE);
end_time=time(0);
fprintf(OutF1, "\n\tsolving the hybrid matrix equation\n");
fprintf(OutF1, "\t\ttime used: %d sec\n", (end_time-mid_time));
mid_time=end_time;
/******* COMPUTATION OF SURFACE FIELDS **************/
EdVector = CMPLX_Vector(HybrdBoundEdgeNum);
JdVector = CMPLX_Vector(HybrdBoundEdgeNum);
JcVector = CMPLX_Vector(TotExtMetalEdgeNum);
Null_CMPLXVector(EdVector,HybrdBoundEdgeNum);
Null_CMPLXVector(JdVector,HybrdBoundEdgeNum);
Null_CMPLXVector(JcVector,TotExtMetalEdgeNum);
SurfaceFieldCompute();
end_time=time(0);
fprintf(OutF1,"\n\tcomputing surface currents:\n");
fprintf(OutF1, "\t\ttime used: %d sec\n", (end_time-mid_time));
mid_time=end_time;
PrintOutput();
end_time=time(0);
fprintf(OutF1, "\n total time used: %d sec\n",
(end_time-start_time) );
fprintf(OutF1, "\n%s\n", log_mesg);
fclose(OutF1);
fclose(InF);
return 0;
} /* end of main() */
/****************************************************************************
Prototype: int SearchNonZeroElement(int RowNum, int ColNum)
Description: To search an element in the global FEM matrix. To conserve
memory, only none-zero elements are stored in the FEM matrix
using row-indexed scheme. If the element is found, return the
column index to the global FEM matrix. Otherwise, return -1.
Input value:
int RowNum --- row number, which is assigned to the observing edge
int ColNum --- column number, which is assigned to the source edge.
Return value: the column index to GBLmat.
Global value used: int *GBLmatColIndex, int *GBLmatColAddr
Global value modified: none
Subroutines called: none
*****************************************************************************/
int SearchNonZeroElement(int RowNum,int ColNum)
{
int Count_i,Count_j;
for(Count_i=0;Count_i<GBLmatColIndex[RowNum];Count_i++)
if( ColNum==(Count_j=GBLmatColAddr[RowNum][Count_i]) ){
return Count_i;
break;
}
return -1;
}
/****************************************************************************
Prototype: void ComputeGaussQuadPoint(int QuadPoint, int *TrngleNode,
double *SrcPointCol)
Description: To compute the coordinates of 7-point Gauss nodes of
a triangular patch.
Input value:
int QuadPoint --- node index, it can be from 0 to 6.
int *TrngleNode --- the three nodes of a tringular patch.
double *SrcPointCol --- where to store the results
Return value: none
Global value used: NodeCord, Qpnt
Global value modified: none
Subroutines called: none
Note: Not very sure. *****************************************************************************/
void ComputeGaussQuadPoint(int QuadPoint, int *TrngleNode, double *SrcPointCol)
{
SrcPointCol[0] = Qpnt[QuadPoint][0]*NodeCord[TrngleNode[0]][0]
+ Qpnt[QuadPoint][1]*NodeCord[TrngleNode[1]][0]
+ Qpnt[QuadPoint][2]*NodeCord[TrngleNode[2]][0];
SrcPointCol[1] = Qpnt[QuadPoint][0]*NodeCord[TrngleNode[0]][1]
+ Qpnt[QuadPoint][1]*NodeCord[TrngleNode[1]][1]
+ Qpnt[QuadPoint][2]*NodeCord[TrngleNode[2]][1];
SrcPointCol[2] = Qpnt[QuadPoint][0]*NodeCord[TrngleNode[0]][2]
+ Qpnt[QuadPoint][1]*NodeCord[TrngleNode[1]][2]
+ Qpnt[QuadPoint][2]*NodeCord[TrngleNode[2]][2];
}
/****************************************************************************
Prototype: void ComputeSingul(int *Node, double *MidPoint, double AREA)
Description: To evaluate the singularity of the MOM integral analytically when
the observing triangle coincides with the source triangle. The
results are stored in global variables RealIno, RealIxi and
RealIeta.
Input value:
int *Node --- the nodes of the triangle
double *MidPoint --- middle point of the triangle
double AREA --- the area of the triangle
Return value: none
Global value used: double **NodeCord
Global value modified: none
Subroutines called: none
*****************************************************************************/
void ComputeSingul(int *Node, double *MidPoint, double AREA,
double *RealIno, double *RealIxi, double *RealIeta)
{
double RIO,RIP,RIQ;
double RO[3],RJ[3],RHO[3],HATL[3][3],SMLRO[3][3];
double X,Y,Z,R,RX,RY,RZ,A,B,C,D,E,F,ABE,ACF,BDF,RJ1,RJ2,RJ3,RJ4,RD;
int I1,I2,I3,L,L1,K,K1;
I1 = Node[0];
I2 = Node[1];
I3 = Node[2];
RIO = 0.0;
RIP = 0.0;
RIQ = 0.0;
RX = MidPoint[0];
RY = MidPoint[1];
RZ = MidPoint[2];
X = (NodeCord[I2][1]-NodeCord[I1][1])*(NodeCord[I3][2]-NodeCord[I1][2])
- (NodeCord[I3][1]-NodeCord[I1][1])*(NodeCord[I2][2]-NodeCord[I1][2]);
Y = (NodeCord[I3][0]-NodeCord[I1][0])*(NodeCord[I2][2]-NodeCord[I1][2])
- (NodeCord[I2][0]-NodeCord[I1][0])*(NodeCord[I3][2]-NodeCord[I1][2]);
Z = (NodeCord[I2][0]-NodeCord[I1][0])*(NodeCord[I3][1]-NodeCord[I1][1])
- (NodeCord[I3][0]-NodeCord[I1][0])*(NodeCord[I2][1]-NodeCord[I1][1]);
D = fabs( X*(RX-NodeCord[I1][0])+Y*(RY-NodeCord[I1][1])
+ Z*(RZ-NodeCord[I1][2]))/sqrt(X*X+Y*Y+Z*Z);
for(L=0;L<=2;L++)
{
K = Node[L];
if( L<2 ) K1 = Node[L+1];
if( L==2 ) K1 = Node[0];
R = sqrt( pow( (NodeCord[K1][0]-NodeCord[K][0]), 2.0)
+ pow( (NodeCord[K1][1]-NodeCord[K][1]), 2.0)
+ pow( (NodeCord[K1][2]-NodeCord[K][2]), 2.0) );
HATL[L][0]=(NodeCord[K1][0]-NodeCord[K][0])/R;
HATL[L][1]=(NodeCord[K1][1]-NodeCord[K][1])/R;
HATL[L][2]=(NodeCord[K1][2]-NodeCord[K][2])/R;
X=(RY-NodeCord[K][1])*(RZ-NodeCord[K1][2])
-(RY-NodeCord[K1][1])*(RZ-NodeCord[K][2]);
Y=(RX-NodeCord[K1][0])*(RZ-NodeCord[K][2])
-(RX-NodeCord[K][0])*(RZ-NodeCord[K1][2]);
Z=(RX-NodeCord[K][0])*(RY-NodeCord[K1][1])
-(RX-NodeCord[K1][0])*(RY-NodeCord[K][1]);
RO[L] = sqrt(X*X+Y*Y+Z*Z)/R;
RHO[L]= sqrt(RO[L]*RO[L]-D*D);
RJ[L] = sqrt( pow( (RX-NodeCord[K][0]), 2.0)
+ pow( (RY-NodeCord[K][1]), 2.0)
+ pow( (RZ-NodeCord[K][2]), 2.0) );
R = (RX-NodeCord[K][0])*HATL[L][0]
+ (RY-NodeCord[K][1])*HATL[L][1]
+ (RZ-NodeCord[K][2])*HATL[L][2];
SMLRO[L][0]=NodeCord[K][0]+R*HATL[L][0];
SMLRO[L][1]=NodeCord[K][1]+R*HATL[L][1];
SMLRO[L][2]=NodeCord[K][2]+R*HATL[L][2];
}
RIO = -2.0*M_PI*D;
for(L=0;L<=2;L++)
{
L1=L+1;
if(L == 2) L1=0;
K=Node[L];
if( L<2 ) K1=Node[L+1];
if( L==2 ) K1=Node[0];
X = HATL[L][0]*(NodeCord[K1][0]-SMLRO[L][0])
+ HATL[L][1]*(NodeCord[K1][1]-SMLRO[L][1])
+ HATL[L][2]*(NodeCord[K1][2]-SMLRO[L][2]);
Y = HATL[L][0]*(NodeCord[K][0]-SMLRO[L][0])
+ HATL[L][1]*(NodeCord[K][1]-SMLRO[L][1])
+ HATL[L][2]*(NodeCord[K][2]-SMLRO[L][2]);
RIO = RIO+(1.0/2.0)*RHO[L]*
log((RJ[L1]+X)/(RJ[L1]-X)*(RJ[L]-Y)/(RJ[L]+Y) )
+ D*asin(D*X/RO[L]/sqrt(RJ[L1]*RJ[L1]-D*D))
- D*asin(D*Y/RO[L]/sqrt(RJ[L] *RJ[L] -D*D));
}
RIO = RIO/(2.0*AREA);
A = pow( (NodeCord[I1][0]-NodeCord[I3][0]), 2.0 )
+ pow( (NodeCord[I1][1]-NodeCord[I3][1]), 2.0 )
+ pow( (NodeCord[I1][2]-NodeCord[I3][2]), 2.0 );
B = pow( (NodeCord[I2][0]-NodeCord[I3][0]), 2.0 )
+ pow( (NodeCord[I2][1]-NodeCord[I3][1]), 2.0 )
+ pow( (NodeCord[I2][2]-NodeCord[I3][2]), 2.0 );
C=-2.0*( (RX-NodeCord[I3][0])*(NodeCord[I1][0]-NodeCord[I3][0])
+(RY-NodeCord[I3][1])*(NodeCord[I1][1]-NodeCord[I3][1])
+(RZ-NodeCord[I3][2])*(NodeCord[I1][2]-NodeCord[I3][2]));
D=-2.0*( (RX-NodeCord[I3][0])*(NodeCord[I2][0]-NodeCord[I3][0])
+(RY-NodeCord[I3][1])*(NodeCord[I2][1]-NodeCord[I3][1])
+(RZ-NodeCord[I3][2])*(NodeCord[I2][2]-NodeCord[I3][2]));
E= 2.0*( (NodeCord[I1][0]-NodeCord[I3][0])*(NodeCord[I2][0]-NodeCord[I3][0])
+(NodeCord[I1][1]-NodeCord[I3][1])*(NodeCord[I2][1]-NodeCord[I3][1])
+(NodeCord[I1][2]-NodeCord[I3][2])*(NodeCord[I2][2]-NodeCord[I3][2]));
F= pow((RX-NodeCord[I3][0]), 2.0)
+ pow((RY-NodeCord[I3][1]), 2.0)
+ pow((RZ-NodeCord[I3][2]), 2.0);
ABE=sqrt(A+B-E);
ACF=sqrt(A+C+F);
BDF=sqrt(B+D+F);
RJ1 = ( (2*B-C+D-E)*BDF+(2*A+C-D-E)*ACF )/4.0/(A+B-E)
+( 4*(A+C)*(B+D+F)+4*F*(B-C-E)-(C+D+E)*(C+D+E) )/(8.0*(A+B-E)*ABE)
*log(fabs( (2*ABE*BDF+(2*B-C+D-E))/(2*ABE*ACF-(2*A+C-D-E)) ));
RJ2 = ( (2*B+D)*BDF-D*sqrt(F) )/(4.0*B)+(4*B*F-D*D)/(8*B*sqrt(B))
*log(fabs( (2*sqrt(B)*BDF+2*B+D)/(2*sqrt(B*F)+D) ));
RJ3= ( (2*A+C-D-E)*ACF+(2*B-C+D-E)*BDF )/(4.0*(A+B-E))
+( 4*(A+C)*(B+D+F)+4*F*(B-C-E)-(C+D+E)*(C+D+E) )/(8.0*(A+B-E)*ABE)
*log(fabs( (2*ABE*ACF+(2*A+C-D-E))/(2*ABE*BDF-(2*B-C+D-E)) ));
RJ4 = ( (2*A+C)*ACF-C*sqrt(F) )/(4.0*A)+(4*A*F-C*C)
/(8*A*sqrt(A))*log(fabs( (2*sqrt(A)*ACF+2*A+C)/(2*sqrt(A*F)+C) ));
RD=4*A*B-E*E;
RIP=( 4*B*(RJ1-RJ2)-2*E*(RJ3-RJ4)-(2*B*C-E*D)*RIO )/RD;
RIQ=( 4*A*(RJ3-RJ4)-2*E*(RJ1-RJ2)-(2*A*D-E*C)*RIO )/RD;
*RealIno = RIO;
*RealIxi = RIP;
*RealIeta = RIQ;
}
/****************************************************************************
Prototype: void ComputeNonSingul(int *TrngleNode, double One,
double *MidPoint, complex *Ino_PQ,
complex *Ixi_PQ,
complex *Ieta_PQ, complex *Izeta_PQ)
Description: To evaluate the MOM surface integral numerically when the
observation triangle and the source triangle are different.
Input value:
int *TrngleNode --- nodes of the source triangle
double One --- a factor
double *MidPoint --- middle point of the observing trianlge
Return value:
complex *Ino_PQ --- to store the value of Ino
complex *Ixi_PQ --- to store the value of Ixi
complex *Ieta_PQ --- to store the value of Ieta
complex *Izeta_PQ --- to store the value of Izeta.
The above values are used in ComputeCMatrix().
Global value used: none
Global value modified: none
Subroutines called: COMplex_add(), COMplex_Cmplx(), COMplex_Expon(),
Real_Mul(), COMplex_Sub()
****************************************************************************/
void ComputeNonSingul(int *TrngleNode, double One, double *MidPoint,
complex *Ino_PQ, complex *Ixi_PQ, complex *Ieta_PQ,
complex *Izeta_PQ)
{
double SrcPoint[3], Rcp;
int QuadPointN;
complex Eprj, V1;
/* initialize the four complex variables to zero */
for(QuadPointN=0;QuadPointN<TotQuadPoint;QuadPointN++)
{
ComputeGaussQuadPoint(QuadPointN,TrngleNode,SrcPoint);
Rcp = VTXmag(MidPoint,SrcPoint);
if( fabs(Rcp) <= 1.0e-07 )
{
Eprj = COMplex_Cmplx(0.0, -(Wght[QuadPointN]*WaveNumber) );
*Ino_PQ = COMplex_Add(*Ino_PQ, Eprj);
*Ixi_PQ = COMplex_Add(*Ixi_PQ, Real_Mul(Qpnt[QuadPointN][0],Eprj));
*Ieta_PQ = COMplex_Add(*Ieta_PQ,
Real_Mul(Qpnt[QuadPointN][1],Eprj) );
} /* if( fabs(Rcp */
else
{
V1= COMplex_Cmplx(0.0,(Rcp*WaveNumber));
V1= COMplex_Expon(-1.0,V1);
V1 = COMplex_Cmplx( (Real(V1) - One),Aimag(V1));
Eprj = Real_Mul( (Wght[QuadPointN]/Rcp),V1);
*Ino_PQ = COMplex_Add(*Ino_PQ, Eprj);
*Ixi_PQ = COMplex_Add(*Ixi_PQ, Real_Mul(Qpnt[QuadPointN][0],Eprj) );
*Ieta_PQ = COMplex_Add(*Ieta_PQ,
Real_Mul(Qpnt[QuadPointN][1],Eprj) );
} /* else */
} /* for(QuadPointN */
*Izeta_PQ = COMplex_Sub(*Ino_PQ,COMplex_Add(*Ixi_PQ , *Ieta_PQ) );
} /* end of ComputeNonSingul( ) */
/****************************************************************************
Prototype: void ComputeNonSingul1(int *TrngleNode, double *MidPoint)
Description: To evaluate the MOM surface integral numerically when the
observation triangle and the source triangle are
different. Results are stored in the global variables JxiPQ,
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -