📄 emap4.c
字号:
a_PML = a_PML_vec[HexNum];
b_PML = b_PML_vec[HexNum];
c_PML = c_PML_vec[HexNum];
Volume = TetVolume; Mult1 = (1296.0*Volume*Volume*Volume);
for(i=0;i<=5;i++) for(j=0;j<=5;j++) {
Length_i = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[i][0]]-1],
NodeCord[TetGlobalNodeNum[t][edge_end[i][1]]-1]);
Length_j = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[j][0]]-1],
NodeCord[TetGlobalNodeNum[t][edge_end[j][1]]-1]);
/* PML calculations (a_PML, b_PML and c_PML are the
diagonal entries in the dielectric matrix).
This formula comes from a matrix-vector dot-product */
S_mat[i][j] = COMplex_Add2(
Real_Mul(
( CoFactor[2][edge_end[i][0]] * CoFactor[3][edge_end[i][1]]
- CoFactor[3][edge_end[i][0]] * CoFactor[2][edge_end[i][1]] )
* ( CoFactor[2][edge_end[j][0]] * CoFactor[3][edge_end[j][1]]
- CoFactor[3][edge_end[j][0]] * CoFactor[2][edge_end[j][1]] ),
COMplex_Div(COMplex_Cmplx(1.0,0.0),a_PML))
, Real_Mul(
( CoFactor[3][edge_end[i][0]] * CoFactor[1][edge_end[i][1]]
- CoFactor[1][edge_end[i][0]] * CoFactor[3][edge_end[i][1]] )
* ( CoFactor[3][edge_end[j][0]] * CoFactor[1][edge_end[j][1]]
- CoFactor[1][edge_end[j][0]] * CoFactor[3][edge_end[j][1]] ),
COMplex_Div(COMplex_Cmplx(1.0,0.0),b_PML))
, Real_Mul(
( CoFactor[1][edge_end[i][0]] * CoFactor[2][edge_end[i][1]]
- CoFactor[2][edge_end[i][0]] * CoFactor[1][edge_end[i][1]] )
* ( CoFactor[1][edge_end[j][0]] * CoFactor[2][edge_end[j][1]]
- CoFactor[2][edge_end[j][0]] * CoFactor[1][edge_end[j][1]] ),
COMplex_Div(COMplex_Cmplx(1.0,0.0),c_PML)));
S_mat[i][j] = Real_Mul((Length_i*Length_j/Mult1),S_mat[i][j]) ;
}/*
printf(" %f %f\n",S_mat[0][0].x,S_mat[0][0].y);
fflush(stdout);*/
}
}
}
/***************************************************************************
* FUNCTION : T_Matrix()
****************************************************************************/
void T_Matrix()
/* This calculate the Fij part of the tetrahedron matrix
See reference [9]*/
{
int i,j;
double Mult1,Length_i,Length_j,T_mat[6][6];
Mult1 = 1.0/(720.0 * TetVolume);
if (Is_PML[HexNum] == 0) /* Not a PML layer */
{
/* (see reference [9] for explanation of the following terms */
T_mat[0][0] = ff(1,1) + ff(2,2) - ff(1,2);
T_mat[0][1] = 2 * ff(2,3) - ff(2,1) - ff(1,3) + ff(1,1);
T_mat[0][2] = 2 * ff(2,4) - ff(2,1) - ff(1,4) + ff(1,1);
T_mat[0][3] =-2 * ff(1,3) - ff(2,2) + ff(2,3) + ff(1,2);
T_mat[0][4] = 2 * ff(1,4) - ff(2,4) - ff(1,2) + ff(2,2);
T_mat[0][5] = ff(2,4) - ff(2,3) - ff(1,4) + ff(1,3);
T_mat[1][1] = ff(1,1) + ff(3,3) - ff(1,3);
T_mat[1][2] = 2 * ff(3,4) - ff(1,3) - ff(1,4) + ff(1,1);
T_mat[1][3] = 2 * ff(1,2) - ff(2,3) - ff(1,3) + ff(3,3);
T_mat[1][4] = ff(2,3) - ff(3,4) - ff(1,2) + ff(1,4);
T_mat[1][5] =-2 * ff(1,4) - ff(3,3) + ff(1,3) + ff(3,4);
T_mat[2][2] = ff(1,1) + ff(4,4) - ff(1,4);
T_mat[2][3] = ff(3,4) - ff(2,4) - ff(1,3) + ff(1,2);
T_mat[2][4] =-2 * ff(1,2) - ff(4,4) + ff(1,4) + ff(2,4);
T_mat[2][5] = 2 * ff(1,3) - ff(3,4) - ff(1,4) + ff(4,4);
T_mat[3][3] = ff(3,3) + ff(2,2) - ff(2,3);
T_mat[3][4] =-2 * ff(3,4) - ff(2,2) + ff(2,3) + ff(2,4);
T_mat[3][5] =-2 * ff(2,4) - ff(3,3) + ff(2,3) + ff(3,4);
T_mat[4][4] = ff(2,2) + ff(4,4) - ff(2,4);
T_mat[4][5] =-2 * ff(2,3) - ff(4,4) + ff(2,4) + ff(3,4);
T_mat[5][5] = ff(3,3) + ff(4,4) - ff(3,4);
for(i=0;i<=5;i++) for(j=0;j<=i-1;j++)
{T_mat[i][j] = T_mat[j][i];}
for(i=0;i<=5;i++) for(j=0;j<=5;j++) {
Length_i = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[i][0]]-1],
NodeCord[TetGlobalNodeNum[t][edge_end[i][1]]-1]);
Length_j = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[j][0]]-1],
NodeCord[TetGlobalNodeNum[t][edge_end[j][1]]-1]);
if(i==j) T_mat[i][j] = 2 * Length_i * Length_j * Mult1 * T_mat[i][j];
else T_mat[i][j] = Length_i * Length_j * Mult1 * T_mat[i][j];
T_mat[i][j] = WaveNumber*WaveNumber * T_mat[i][j] ;
T_mat1[i][j] = Real_Mul(T_mat[i][j],RELPerm[HexNum]);
}
}
else /* Part of PML layer */
/* For PMLs, the matrix calculation is different because the
dielectric matrix is anisotropic. The anisotropicity is used in the
ff1 calculations. See function above the main routine.
The calculation is same as above but uses complex numbers
and includes the PML dielectric matrix */
{
a_PML = a_PML_vec[HexNum];
b_PML = b_PML_vec[HexNum];
c_PML = c_PML_vec[HexNum];
T_mat1[0][0] = COMplex_Add(ff1(1,1),COMplex_Sub(ff1(2,2),ff1(1,2)));
T_mat1[1][1] = COMplex_Add(ff1(1,1),COMplex_Sub(ff1(3,3),ff1(1,3)));
T_mat1[2][2] = COMplex_Add(ff1(1,1),COMplex_Sub(ff1(4,4),ff1(1,4)));
T_mat1[3][3] = COMplex_Add(ff1(3,3),COMplex_Sub(ff1(2,2),ff1(2,3)));
T_mat1[4][4] = COMplex_Add(ff1(2,2),COMplex_Sub(ff1(4,4),ff1(2,4)));
T_mat1[5][5] = COMplex_Add(ff1(3,3),COMplex_Sub(ff1(4,4),ff1(3,4)));
T_mat1[0][5] = COMplex_Sub(COMplex_Sub(ff1(2,4),ff1(2,3)),
COMplex_Sub(ff1(1,4),ff1(1,3)));
T_mat1[1][4] = COMplex_Sub(COMplex_Sub(ff1(2,3),ff1(3,4)),
COMplex_Sub(ff1(1,2),ff1(1,4)));
T_mat1[2][3] = COMplex_Sub(COMplex_Sub(ff1(3,4),ff1(2,4)),
COMplex_Sub(ff1(1,3),ff1(1,2)));
T_mat1[0][1] = COMplex_Sub(COMplex_Sub(
COMplex_Add(ff1(2,3),ff1(2,3)),ff1(2,1)),
COMplex_Sub(ff1(1,3),ff1(1,1)));
T_mat1[0][2] = COMplex_Sub(COMplex_Sub(
COMplex_Add(ff1(2,4),ff1(2,4)),ff1(2,1)),
COMplex_Sub(ff1(1,4),ff1(1,1)));
T_mat1[0][4] = COMplex_Sub(COMplex_Sub(
COMplex_Add(ff1(1,4),ff1(1,4)),ff1(2,4)),
COMplex_Sub(ff1(1,2),ff1(2,2)));
T_mat1[1][2] = COMplex_Sub(COMplex_Sub(
COMplex_Add(ff1(3,4),ff1(3,4)),ff1(1,3)),
COMplex_Sub(ff1(1,4),ff1(1,1)));
T_mat1[1][3] = COMplex_Sub(COMplex_Sub(
COMplex_Add(ff1(1,2),ff1(1,2)),ff1(2,3)),
COMplex_Sub(ff1(1,3),ff1(3,3)));
T_mat1[2][5] = COMplex_Sub(COMplex_Sub(
COMplex_Add(ff1(1,3),ff1(1,3)),ff1(3,4)),
COMplex_Sub(ff1(1,4),ff1(4,4)));
T_mat1[0][3] = COMplex_Sub(COMplex_Add(ff1(2,3),ff1(1,2)),
COMplex_Add2(ff1(1,3),ff1(1,3),ff1(2,2)));
T_mat1[1][5] = COMplex_Sub(COMplex_Add(ff1(1,3),ff1(3,4)),
COMplex_Add2(ff1(1,4),ff1(1,4),ff1(3,3)));
T_mat1[2][4] = COMplex_Sub(COMplex_Add(ff1(1,4),ff1(2,4)),
COMplex_Add2(ff1(1,2),ff1(1,2),ff1(4,4)));
T_mat1[3][4] = COMplex_Sub(COMplex_Add(ff1(2,3),ff1(2,4)),
COMplex_Add2(ff1(3,4),ff1(3,4),ff1(2,2)));
T_mat1[3][5] = COMplex_Sub(COMplex_Add(ff1(2,3),ff1(3,4)),
COMplex_Add2(ff1(2,4),ff1(2,4),ff1(3,3)));
T_mat1[4][5] = COMplex_Sub(COMplex_Add(ff1(2,4),ff1(3,4)),
COMplex_Add2(ff1(2,3),ff1(2,3),ff1(4,4)));
for(i=0;i<=5;i++) for(j=0;j<=i-1;j++)
{T_mat1[i][j] = T_mat1[j][i];}
for(i=0;i<=5;i++) for(j=0;j<=5;j++) {
Length_i = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[i][0]]-1],
NodeCord[TetGlobalNodeNum[t][edge_end[i][1]]-1]);
Length_j = VTXmag(NodeCord[TetGlobalNodeNum[t][edge_end[j][0]]-1],
NodeCord[TetGlobalNodeNum[t][edge_end[j][1]]-1]);
if(i==j) T_mat1[i][j] = Real_Mul(2*Length_i*Length_j*Mult1,T_mat1[i][j]);
else T_mat1[i][j] = Real_Mul(Length_i*Length_j*Mult1,T_mat1[i][j]);
T_mat1[i][j] = Real_Mul(WaveNumber*WaveNumber,T_mat1[i][j]) ;
}
}
}
/***************************************************************************
* FUNCTION : TetraSubMatrix()
****************************************************************************/
void TetraSubMatrix()
/* Combine the Eij and Fij matrices into the A matrix */
{
int i,j;
for(i=0;i<=5;i++) for(j=0;j<=5;j++)
{A_mat[t][i][j]=COMplex_Null();}
for(i=0;i<=5;i++) for(j=0;j<=5;j++)
{A_mat[t][i][j] = COMplex_Sub(Real_Mul(4.0,S_mat[i][j]),T_mat1[i][j]);}
for(i=0;i<=5;i++) for(j=0;j<=5;j++)
{A_mat[t][i][j] = Real_Mul(Sign(TetEdgeNum[t][j]),A_mat[t][i][j]);}
for(i=0;i<=5;i++) for(j=0;j<=5;j++)
{A_mat[t][i][j] = Real_Mul(Sign(TetEdgeNum[t][i]),A_mat[t][i][j]);}
}
/***************************************************************************
* FUNCTION : FindGlobalMatrix()
****************************************************************************/
void FindGlobalMatrix()
/* Assembles the tetrahedron matrices into the global matrices
and stores the global matrix in the row-indexed format */
{
int m,n,k,l,rval;
{
for(m=0;m<=5;m++) for(n=0;n<=5;n++) {
k=absol(TetEdgeNum[t][m])-1; l=absol(TetEdgeNum[t][n])-1;
if(COMplex_Abs(A_mat[t][m][n])>1.0E-10) {
rval=Srch_NZ_Element_Column(k,l); /* Search if there is such an element */
/* Matrix term already exists, so add the new term to the old one */
if(rval>=0) {
GBL_Matrix_Data[k][rval]=COMplex_Add(
GBL_Matrix_Data[k][rval],A_mat[t][m][n]);
}
/* Create new global matrix element */
else {
GBL_Matrix_Data[k][GBL_Matrix_Index[k]]=A_mat[t][m][n];
GBL_Matrix_ColNos[k][GBL_Matrix_Index[k]] = l;
GBL_Matrix_Index[k]++;
}
}
}
}
}
/***************************************************************************
* FUNCTION : PartitionGlobalMatrix()
****************************************************************************/
void PartitionGlobalMatrix()
{
int loop_i,loop_j,No_cols,k_row=0,II,JJ,flag=0;
/* Just for debugging
for (II=0;II<TotEdgeNum;II++)
{
for(JJ=0;JJ<GBL_Matrix_Index[II];JJ++)
{
if ( GBL_Matrix_Data[II][JJ].y > 1000)
{
printf("\n %d \n", II);
printf("%f ",NodeCord[TetGlobalEdgeEnds[II][0]][0]);
printf("%f ",NodeCord[TetGlobalEdgeEnds[II][0]][1]);
printf("%f ",NodeCord[TetGlobalEdgeEnds[II][0]][2]);
printf("%f ",NodeCord[TetGlobalEdgeEnds[II][1]][0]);
printf("%f ",NodeCord[TetGlobalEdgeEnds[II][1]][1]);
printf("%f\n",NodeCord[TetGlobalEdgeEnds[II][1]][2]);
}
}
}
*/
for (loop_i=0;loop_i<TotEdgeNum;loop_i++) /* Loop over all the edges */
{
flag=0;
if (ForceStat[loop_i] == 0) /* Inner edge */
{
No_cols = GBL_Matrix_Index[loop_i];
/* printf("%d \n",No_cols); */
for (loop_j=0;loop_j<No_cols;loop_j++) /*Loop over all non-zero elements
in row loop_i */
{
if (ForceStat[GBL_Matrix_ColNos[loop_i][loop_j]] == 0)
/* Inner edge */
{
/* Construct the LHS matrix from the elements
of the global matrix */
LHS_Matrix_Data[k_row][LHS_Matrix_Index[k_row]]=
GBL_Matrix_Data[loop_i][loop_j];
if ((loop_i == GBL_Matrix_ColNos[loop_i][loop_j])
&& (COMplex_Abs(Resist[loop_i])>TOLERANCE))
{
/* Add the term for the resistor */
/* printf("Inside this resistor loop"); */
/* fflush(stdout); */
LHS_Matrix_Data[k_row][LHS_Matrix_Index[k_row]] =
COMplex_Add(LHS_Matrix_Data[k_row][LHS_Matrix_Index[k_row]],
Resist[loop_i]);
}
LHS_Matrix_ColNos[k_row][LHS_Matrix_Index[k_row]]=
InnerEdgeStat1[GBL_Matrix_ColNos[loop_i][loop_j]];
LHS_Matrix_Index[k_row]++;
}
else
{
/* Construct the RHS Vector */
RHSVector[k_row]=COMplex_Add(RHSVector[k_row],COMplex_Mul(
ForcdValue[ForcdEdgeStat1[GBL_Matrix_ColNos[loop_i][loop_j]]],
GBL_Matrix_Data[loop_i][loop_j]));
if ((COMplex_Abs(Isource[loop_i])>TOLERANCE) && (flag==0))
{
/* Add "isource" term to the RHS vector */
flag=1; /* Add the "isource" term only once */
/*printf("Inside this isource loop %d %d\n",loop_i,No_cols);*/
RHSVector[k_row]=COMplex_Add(RHSVector[k_row],Isource[loop_i]);
}
}/* end of else*/
} /* end of for loop_j*/
fflush(stdout);
k_row++;
} /* end of if */
} /* end of for loop_i*/
printf("%d %d\n",TotInnerEdgeNum,k_row);
}
/*************************************************************************
* SUBROUTINE : Conjugatesolver()
*************************************************************************/
/**********************************************
FUNCTION :Vnorm()
Computes the Euclidean norm of a vector: |*|^2
***********************************************/
double Vnorm(Vec,Size)
int Size; complex *Vec;
{
int II; double sum=0;
for(II=0;II<Size;II++)
{
sum += COMplex_Abs(Vec[II])*COMplex_Abs(Vec[II]);
}
return sum;
}
/**********************************************
FUNCTION :Inner_Prod()
Computes the Inner product of two vectors : <A*,B>
***********************************************/
complex Inner_Prod(AVec,BVec,Size)
int Size; complex *AVec,*BVec;
{
int II; complex I_P;
I_P = COMplex_Null();
for(II=0;II<Size;II++)
{
I_P = COMplex_Add(I_P,COMplex_Mul(COMplex_Conjg(AVec[II]),BVec[II]));
}
return I_P;
}
/**********************************************
FUNCTION :MaxVecProd()
Multiply a complex Matrix with a vector
***********************************************/
MatrixVectorProd(S,Size,XVec,AVec,RIIndex,RICol,RIDat)
char S;
int Size,**RICol,*RIIndex;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -