📄 emap5.c
字号:
} /*for(IEdge) */
} /* if (!strcmp( */
/*******************************************************/
TotExtMetalEdgeNum = TotBoundEdgeNum - HybrdBoundEdgeNum;
/* revise here on Feb 2 */
TotMatrixEqnNum = TotInnerEdgeNum + HybrdBoundEdgeNum ;
/*********************************************************************
* BOUNDARY ELEMENT PART OF THE HYBRID CODE *
*********************************************************************/
ObsPointArray = double_Matrix(TotQuadPoint,3);
SrcPointArray = double_Matrix(TotQuadPoint,3);
Qpnt = double_Matrix(TotQuadPoint,3);
Wght = double_Vector(TotQuadPoint);
RhoCtrd = double_Matrix2(TotTrngleNum,3,3);
Wght[0]=0.112500000000000;
Wght[2]=Wght[1]= Wght[3]=0.062969590272413;
Wght[6]=Wght[5]= Wght[4]=0.066197076394253;
Qpnt[0][1]=Qpnt[0][0]=0.333333333333333;
Qpnt[2][1]=Qpnt[1][0]=0.797426985353087;
Qpnt[3][0]=Qpnt[3][1]=Qpnt[2][0]=Qpnt[1][1]=0.101286507323456;
Qpnt[6][1]=Qpnt[5][0]=Qpnt[4][1]=Qpnt[4][0]=0.470142064105115;
Qpnt[6][0]=Qpnt[5][1]=0.059715871789770;
for(QuadPointM=0;QuadPointM<TotQuadPoint;QuadPointM++)
Qpnt[QuadPointM][2] = 1.0-Qpnt[QuadPointM][0]-Qpnt[QuadPointM][1];
/*********** malloc menory *********************/
Cdd = CMPLX_Matrix(HybrdBoundEdgeNum,HybrdBoundEdgeNum);
Cdc = CMPLX_Matrix(HybrdBoundEdgeNum,TotExtMetalEdgeNum);
Ccd = CMPLX_Matrix(TotExtMetalEdgeNum,HybrdBoundEdgeNum);
Ccc = CMPLX_Matrix(TotExtMetalEdgeNum,TotExtMetalEdgeNum);
Null_CMPLXMatrix(Cdd,HybrdBoundEdgeNum,HybrdBoundEdgeNum);
Null_CMPLXMatrix(Cdc,HybrdBoundEdgeNum,TotExtMetalEdgeNum);
Null_CMPLXMatrix(Ccd,TotExtMetalEdgeNum,HybrdBoundEdgeNum);
Null_CMPLXMatrix(Ccc,TotExtMetalEdgeNum,TotExtMetalEdgeNum);
Ddd = CMPLX_Matrix(HybrdBoundEdgeNum,HybrdBoundEdgeNum);
Dcd = CMPLX_Matrix(TotExtMetalEdgeNum,HybrdBoundEdgeNum);
Bdd = double_Matrix(HybrdBoundEdgeNum,HybrdBoundEdgeNum);
Null_CMPLXMatrix(Ddd,HybrdBoundEdgeNum,HybrdBoundEdgeNum);
Null_CMPLXMatrix(Dcd,TotExtMetalEdgeNum,HybrdBoundEdgeNum);
Null_double_Matrix(Bdd,HybrdBoundEdgeNum,HybrdBoundEdgeNum);
if( SourceType=='P' )
{
Einc = CMPLX_Vector(TotBoundEdgeNum);
Null_CMPLXVector(Einc,TotBoundEdgeNum);
}
ComputeParameters();
mid_time=start_time=time(0); /* start timer */
for(ObservTrngle=0;ObservTrngle<TotTrngleNum;ObservTrngle++)
{
int RowNum, MultRow,RowRootIndxNum,RowAddlIndxNum, RowCurrIndxNum,
RowStart, RowEnd;
double SgnM;
double SrcTrngleEdgeLen[3], ObsTrngleEdgeLen[3], SrcTrngleNorm[3];
int SourceTrngleNode[3], SourceTrngleEdge[3], ObservTrngleNode[3],
ObservTrngleEdge[3];
for(IEdge=0;IEdge<=2;IEdge++)
{
double Buff[3];
ObservTrngleNode[IEdge] = TrngleNode[ObservTrngle][IEdge]-1;
ObservTrngleEdge[IEdge] = TrngleEdge[ObservTrngle][IEdge];
ObsTrngleEdgeLen[IEdge] =
EdgeLength[abs(ObservTrngleEdge[IEdge])-1];
VTXsub(TrngleCenTroid[ObservTrngle],
NodeCord[ObservTrngleNode[IEdge]], Buff);
RhoCtrd[ObservTrngle][IEdge][0] = Buff[0];
RhoCtrd[ObservTrngle][IEdge][1] = Buff[1];
RhoCtrd[ObservTrngle][IEdge][2] = Buff[2];
}
for(SourceTrngle=0;SourceTrngle<TotTrngleNum;SourceTrngle++)
{
double BE[3][3];
for(JEdge=0;JEdge<=2;JEdge++)
{
SourceTrngleNode[JEdge]=TrngleNode[SourceTrngle][JEdge]-1;
SourceTrngleEdge[JEdge]=TrngleEdge[SourceTrngle][JEdge];
SrcTrngleEdgeLen[JEdge]
=EdgeLength[abs(SourceTrngleEdge[JEdge])-1];
SrcTrngleNorm[JEdge]=TrngleNormal[SourceTrngle][JEdge];
}
for(QuadPointM=0;QuadPointM<=TotQuadPoint-1;QuadPointM++)
ComputeGaussQuadPoint( QuadPointM, ObservTrngleNode,
ObsPointArray[QuadPointM]);
for(QuadPointN=0;QuadPointN<=TotQuadPoint-1;QuadPointN++)
ComputeGaussQuadPoint( QuadPointN,SourceTrngleNode,
SrcPointArray[QuadPointN]);
ComputeCMatrix(ObservTrngle, SourceTrngle,
ObsTrngleEdgeLen, SrcTrngleEdgeLen,
ObservTrngleNode, SourceTrngleNode,
ObservTrngleEdge, SourceTrngleEdge,
ObsPointArray);
ComputeBddMatrix(ObservTrngle, SourceTrngle,
ObsTrngleEdgeLen, SrcTrngleEdgeLen,
ObservTrngleNode, SourceTrngleNode,
ObservTrngleEdge, SourceTrngleEdge,
SrcTrngleNorm, BE);
ComputeDMatrix(ObservTrngle, SourceTrngle,
ObsTrngleEdgeLen, SrcTrngleEdgeLen,
ObservTrngleNode, SourceTrngleNode,
ObservTrngleEdge, SourceTrngleEdge,
ObsPointArray, BE);
} /* end of SourceTrngle --- the big loop */
/****************Compute Eincident field Vector *************/
if( SourceType=='P' )
{
complex QE[3], EatCentroid[3];
double MidPoint[3];
MidPoint[0] = TrngleCenTroid[ObservTrngle][0];
MidPoint[1] = TrngleCenTroid[ObservTrngle][1];
MidPoint[2] = TrngleCenTroid[ObservTrngle][2];
ComputeSourceVector(WaveNumber,MidPoint,EatCentroid);
for(ObservEdge=0;ObservEdge<=2;ObservEdge++)
{
complex V1;
SgnM = Sign(ObservTrngleEdge[ObservEdge]);
V1 = COMplex_Add2(
Real_Mul(RhoCtrd[ObservTrngle][ObservEdge][0],
EatCentroid[0]),
Real_Mul(RhoCtrd[ObservTrngle][ObservEdge][1],
EatCentroid[1]),
Real_Mul(RhoCtrd[ObservTrngle][ObservEdge][2],
EatCentroid[2]));
QE[ObservEdge]= Real_Mul(
(ObsTrngleEdgeLen[ObservEdge]/2.0), V1);
QE[ObservEdge]=Real_Mul((SgnM),QE[ObservEdge]);
}
for(ObservEdge=0;ObservEdge<=2;ObservEdge++)
{
RowNum = abs(ObservTrngleEdge[ObservEdge]);
SgnM = Sign(ObservTrngleEdge[ObservEdge]);
MultRow= GlobalEdgeEnds[RowNum-1][2];
RowRootIndxNum= GlobalEdgeEnds[RowNum-1][3];
RowAddlIndxNum= GlobalEdgeEnds[RowNum-1][5];
RowCurrIndxNum= RowRootIndxNum+RowAddlIndxNum;
/* don't need to compute this edge. Thus jump to the next */
if(MultRow<=0) continue;
/* the following codes are for MultRow >0 */
if( SgnM >=0.0 )
{
for(IEdge=0;IEdge<PlusTrngleIndex[RowNum-1];IEdge++)
if ((PlusTrngleDetect[RowNum-1][IEdge]-1)==ObservTrngle)
{
RowStart = RowCurrIndxNum+IEdge;
/* revise on Jan 27 */
/*RowEnd = RowStart + MultRow;*/
RowEnd = RowStart +1;
break;
}
else continue;
}
else
{
for(IEdge=0;IEdge<MinusTrngleIndex[RowNum-1];IEdge++)
if ((MinusTrngleDetect[RowNum-1][IEdge]-1)
==ObservTrngle )
{
RowStart = RowCurrIndxNum;
RowEnd = RowStart + 1; break;
}
else continue;
}
if( MultRow>0 )
for(IEdge=RowStart;IEdge<RowEnd;IEdge++)
{
Einc[IEdge] = COMplex_Add(Einc[IEdge],QE[ObservEdge]);
Einc[IEdge] = Real_Mul(PlaneWaveE_Mag, Einc[IEdge] );
} /* for(IEdge */
} /* for ObservEdge */
} /* if(!strcmp) */
} /* for(ObservTrngle=0 */
end_time=time(0);
fprintf(OutF1, "Time Usage Report:\n\n");
fprintf(OutF1, "\tcomputing Bdd, C, D matrices\n");
fprintf(OutF1, "\t\ttime used: %d sec\n\n", (end_time-mid_time));
mid_time=end_time;
VcRowNum = INT_Vector(VSourceNum);
VdVector = CMPLX_Vector(HybrdBoundEdgeNum);
GdVector = CMPLX_Vector(HybrdBoundEdgeNum);
UdVector = CMPLX_Vector(HybrdBoundEdgeNum);
MdVector = CMPLX_Vector(HybrdBoundEdgeNum);
McVector = CMPLX_Vector(TotExtMetalEdgeNum);
VcVector = CMPLX_Vector(TotExtMetalEdgeNum);
Null_CMPLXVector(VdVector,HybrdBoundEdgeNum);
Null_CMPLXVector(GdVector,HybrdBoundEdgeNum);
Null_CMPLXVector(UdVector,HybrdBoundEdgeNum);
Null_CMPLXVector(MdVector,HybrdBoundEdgeNum);
Null_CMPLXVector(McVector,TotExtMetalEdgeNum);
Null_CMPLXVector(VcVector,TotExtMetalEdgeNum);
if( SourceType=='P' )
{
for(IEdge=0;IEdge<HybrdBoundEdgeNum;IEdge++)
VdVector[IEdge] = Einc[IEdge];
for(IEdge=HybrdBoundEdgeNum;IEdge<TotBoundEdgeNum;IEdge++)
VcVector[IEdge-HybrdBoundEdgeNum] = Einc[IEdge];
}
if( SourceType=='V' )
{
int flag=0;
for(IEdge=0;IEdge<VSourceNum;IEdge++)
{
for(JEdge=0;JEdge<TotBoundEdgeNum;JEdge++)
{
if( BoundEdgeStat[JEdge]==VSourceEdge[IEdge] )
{
/* vsource can not be defined on the hybridboudary */
if( JEdge<HybrdBoundEdgeNum )
{
printf("\nvsource can not be defined \
on the hybrid boudary\n Sorry, exit!!");
exit(1);
} /* if (JEdge */
VcRowNum[IEdge]=JEdge-HybrdBoundEdgeNum;
flag=1;
break;
} /* if( BoundEdgeStat */
} /* for(JEdge) */
/* Can't find the Vsource edge */
if( flag==0 )
{
printf("Sorry: Check your Vsource edge\n");
exit(1);
} /* if(flag) */
} /* for(IEdge) */
for(IEdge=0;IEdge<VSourceNum;IEdge++)
VcVector[VcRowNum[IEdge]]= COMplex_Cmplx(
VSourceMag[IEdge]*EdgeLength[VSourceEdge[IEdge]-1],0.0);
} /* if( !strcmp */
ComputeCorrectionTerm();
end_time=time(0);
fprintf(OutF1, "\tcomputing correction term\n");
fprintf(OutF1, "\t\ttime used: %d sec\n", (end_time-mid_time));
mid_time=end_time;
/**************************************************************/
/********* FINITE ELEMENT PART OF THE HYBRID CODE *************/
/**************************************************************/
GBLmatColIndex=INT_Vector(TotEdgeNum);
GBLmat= CMPLX_Matrix(TotEdgeNum,MaxAmatElementNum);
GBLmatColAddr=INT_Matrix(TotEdgeNum,MaxAmatElementNum);
Null_INTVector(GBLmatColIndex,TotEdgeNum);
Null_CMPLXMatrix(GBLmat,TotEdgeNum,MaxAmatElementNum);
for(IEdge=0;IEdge<TotEdgeNum;IEdge++)
for(JEdge=0;JEdge<MaxAmatElementNum;JEdge++)
GBLmatColAddr[IEdge][JEdge] = -1;
FemMatrixCompute(TotTetElement, MaxAmatElementNum, TetEdge,
GlobalEdgeEnds, NodeCord, TParam, Epsilon, GBLmatColAddr,
GBLmat, GBLmatColIndex,TetVolume);
/* consider the resistor */
/* the general form is L^2/ZL */
{
int i;
for(i=0; i<ResistorEdgeNum; i++)
{
int j,flag, RowIndex, ColIndex;
/* find the row index to the FEM matrix of the resistive edge */
flag=0;
for(j=0;j<TotInnerEdgeNum;j++)
if( InnerEdgeStat[j]==ResistorEdgeStat[i] )
{
flag=1;
RowIndex=j; /* row index to the FEM matrix */
break;
}
if( !flag )
{ fprintf(stderr, "The resistive edge can be only \
defined within the FEM region\n");
exit(1);
}
/* the column index */
ColIndex=SearchNonZeroElement(RowIndex, RowIndex);
GBLmat[RowIndex][ColIndex]=COMplex_Add(GBLmat[RowIndex][ColIndex],
COMplex_Cmplx(
pow( EdgeLength[ResistorEdgeStat[i]-1],2.0) /ResistorEdgeValue[i], 0.0) );
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -