📄 emap4.c
字号:
TetGlobalNodeNum[3][2]=HexNode[HexNum][4];
TetGlobalNodeNum[2][3]=HexNode[HexNum][5];
TetGlobalNodeNum[3][3]=HexNode[HexNum][5];
TetGlobalNodeNum[4][0]=HexNode[HexNum][0];
TetGlobalNodeNum[4][2]=HexNode[HexNum][3];
TetGlobalNodeNum[4][1]=HexNode[HexNum][5];
TetGlobalNodeNum[4][3]=HexNode[HexNum][6];
}
/***************************************************************************
* FUNCTION : HexHedraSubDiv2()
****************************************************************************/
void HexHedraSubDiv2()
/* Divides the 8 nodes in the hexahedron into five groups of 4 each
in another way*/
{
TetGlobalNodeNum[0][0]=HexNode[HexNum][0];
TetGlobalNodeNum[1][0]=HexNode[HexNum][1];
TetGlobalNodeNum[0][1]=HexNode[HexNum][1];
TetGlobalNodeNum[1][1]=HexNode[HexNum][7];
TetGlobalNodeNum[0][2]=HexNode[HexNum][2];
TetGlobalNodeNum[1][2]=HexNode[HexNum][4];
TetGlobalNodeNum[0][3]=HexNode[HexNum][4];
TetGlobalNodeNum[1][3]=HexNode[HexNum][5];
TetGlobalNodeNum[2][0]=HexNode[HexNum][1];
TetGlobalNodeNum[3][0]=HexNode[HexNum][2];
TetGlobalNodeNum[2][1]=HexNode[HexNum][7];
TetGlobalNodeNum[3][1]=HexNode[HexNum][6];
TetGlobalNodeNum[2][2]=HexNode[HexNum][3];
TetGlobalNodeNum[3][2]=HexNode[HexNum][4];
TetGlobalNodeNum[2][3]=HexNode[HexNum][2];
TetGlobalNodeNum[3][3]=HexNode[HexNum][7];
TetGlobalNodeNum[4][0]=HexNode[HexNum][1];
TetGlobalNodeNum[4][2]=HexNode[HexNum][4];
TetGlobalNodeNum[4][1]=HexNode[HexNum][2];
TetGlobalNodeNum[4][3]=HexNode[HexNum][7];
}
/***************************************************************************
* FUNCTION : EdgeSubDiv1()
****************************************************************************/
void EdgeSubDiv1()
/* Divides the 18 edges in the hexahedron into five groups of 6 each */
{
TetEdgeNum[0][0] = HexEdgeNum[HexNum][2];
TetEdgeNum[0][1] = HexEdgeNum[HexNum][1];
TetEdgeNum[0][2] = HexEdgeNum[HexNum][8];
TetEdgeNum[0][3] =-HexEdgeNum[HexNum][4];
TetEdgeNum[0][4] =-HexEdgeNum[HexNum][11];
TetEdgeNum[0][5] = HexEdgeNum[HexNum][10];
TetEdgeNum[1][0] = HexEdgeNum[HexNum][9];
TetEdgeNum[1][1] = HexEdgeNum[HexNum][12];
TetEdgeNum[1][2] = HexEdgeNum[HexNum][11];
TetEdgeNum[1][3] = HexEdgeNum[HexNum][16];
TetEdgeNum[1][4] =-HexEdgeNum[HexNum][15];
TetEdgeNum[1][5] =-HexEdgeNum[HexNum][17];
TetEdgeNum[2][0] = HexEdgeNum[HexNum][0];
TetEdgeNum[2][1] = HexEdgeNum[HexNum][2];
TetEdgeNum[2][2] = HexEdgeNum[HexNum][6];
TetEdgeNum[2][3] = HexEdgeNum[HexNum][3];
TetEdgeNum[2][4] =-HexEdgeNum[HexNum][7];
TetEdgeNum[2][5] = HexEdgeNum[HexNum][9];
TetEdgeNum[3][0] = HexEdgeNum[HexNum][8];
TetEdgeNum[3][1] = HexEdgeNum[HexNum][5];
TetEdgeNum[3][2] = HexEdgeNum[HexNum][6];
TetEdgeNum[3][3] =-HexEdgeNum[HexNum][14];
TetEdgeNum[3][4] = HexEdgeNum[HexNum][15];
TetEdgeNum[3][5] = HexEdgeNum[HexNum][13];
TetEdgeNum[4][0] = HexEdgeNum[HexNum][6];
TetEdgeNum[4][1] = HexEdgeNum[HexNum][2];
TetEdgeNum[4][2] = HexEdgeNum[HexNum][8];
TetEdgeNum[4][3] =-HexEdgeNum[HexNum][9];
TetEdgeNum[4][4] =-HexEdgeNum[HexNum][15];
TetEdgeNum[4][5] = HexEdgeNum[HexNum][11];
}
/***************************************************************************
* FUNCTION : EdgeSubDiv2()
****************************************************************************/
void EdgeSubDiv2()
/* Divides the 18 edges in the hexahedron into five groups of 6 each
in another way */
{
TetEdgeNum[0][0] = HexEdgeNum[HexNum][0];
TetEdgeNum[0][1] = HexEdgeNum[HexNum][1];
TetEdgeNum[0][2] = HexEdgeNum[HexNum][5];
TetEdgeNum[0][3] = HexEdgeNum[HexNum][2];
TetEdgeNum[0][4] =-HexEdgeNum[HexNum][6];
TetEdgeNum[0][5] = HexEdgeNum[HexNum][8];
TetEdgeNum[1][0] = HexEdgeNum[HexNum][9];
TetEdgeNum[1][1] = HexEdgeNum[HexNum][6];
TetEdgeNum[1][2] = HexEdgeNum[HexNum][7];
TetEdgeNum[1][3] =-HexEdgeNum[HexNum][15];
TetEdgeNum[1][4] = HexEdgeNum[HexNum][16];
TetEdgeNum[1][5] = HexEdgeNum[HexNum][13];
TetEdgeNum[2][0] = HexEdgeNum[HexNum][9];
TetEdgeNum[2][1] = HexEdgeNum[HexNum][3];
TetEdgeNum[2][2] = HexEdgeNum[HexNum][2];
TetEdgeNum[2][3] =-HexEdgeNum[HexNum][12];
TetEdgeNum[2][4] = HexEdgeNum[HexNum][11];
TetEdgeNum[2][5] =-HexEdgeNum[HexNum][4];
TetEdgeNum[3][0] = HexEdgeNum[HexNum][10];
TetEdgeNum[3][1] = HexEdgeNum[HexNum][8];
TetEdgeNum[3][2] = HexEdgeNum[HexNum][11];
TetEdgeNum[3][3] =-HexEdgeNum[HexNum][14];
TetEdgeNum[3][4] =-HexEdgeNum[HexNum][17];
TetEdgeNum[3][5] = HexEdgeNum[HexNum][15];
TetEdgeNum[4][0] = HexEdgeNum[HexNum][2];
TetEdgeNum[4][1] = HexEdgeNum[HexNum][6];
TetEdgeNum[4][2] = HexEdgeNum[HexNum][9];
TetEdgeNum[4][3] = HexEdgeNum[HexNum][8];
TetEdgeNum[4][4] =-HexEdgeNum[HexNum][11];
TetEdgeNum[4][5] = HexEdgeNum[HexNum][15];
}
/***************************************************************************
* FUNCTION : GlobalEdgeEndsDiv1()
****************************************************************************/
void GlobalEdgeEndsDiv1()
/* Specifies the end-nodes for each edge for the first division
See the "Division of hexahedra into five tetrahedra" figure*/
{
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][1]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][0]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][1]=HexNode[HexNum][7];
}
/***************************************************************************
* FUNCTION : GlobalEdgeEndsDiv2()
****************************************************************************/
void GlobalEdgeEndsDiv2()
/* Specifies the end-nodes for each edge for the second division */
{
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][0]-1][1]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][1]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][2]-1][1]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][3]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][4]-1][1]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][0]=HexNode[HexNum][0];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][5]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][6]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][7]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][8]-1][1]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][0]=HexNode[HexNum][1];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][9]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][10]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][0]=HexNode[HexNum][2];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][11]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][0]=HexNode[HexNum][3];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][12]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][13]-1][1]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][14]-1][1]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][0]=HexNode[HexNum][4];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][15]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][0]=HexNode[HexNum][5];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][16]-1][1]=HexNode[HexNum][7];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][0]=HexNode[HexNum][6];
TetGlobalEdgeEnds[HexEdgeNum[HexNum][17]-1][1]=HexNode[HexNum][7];
}
/***************************************************************************
* FUNCTION : Find_Volume_Determinent()
****************************************************************************/
void FindTetraHedronVolume()
/* Finds the volume of each tetrahedron to be used further in calculations */
{
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]=NodeCord[TetGlobalNodeNum[t][Count_i]-1][Count_j];
VTXsub1(1,0,Buff,Buff1);
VTXsub1(2,0,Buff,Buff2);
VTXsub1(3,0,Buff,Buff3);
VTXcross(Buff1,Buff2,Buff4);
TetVolume=(1.0/6.0)*(Buff3[0]*Buff4[0]+Buff3[1]*Buff4[1]+Buff3[2]*Buff4[2]);
}
/***************************************************************************
* FUNCTION : ComputeTetraCoFactor()
****************************************************************************/
void ComputeTetraCoFactor()
/* This calculates the a's, b's, c's and d's in reference [9] */
/* This just involves some basic matrix simulations */
{
int i,j,RowNum,ColNum,Row,Col,Count=0;
double CoF[4][4],Buff[4][4],CoFs[4][4];
{
for(i=0;i<=3;i++) for(j=0;j<=3;j++) {
if(j==3) Buff[i][j]= 1.0;
else Buff[i][j]= NodeCord[TetGlobalNodeNum[t][i]-1][j];
}
/* The main matrix is made up of a column of 1's and
the other columns are made of the node coordinates */
/* The co-factor matrices are formed from the main matrix
by taking part of the matrix */
for(RowNum=0;RowNum<=3;RowNum++)
for(ColNum=0;ColNum<=3;ColNum++) {
Row=0;Col=0;Count=0;
for(i=0;i<=3;i++)
for(j=0;j<=3;j++) {
if((i!=RowNum) && (j!=ColNum)) {
Count++;
if(Count<4) Row=0;
if((Count>3)&&(Count<7)) Row=1;
else if(Count>6) Row=2;
CoF[Row][Col]=Buff[i][j];
Col++; Col=Col%3; }
}
/* The next few steps calculate the co-factors
by taking the determinant of the sub-matrix */
Row=0;Col=0;Count=0;
CoFs[RowNum][ColNum]
=CoF[0][0]*((CoF[1][1]*CoF[2][2])-(CoF[1][2]*CoF[2][1]))
-CoF[0][1]*((CoF[1][0]*CoF[2][2])-(CoF[2][0]*CoF[1][2]))
+CoF[0][2]*((CoF[1][0]*CoF[2][1])-(CoF[2][0]*CoF[1][1]));
if((RowNum+ColNum)%2!=0) CoFs[RowNum][ColNum]*=-1.0;
}
for(RowNum=0;RowNum<=3;RowNum++) {
if (RowNum==0) Col = 3;
else if(RowNum==1) Col = 0;
else if(RowNum==2) Col = 1;
else if(RowNum==3) Col = 2;
for(ColNum=0;ColNum<=3;ColNum++) {
CoFactor[RowNum][ColNum] = -CoFs[ColNum][Col];
}}
}
}
/***************************************************************************
* FUNCTION : S_Matrix()
****************************************************************************/
void S_Matrix()
/* This calculates the Eij part of the tetrahedron matrix
All the formulas used are in reference [9] or see reference [1]*/
{
int i,j; double Length_i,Length_j,Mult1,Volume;
double S_mat1[6][6];
{
if (Is_PML[HexNum] == 0) /* Not a PML */
{
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]);
S_mat1[i][j] = (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]] )
+ ( 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]] )
+ ( 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]] );
S_mat1[i][j] = (Length_i * Length_j * S_mat1[i][j])/ Mult1 ;
S_mat[i][j] = COMplex_Cmplx(S_mat1[i][j],0.0); /* Make it complex */
}/*
printf(" %f %f\n",S_mat[0][0].x,S_mat[0][0].y);
*/ }
else /* The hexahedron is a PML */
/* For PMLs, the matrix calculation is different because the
dielectric matrix is anisotropic */
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -