📄 emap4.c
字号:
main(argc, argv)
/* MAIN PROGRAM */
int argc;
char *argv[];
{
int con,i,j,k;
double FreeSpaceVel,AbsPermeable,AbsPermitt,WaveLength;
/*Print error messages if the command line is not correct */
if (argc > 2) {
fprintf(stderr,"Usage: Emap4 [input_file] \n"); exit(1);}
if (argc < 2) { argv[1] = "emap4.in"; strcpy(Ifile, argv[1]); }
else strcpy(Ifile, argv[1]);
fprintf(stdout, "\nEMAP4\n\n");
fprintf(stdout, "The input will be read from the file: %s\n", Ifile);
InF = fopen(Ifile, "r");
if (InF == NULL) {
fprintf(stderr, " Error: Can't open input file.\n"); exit(1); }
/* Pre-processing begins */
Read_Input_Pass_1(); /* check for errors and inconsistencies
in the input file */
ParameterInfo(); /* Take in the values of many global parameters */
/* A little bit of initialization */
HexNode = INT_Matrix(TotNumHexHedron,8);
HexEdgeNum = INT_Matrix(TotNumHexHedron,19);
RELPerm = CMPLX_Vector(TotNumHexHedron,2);
NodeCord = FLOAT_Matrix(TotNodeNum,3);
TetGlobalNodeNum = INT_Matrix(5,4);
TetEdgeNum = INT_Matrix(5,6);
TetGlobalEdgeEnds = INT_Matrix(TotEdgeNum,2);
for(i=0;i<TotNumHexHedron;i++) {
for(j=0;j<8;j++) {HexNode[i][j] = 0;}
for(j=0;j<19;j++) {HexEdgeNum[i][j] = 0;}
RELPerm[i] = COMplex_Null();
}
for(i=0;i<TotNodeNum;i++) {
for(j=0;j<3;j++) {NodeCord[i][j] = 0.0;}
}
for(i=0;i<TotEdgeNum;i++) {
for(j=0;j<2;j++) {TetGlobalEdgeEnds[i][j] = 0;}
}
for(i=0;i<5;i++) {
for(j=0;j<4;j++) {TetGlobalNodeNum[i][j] = 0;}
for(j=0;j<6;j++) {TetEdgeNum[i][j] = 0;}
}
/* Mesh Generation begins*/
Read_Input_Pass_2(); /* Complete the reading of the input file and
complete the initialization */
AssignGlobalCoorDinates(); /* Assign global coordinates
based on the input file */
AssignHexHedronNodeNumbering(); /* Global node and edge numbering */
AssignHexHedronEdgeNumbering();
/* More initialization */
GBL_Matrix_Data=CMPLX_Matrix(TotEdgeNum,GBL_Matrix_MaxColNos);
GBL_Matrix_ColNos=INT_Matrix(TotEdgeNum,GBL_Matrix_MaxColNos);
GBL_Matrix_Index=INT_Vector(TotEdgeNum);
LHS_Matrix_Data=CMPLX_Matrix(TotInnerEdgeNum,GBL_Matrix_MaxColNos);
LHS_Matrix_ColNos=INT_Matrix(TotInnerEdgeNum,GBL_Matrix_MaxColNos);
LHS_Matrix_Index=INT_Vector(TotInnerEdgeNum);
for(freq_step=0;freq_step<NUM_OF_STEPS;freq_step++){
/* For Simulations over a range of frequencies
(All the routines above are done only the first step) */
for(i=0;i<TotEdgeNum;i++)
{
GBL_Matrix_Index[i] = 0;
for(j=0;j<GBL_Matrix_MaxColNos;j++)
{
GBL_Matrix_Data[i][j] = COMplex_Null();
/* Recreate matrix for every iteration */
GBL_Matrix_ColNos[i][j] = 0;
}
}
for(i=0;i<TotInnerEdgeNum;i++)
{
LHS_Matrix_Index[i] = 0;
for(j=0;j<GBL_Matrix_MaxColNos;j++)
{
LHS_Matrix_Data[i][j] = COMplex_Null();
LHS_Matrix_ColNos[i][j] = 0;
}
}
/* Step for every hexahedron
(XdiM, YdiM and ZdiM are the total number of bricks in each direction)*/
for(k=0;k<=ZdiM-1;k++)
for(j=0;j<=YdiM-1;j++)
for(i=0;i<=XdiM-1;i++)
{
HexNum = k*YdiM*XdiM + j*XdiM + i; /*Hexahedron number */
con = pow(-1.0,(double)i)*pow(-1.0,(double)j)*pow(-1.0,(double)k);
/* Divide hexahedron into 5 tetrahedra in two different ways */
/* Adjacent hexahedra are divided in different ways (using "con")*/
if(con==1) {HexHedraSubDiv1(); EdgeSubDiv1(); GlobalEdgeEndsDiv1();}
else {HexHedraSubDiv2(); EdgeSubDiv2(); GlobalEdgeEndsDiv2();}
/* Compute the tetrahedron matrices and
assemble them into the final matrix */
for(t=0;t<=4;t++)
{
FindTetraHedronVolume(); /* Finds volume (see user's guide) */
ComputeTetraCoFactor(); /* Computes a's, b's c's d's */
S_Matrix(); /* Computes the Eij part */
T_Matrix(); /* COmputes the Fij part */
TetraSubMatrix(); /* Complete Tetrahedron matrix */
FindGlobalMatrix();
/* Assemble tetrahedron matrices into global matrix */
}
}
if(freq_step==0) /* This comes frequently :
Call Malloc only when it is the first frequency step */
{
RHSVector = CMPLX_Vector(TotInnerEdgeNum);
Old_Solution = CMPLX_Vector(TotInnerEdgeNum);
for(i=0;i<TotInnerEdgeNum;i++)
{
Old_Solution[i]=COMplex_Null();
}
}
for(i=0;i<TotInnerEdgeNum;i++)
{
RHSVector[i]=COMplex_Null();
}
/* printf("Reached here \n"); */
/* fflush(stdout); */
/* Partition the global matrix into the "inner" and "forced" edges */
PartitionGlobalMatrix();
/* printf("Reached here \n"); */
/* fflush(stdout); */
/* Solve the final matrix equation */
ConjugateSolver(TotInnerEdgeNum,LHS_Matrix_Index,LHS_Matrix_ColNos,LHS_Matrix_Data,RHSVector,MAX_ITERATIONS,TOLERANCE);
/* printf("Reached here \n"); */
/* fflush(stdout); */
/* The following routine takes the "inner edge" vector along with the
"forced edge" vector and calculates and
writes out to all the required files */
Produce_Output();
OperateFreq +=FREQ_INC; /* Goto next frequency in the frequency range */
/* Update all the variables used */
AbsPermeable = 1.25663706144E-06;
AbsPermitt = 8.8542E-12;
FreeSpaceVel = 1.0/sqrt(AbsPermeable*AbsPermitt);
WaveLength = FreeSpaceVel/(OperateFreq*1.0E+06);
WaveNumber = 2.0*M_PI/WaveLength;
printf("OperateFrequency=%g MHz\n",OperateFreq);
printf("WaveLength=%f m\n",WaveLength );
printf("FreeSpaceVelocity=%E m/s\n",FreeSpaceVel );
printf("WaveNumber=%g\n",WaveNumber);}
return(0);
}
/***************************************************************************/
void ParameterInfo() /* This is done only for the first frequency step */
{
AbsPermeable = 1.25663706144E-06;
AbsPermitt = 8.8542E-12;
FreeSpaceVel = 1.0/sqrt(AbsPermeable*AbsPermitt);
WaveLength = FreeSpaceVel/(OperateFreq*1.0E+06);
WaveNumber = 2.0*M_PI/WaveLength;
Omega = WaveNumber*FreeSpaceVel;
printf("OperateFrequency=%g MHz\n",OperateFreq);
printf("WaveLength=%f m\n",WaveLength );
printf("FreeSpaceVelocity=%E m/s\n",FreeSpaceVel );
printf("WaveNumber=%g\n",WaveNumber);
/* Calculate the necessary parameters */
TotEdgeNum = 6*XdiM*YdiM*ZdiM + (XdiM+YdiM+ZdiM)
+ 3*(XdiM*YdiM+YdiM*ZdiM+ZdiM*XdiM);
TotNodeNum = (XdiM+1)*(YdiM+1)*(ZdiM+1);
TotNumHexHedron = XdiM*YdiM*ZdiM;
TotTetElement = 5 * TotNumHexHedron;
TotTrngleNum = 4*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM);
TotBoundEdgeNum = 6*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM);
TotBoundNodeNum = 2*(XdiM*YdiM + YdiM*ZdiM + ZdiM*XdiM+1);
printf("Total Edge Number = %d \n",TotEdgeNum);
printf("Total Node Number = %d \n",TotNodeNum);
printf("Total Cube Number = %d \n",TotNumHexHedron);
printf("Total Tetrahedron Number = %d \n",TotTetElement);
printf("Total Boundary Edge Number = %d \n",TotBoundEdgeNum);
printf("Total Boundary Node Number = %d \n",TotBoundNodeNum);
/* Local edge orientation of tetrahedron edges */
edge_end[0][0] = 0; edge_end[0][1] = 1;
edge_end[1][0] = 0; edge_end[1][1] = 2;
edge_end[2][0] = 0; edge_end[2][1] = 3;
edge_end[3][0] = 1; edge_end[3][1] = 2;
edge_end[4][0] = 3; edge_end[4][1] = 1;
edge_end[5][0] = 2; edge_end[5][1] = 3;
}
/***************************************************************************
* FUNCTION : AssignGlobalCoorDinates()
****************************************************************************/
void AssignGlobalCoorDinates()
/* Calculate the node coordinates for each global node */
{
int i,j,k,NodeNum; double Xtemp,Ytemp,Ztemp;
{
NodeNum=0; Xtemp=0.0; Ytemp=0.0; Ztemp=0.0;
for(k=0;k<=ZdiM;k++) {
for(j=0;j<=YdiM;j++) {
for(i=0;i<=XdiM;i++) {
NodeCord[NodeNum][0] = Xtemp;
NodeCord[NodeNum][1] = Ytemp;
NodeCord[NodeNum][2] = Ztemp;
/* printf(" \n %f %f %f",NodeCord[NodeNum][0],NodeCord[NodeNum][1],NodeCord[NodeNum][2]);*/
NodeNum++; Xtemp += DivisorX[i+1];
}
Xtemp = 0.0; Ytemp = Ytemp + DivisorY[j+1];
}
Ytemp = 0.0; Ztemp = Ztemp + DivisorZ[k+1];
}
}
}
/***************************************************************************
* FUNCTION : AssignHexHedronNodeNumbering()
****************************************************************************/
void AssignHexHedronNodeNumbering()
/* Assign global node numbering
(see "Mesh generation" chapter in the user's guide
to verify these formulas )*/
{
int i;
{
start=1;
for(i=0;i<=TotNumHexHedron-1;i++)
{
HexNode[i][0]=start;
HexNode[i][1]=HexNode[i][0]+1;
HexNode[i][2]=HexNode[i][0]+XdiM+1;
HexNode[i][3]=HexNode[i][2]+1;
HexNode[i][4]=HexNode[i][0]+(XdiM+1)*(YdiM+1);
HexNode[i][5]=HexNode[i][4]+1;
HexNode[i][6]=HexNode[i][4]+(XdiM+1);
HexNode[i][7]=HexNode[i][6]+1;
if(((i+1)%((YdiM*XdiM))==0)) {start=start+(XdiM+3);}
else {start=start+1;}
if((start%(XdiM+1))==0) { start++; }
}
}
}
/***************************************************************************
* FUNCTION : AssignHexHedronEdgeNumbering()
****************************************************************************/
void AssignHexHedronEdgeNumbering()
/* Assign global edge numbering
(see "Mesh generation" chapter in the user's guide
to verify these formulas and try to work it out by hand) */
{
int EdgeNum,xx1,xx2,xx3,kk,dy1,dy2,divz,divz1,i,j,k;
{
for(k=0;k<=ZdiM-1;k++)
for(j=0;j<=YdiM-1;j++)
for(i=0;i<=XdiM-1;i++) {
EdgeNum = k*YdiM*XdiM + j*XdiM + i;
for(kk=0;kk<=17;kk++) {HexEdgeNum[EdgeNum][kk] = 0;}
divz=(6*XdiM*YdiM+3*XdiM+3*YdiM+1)*k;
divz1= 6*XdiM*YdiM + 3*XdiM + 3*YdiM + 1 ;
dy1 = (3*XdiM+1)*j; dy2 = (3*XdiM+2)*j;
xx1 = 1 + (YdiM+1)*XdiM + (2*XdiM+1)*YdiM;
xx2 = xx1 + (2*XdiM+1);
xx3 = xx2 + (XdiM + 1);
HexEdgeNum[EdgeNum][0] = 1 + i + dy1 + divz;
HexEdgeNum[EdgeNum][1] = 1 + XdiM + 2*i + dy1 + divz;
HexEdgeNum[EdgeNum][2] = HexEdgeNum[EdgeNum][1] + 1;
HexEdgeNum[EdgeNum][3] = HexEdgeNum[EdgeNum][2] + 1;
HexEdgeNum[EdgeNum][4] = 1 + 3*XdiM + 1 + i + dy1 + divz;
HexEdgeNum[EdgeNum][5] = xx1 + 2*i + dy2 + divz;
HexEdgeNum[EdgeNum][6] = HexEdgeNum[EdgeNum][5] + 1;
HexEdgeNum[EdgeNum][7] = HexEdgeNum[EdgeNum][6] + 1;
HexEdgeNum[EdgeNum][8] = xx2 + i + dy2 + divz;
HexEdgeNum[EdgeNum][9] = HexEdgeNum[EdgeNum][8] + 1;
HexEdgeNum[EdgeNum][10]= xx3 + 2*i + dy2 + divz;
HexEdgeNum[EdgeNum][11]= HexEdgeNum[EdgeNum][10] + 1;
HexEdgeNum[EdgeNum][12]= HexEdgeNum[EdgeNum][11] + 1;
HexEdgeNum[EdgeNum][13]= HexEdgeNum[EdgeNum][0] + divz1;
HexEdgeNum[EdgeNum][14]= HexEdgeNum[EdgeNum][1] + divz1;
HexEdgeNum[EdgeNum][15]= HexEdgeNum[EdgeNum][2] + divz1;
HexEdgeNum[EdgeNum][16]= HexEdgeNum[EdgeNum][3] + divz1;
HexEdgeNum[EdgeNum][17]= HexEdgeNum[EdgeNum][4] + divz1;
}
}
}
/***************************************************************************
* FUNCTION : HexHedraSubDiv1()
****************************************************************************/
void HexHedraSubDiv1()
/* Divides the 8 nodes in the hexahedron into five groups of 4 each */
{
TetGlobalNodeNum[0][0]=HexNode[HexNum][0];
TetGlobalNodeNum[1][0]=HexNode[HexNum][3];
TetGlobalNodeNum[0][1]=HexNode[HexNum][3];
TetGlobalNodeNum[1][1]=HexNode[HexNum][5];
TetGlobalNodeNum[0][2]=HexNode[HexNum][2];
TetGlobalNodeNum[1][2]=HexNode[HexNum][7];
TetGlobalNodeNum[0][3]=HexNode[HexNum][6];
TetGlobalNodeNum[1][3]=HexNode[HexNum][6];
TetGlobalNodeNum[2][0]=HexNode[HexNum][0];
TetGlobalNodeNum[3][0]=HexNode[HexNum][0];
TetGlobalNodeNum[2][1]=HexNode[HexNum][1];
TetGlobalNodeNum[3][1]=HexNode[HexNum][6];
TetGlobalNodeNum[2][2]=HexNode[HexNum][3];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -