📄 radarfdtd.cpp
字号:
lowerY = hz[i][j-1][k];
else
lowerY = pmlHZXBottomXZ[i][pmlWidth-1][k] + pmlHZYBottomXZ[i][pmlWidth-1][k];
if (k > 0)
lowerZ = hy[i][j][k-1];
else
lowerZ = pmlHYXBottomXY[i][j+pmlWidth][pmlWidth-1] + pmlHYZBottomXY[i][j+pmlWidth][pmlWidth-1];
switch(material[i][j][k])
{
case 0:
ex[i][j][k] = ex[i][j][k] + materialConstants[0][DTEPSDY] * (hz[i][j][k] - lowerY) - materialConstants[0][DTEPSDZ] * (hy[i][j][k] - lowerZ);
break;
case 1:
ex[i][j][k] = 0;
break;
default:
ex[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ex[i][j][k] + materialConstants[material[i][j][k]][DTEPSDY]*(hz[i][j][k] - lowerY) - materialConstants[material[i][j][k]][DTEPSDZ]*(hy[i][j][k] - lowerZ);
}
}
}
/* Update the ey values: */
void CradarFDTD::UpdateInnerEy(void)
{
int i,j,k;
PRECISION lowerX, lowerZ;
for(i=0; i<(nx); i++)
for(j=0; j<(ny); j++)
for(k=0; k<(nz); k++)
{
if (i > 0)
lowerX = hz[i-1][j][k];
else
lowerX = pmlHZXBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth] + pmlHZYBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth];
if (k > 0)
lowerZ = hx[i][j][k-1];
else
lowerZ = pmlHXYBottomXY[i][j+pmlWidth][pmlWidth-1] + pmlHXZBottomXY[i][j+pmlWidth][pmlWidth-1];
switch(material[i][j][k])
{
case 0:
ey[i][j][k] = ey[i][j][k] + materialConstants[0][DTEPSDZ] * (hx[i][j][k] - lowerZ) - materialConstants[0][DTEPSDX] * (hz[i][j][k] - lowerX);
break;
case 1:
ey[i][j][k] = 0;
break;
default:
ey[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ey[i][j][k] + materialConstants[material[i][j][k]][DTEPSDZ]*(hx[i][j][k] - lowerZ) - materialConstants[material[i][j][k]][DTEPSDX]*(hz[i][j][k] - lowerX);
}
}
}
/* Update the ez values: */
void CradarFDTD::UpdateInnerEz(void)
{
int i,j,k;
PRECISION lowerX, lowerY;
for(i=0; i<(nx); i++)
for(j=0; j<(ny); j++)
for(k=0; k<(nz); k++)
{
if (i > 0)
lowerX = hy[i-1][j][k];
else
lowerX = pmlHYXBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth] + pmlHYZBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth];
if (j > 0)
lowerY = hx[i][j-1][k];
else
lowerY = pmlHXYBottomXZ[i][pmlWidth-1][k] + pmlHXZBottomXZ[i][pmlWidth-1][k];
switch(material[i][j][k])
{
case 0:
ez[i][j][k] = ez[i][j][k] + materialConstants[0][DTEPSDX] * (hy[i][j][k] - lowerX) - materialConstants[0][DTEPSDY] * (hx[i][j][k] - lowerY);
break;
case 1:
ez[i][j][k] = 0;
break;
default:
ez[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ez[i][j][k] + materialConstants[material[i][j][k]][DTEPSDX]*(hy[i][j][k] - lowerX) - materialConstants[material[i][j][k]][DTEPSDY]*(hx[i][j][k] - lowerY);
}
}
}
void CradarFDTD::UpdateInnerHx(void)
{
int i,j,k;
PRECISION nextY, nextZ; /* the upper neighbour-values in y and z direction */
for(i=0; i<(nx); i++)
for(j=0; j<(ny); j++)
for(k=0; k<(nz); k++)
{
if (j < ny-1)
nextY = ez[i][j+1][k];
else
nextY = pmlEZXTopXZ[i][0][k] + pmlEZYTopXZ[i][0][k];
if (k < nz-1)
nextZ = ey[i][j][k+1];
else
nextZ = pmlEYXTopXY[i][j+pmlWidth][0] + pmlEYZTopXY[i][j+pmlWidth][0];
hx[i][j][k] = hx[i][j][k] + (materialConstants[material[i][j][k]][DTMUDZ]*(nextZ - ey[i][j][k]) - materialConstants[material[i][j][k]][DTMUDY]*(nextY - ez[i][j][k]));
}
}
/* Update the hy values: */
void CradarFDTD::UpdateInnerHy(void)
{
int i,j,k;
PRECISION nextX, nextZ; /* the upper neighbour-values in y and z direction */
for(i=0; i<(nx); i++)
for(j=0; j<(ny); j++)
for(k=0; k<(nz); k++)
{
if (i < nx-1)
nextX = ez[i+1][j][k];
else
nextX = pmlEZXTopYZ[0][j+pmlWidth][k+pmlWidth] + pmlEZYTopYZ[0][j+pmlWidth][k+pmlWidth];
if (k < nz-1)
nextZ = ex[i][j][k+1];
else
nextZ = pmlEXYTopXY[i][j+pmlWidth][0] + pmlEXZTopXY[i][j+pmlWidth][0];
hy[i][j][k] = hy[i][j][k] + (materialConstants[material[i][j][k]][DTMUDX]*(nextX - ez[i][j][k]) - materialConstants[material[i][j][k]][DTMUDZ]*(nextZ - ex[i][j][k]));
}
}
/* Update the hz values: */
void CradarFDTD::UpdateInnerHz(void)
{
int i,j,k;
PRECISION nextX, nextY;
for(i=0; i<(nx); i++)
for(j=0; j<(ny); j++)
for(k=0; k<(nz); k++)
{
if (i < nx-1)
nextX = ey[i+1][j][k];
else
nextX = pmlEYXTopYZ[0][j+pmlWidth][k+pmlWidth] + pmlEYZTopYZ[0][j+pmlWidth][k+pmlWidth];
if (j < ny-1)
nextY = ex[i][j+1][k];
else
nextY = pmlEXYTopXZ[i][0][k] + pmlEXZTopXZ[i][0][k];
hz[i][j][k] = hz[i][j][k] + (materialConstants[material[i][j][k]][DTMUDY]*(nextY - ex[i][j][k]) - materialConstants[material[i][j][k]][DTMUDX]*(nextX - ey[i][j][k]));
}
}
void CradarFDTD::UpdatePMLComplexE(void)
{
/* speed up variables */
int ny_plus_2pmlWidth = ny + 2*pmlWidth;
int nz_plus_2pmlWidth = nz + 2*pmlWidth;
/* Bottom YZ */
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny_plus_2pmlWidth; j++)
for(k=0; k<nz_plus_2pmlWidth; k++)
{
pmlComplexEXYBottomYZ[i][j][k] += dt * pmlEXYBottomYZ[i][j][k];
pmlComplexEXZBottomYZ[i][j][k] += dt * pmlEXZBottomYZ[i][j][k];
pmlComplexEYXBottomYZ[i][j][k] += dt * pmlEYXBottomYZ[i][j][k];
pmlComplexEYZBottomYZ[i][j][k] += dt * pmlEYZBottomYZ[i][j][k];
pmlComplexEZXBottomYZ[i][j][k] += dt * pmlEZXBottomYZ[i][j][k];
pmlComplexEZYBottomYZ[i][j][k] += dt * pmlEZYBottomYZ[i][j][k];
}
/* Top YZ */
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny_plus_2pmlWidth; j++)
for(k=0; k<nz_plus_2pmlWidth; k++)
{
pmlComplexEXYTopYZ[i][j][k] += dt * pmlEXYTopYZ[i][j][k];
pmlComplexEXZTopYZ[i][j][k] += dt * pmlEXZTopYZ[i][j][k];
pmlComplexEYXTopYZ[i][j][k] += dt * pmlEYXTopYZ[i][j][k];
pmlComplexEYZTopYZ[i][j][k] += dt * pmlEYZTopYZ[i][j][k];
pmlComplexEZXTopYZ[i][j][k] += dt * pmlEZXTopYZ[i][j][k];
pmlComplexEZYTopYZ[i][j][k] += dt * pmlEZYTopYZ[i][j][k];
}
/* Bottom XY */
for(i=0; i<nx; i++)
for(j=0; j<ny_plus_2pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlComplexEXYBottomXY[i][j][k] += dt * pmlEXYBottomXY[i][j][k];
pmlComplexEXZBottomXY[i][j][k] += dt * pmlEXZBottomXY[i][j][k];
pmlComplexEYXBottomXY[i][j][k] += dt * pmlEYXBottomXY[i][j][k];
pmlComplexEYZBottomXY[i][j][k] += dt * pmlEYZBottomXY[i][j][k];
pmlComplexEZXBottomXY[i][j][k] += dt * pmlEZXBottomXY[i][j][k];
pmlComplexEZYBottomXY[i][j][k] += dt * pmlEZYBottomXY[i][j][k];
}
/* Top XY */
for(i=0; i<nx; i++)
for(j=0; j<ny_plus_2pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlComplexEXYTopXY[i][j][k] += dt * pmlEXYTopXY[i][j][k];
pmlComplexEXZTopXY[i][j][k] += dt * pmlEXZTopXY[i][j][k];
pmlComplexEYXTopXY[i][j][k] += dt * pmlEYXTopXY[i][j][k];
pmlComplexEYZTopXY[i][j][k] += dt * pmlEYZTopXY[i][j][k];
pmlComplexEZXTopXY[i][j][k] += dt * pmlEZXTopXY[i][j][k];
pmlComplexEZYTopXY[i][j][k] += dt * pmlEZYTopXY[i][j][k];
}
/* BottomXZ */
for(i=0; i<nx; i++)
for(j=0; j<pmlWidth; j++)
for(k=0; k<nz; k++)
{
pmlComplexEXYBottomXZ[i][j][k] += dt * pmlEXYBottomXZ[i][j][k];
pmlComplexEXZBottomXZ[i][j][k] += dt * pmlEXZBottomXZ[i][j][k];
pmlComplexEYXBottomXZ[i][j][k] += dt * pmlEYXBottomXZ[i][j][k];
pmlComplexEYZBottomXZ[i][j][k] += dt * pmlEYZBottomXZ[i][j][k];
pmlComplexEZXBottomXZ[i][j][k] += dt * pmlEZXBottomXZ[i][j][k];
pmlComplexEZYBottomXZ[i][j][k] += dt * pmlEZYBottomXZ[i][j][k];
}
/* TopXZ */
for(i=0; i<nx; i++)
for(j=0; j<pmlWidth; j++)
for(k=0; k<nz; k++)
{
pmlComplexEXYTopXZ[i][j][k] += dt * pmlEXYTopXZ[i][j][k];
pmlComplexEXZTopXZ[i][j][k] += dt * pmlEXZTopXZ[i][j][k];
pmlComplexEYXTopXZ[i][j][k] += dt * pmlEYXTopXZ[i][j][k];
pmlComplexEYZTopXZ[i][j][k] += dt * pmlEYZTopXZ[i][j][k];
pmlComplexEZXTopXZ[i][j][k] += dt * pmlEZXTopXZ[i][j][k];
pmlComplexEZYTopXZ[i][j][k] += dt * pmlEZYTopXZ[i][j][k];
}
}
void CradarFDTD::UpdatePMLHx(void)
{
int i,j,k;
PRECISION mudt = MU_0 / dt;
PRECISION topY1, topY2, topZ1, topZ2; /* components which may be in different memory blocks than the one worked on */
PRECISION commonFactor;
/* speed up variables */
int ny_plus_2pmlWidth = ny + 2*pmlWidth;
int nz_plus_2pmlWidth = nz + 2*pmlWidth;
int ny_plus_pmlWidth = ny + pmlWidth;
int pmlWidth_minus_1 = pmlWidth-1;
int nz_minus_1 = nz-1;
/* start with the YZ planes */
/* there are no problems at any boundary ! */
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny_plus_2pmlWidth; j++)
for(k=0; k<nz_plus_2pmlWidth; k++)
{ /* bottom */
commonFactor = mudt + 0.5 * pmlSigmaYStarBottomYZ[i][j][k];
pmlHXYBottomYZ[i][j][k] = pmlHXYBottomYZ[i][j][k] *
(mudt - 0.5 * pmlSigmaYStarBottomYZ[i][j][k]) / commonFactor -
(pmlEZXBottomYZ[i][j+1][k] - pmlEZXBottomYZ[i][j][k] +
pmlEZYBottomYZ[i][j+1][k] - pmlEZYBottomYZ[i][j][k]) /
(dy * pmlSYBottomYZ[i][j][k] * commonFactor);
commonFactor = mudt + 0.5 * pmlSigmaZStarBottomYZ[i][j][k];
pmlHXZBottomYZ[i][j][k] = pmlHXZBottomYZ[i][j][k] *
(mudt - 0.5 * pmlSigmaZStarBottomYZ[i][j][k]) / commonFactor +
(pmlEYZBottomYZ[i][j][k+1] - pmlEYZBottomYZ[i][j][k] +
pmlEYXBottomYZ[i][j][k+1] - pmlEYXBottomYZ[i][j][k]) /
(dz * pmlSZBottomYZ[i][j][k] * commonFactor);
}
/* top, there are no problems at any boundary */
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny_plus_2pmlWidth; j++)
for(k=0; k<nz_plus_2pmlWidth; k++)
{
commonFactor = mudt + 0.5 * pmlSigmaYStarBottomYZ[i][j][k];
pmlHXYTopYZ[i][j][k] = pmlHXYTopYZ[i][j][k] *
(mudt - 0.5 * pmlSigmaYStarTopYZ[i][j][k]) / commonFactor -
(pmlEZXTopYZ[i][j+1][k] - pmlEZXTopYZ[i][j][k] +
pmlEZYTopYZ[i][j+1][k] - pmlEZYTopYZ[i][j][k]) /
(dy * pmlSYTopYZ[i][j][k] * commonFactor);
commonFactor = mudt + 0.5 * pmlSigmaZStarTopYZ[i][j][k];
pmlHXZTopYZ[i][j][k] = pmlHXZTopYZ[i][j][k] *
(mudt - 0.5 * pmlSigmaZStarTopYZ[i][j][k]) / commonFactor +
(pmlEYZTopYZ[i][j][k+1] - pmlEYZTopYZ[i][j][k] +
pmlEYXTopYZ[i][j][k+1] - pmlEYXTopYZ[i][j][k]) /
(dz * pmlSZTopYZ[i][j][k] * commonFactor);
}
/* now the XY planes */
for(i=0; i<nx; i++)
for(j=0; j<ny_plus_2pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{ /* bottom */
/* there are three boundary-problems:
0 <= y < pmlWidth: XZ - plane (bottom)
pmlWid⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -