remove old and unused split-field pml
Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>pull/1/head
Thorsten Liebig
parent
5f0261b3c5
commit
b8f1184071
|
|
@ -1,262 +0,0 @@
|
||||||
/*
|
|
||||||
* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
|
||||||
*
|
|
||||||
* This program is free software: you can redistribute it and/or modify
|
|
||||||
* it under the terms of the GNU General Public License as published by
|
|
||||||
* the Free Software Foundation, either version 3 of the License, or
|
|
||||||
* (at your option) any later version.
|
|
||||||
*
|
|
||||||
* This program is distributed in the hope that it will be useful,
|
|
||||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
* GNU General Public License for more details.
|
|
||||||
*
|
|
||||||
* You should have received a copy of the GNU General Public License
|
|
||||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include "engine_ext_pml_sf.h"
|
|
||||||
#include "operator_ext_pml_sf.h"
|
|
||||||
#include "FDTD/engine_sse.h"
|
|
||||||
#include "tools/array_ops.h"
|
|
||||||
|
|
||||||
/************************************************ Engine_Ext_PML_SF **************************************************************************/
|
|
||||||
Engine_Ext_PML_SF::Engine_Ext_PML_SF(Operator_Ext_PML_SF* op_ext) : Engine_Extension(op_ext)
|
|
||||||
{
|
|
||||||
m_Op_PML_SF = op_ext;
|
|
||||||
|
|
||||||
volt[0] = Create_N_3DArray<FDTD_FLOAT>(m_Op_PML_SF->m_numLines);
|
|
||||||
volt[1] = Create_N_3DArray<FDTD_FLOAT>(m_Op_PML_SF->m_numLines);
|
|
||||||
curr[0] = Create_N_3DArray<FDTD_FLOAT>(m_Op_PML_SF->m_numLines);
|
|
||||||
curr[1] = Create_N_3DArray<FDTD_FLOAT>(m_Op_PML_SF->m_numLines);
|
|
||||||
}
|
|
||||||
|
|
||||||
Engine_Ext_PML_SF::~Engine_Ext_PML_SF()
|
|
||||||
{
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(volt[0],m_Op_PML_SF->m_numLines);
|
|
||||||
volt[0]=NULL;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(volt[1],m_Op_PML_SF->m_numLines);
|
|
||||||
volt[1]=NULL;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(curr[0],m_Op_PML_SF->m_numLines);
|
|
||||||
curr[0]=NULL;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(curr[1],m_Op_PML_SF->m_numLines);
|
|
||||||
curr[1]=NULL;
|
|
||||||
}
|
|
||||||
|
|
||||||
void Engine_Ext_PML_SF::UpdateVoltages(unsigned int startX, unsigned int numX)
|
|
||||||
{
|
|
||||||
unsigned int pos[3];
|
|
||||||
bool shift[3];
|
|
||||||
|
|
||||||
pos[0] = startX;
|
|
||||||
//voltage updates
|
|
||||||
for (unsigned int posX=0; posX<numX; ++posX)
|
|
||||||
{
|
|
||||||
shift[0]=pos[0];
|
|
||||||
for (pos[1]=0; pos[1]<m_Op_PML_SF->m_numLines[1]; ++pos[1])
|
|
||||||
{
|
|
||||||
shift[1]=pos[1];
|
|
||||||
for (pos[2]=0; pos[2]<m_Op_PML_SF->m_numLines[2]; ++pos[2])
|
|
||||||
{
|
|
||||||
shift[2]=pos[2];
|
|
||||||
//do the updates here
|
|
||||||
//for x
|
|
||||||
volt[0][0][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->vv[0][0][pos[0]][pos[1]][pos[2]];
|
|
||||||
volt[0][0][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->vi[0][0][pos[0]][pos[1]][pos[2]] * ( curr[0][2][pos[0]][pos[1]][pos[2]] + curr[1][2][pos[0]][pos[1]][pos[2]] - curr[0][2][pos[0]][pos[1]-shift[1]][pos[2]] - curr[1][2][pos[0]][pos[1]-shift[1]][pos[2]]);
|
|
||||||
|
|
||||||
volt[1][0][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->vv[1][0][pos[0]][pos[1]][pos[2]];
|
|
||||||
volt[1][0][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->vi[1][0][pos[0]][pos[1]][pos[2]] * ( curr[0][1][pos[0]][pos[1]][pos[2]-shift[2]] + curr[1][1][pos[0]][pos[1]][pos[2]-shift[2]] - curr[0][1][pos[0]][pos[1]][pos[2]] - curr[1][1][pos[0]][pos[1]][pos[2]]);
|
|
||||||
|
|
||||||
//for y
|
|
||||||
volt[0][1][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->vv[0][1][pos[0]][pos[1]][pos[2]];
|
|
||||||
volt[0][1][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->vi[0][1][pos[0]][pos[1]][pos[2]] * ( curr[0][0][pos[0]][pos[1]][pos[2]] + curr[1][0][pos[0]][pos[1]][pos[2]] - curr[0][0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[1][0][pos[0]][pos[1]][pos[2]-shift[2]]);
|
|
||||||
|
|
||||||
volt[1][1][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->vv[1][1][pos[0]][pos[1]][pos[2]];
|
|
||||||
volt[1][1][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->vi[1][1][pos[0]][pos[1]][pos[2]] * ( curr[0][2][pos[0]-shift[0]][pos[1]][pos[2]] + curr[1][2][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][2][pos[0]][pos[1]][pos[2]] - curr[1][2][pos[0]][pos[1]][pos[2]]);
|
|
||||||
|
|
||||||
//for z
|
|
||||||
volt[0][2][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->vv[0][2][pos[0]][pos[1]][pos[2]];
|
|
||||||
volt[0][2][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->vi[0][2][pos[0]][pos[1]][pos[2]] * ( curr[0][1][pos[0]][pos[1]][pos[2]] + curr[1][1][pos[0]][pos[1]][pos[2]] - curr[0][1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[1][1][pos[0]-shift[0]][pos[1]][pos[2]]);
|
|
||||||
|
|
||||||
volt[1][2][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->vv[1][2][pos[0]][pos[1]][pos[2]];
|
|
||||||
volt[1][2][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->vi[1][2][pos[0]][pos[1]][pos[2]] * ( curr[0][0][pos[0]][pos[1]-shift[1]][pos[2]] + curr[1][0][pos[0]][pos[1]-shift[1]][pos[2]] - curr[0][0][pos[0]][pos[1]][pos[2]] - curr[1][0][pos[0]][pos[1]][pos[2]]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
++pos[0];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void Engine_Ext_PML_SF::UpdateCurrents(unsigned int startX, unsigned int numX)
|
|
||||||
{
|
|
||||||
unsigned int pos[3];
|
|
||||||
pos[0] = startX;
|
|
||||||
for (unsigned int posX=0; posX<numX; ++posX)
|
|
||||||
{
|
|
||||||
for (pos[1]=0; pos[1]<m_Op_PML_SF->m_numLines[1]-1; ++pos[1])
|
|
||||||
{
|
|
||||||
for (pos[2]=0; pos[2]<m_Op_PML_SF->m_numLines[2]-1; ++pos[2])
|
|
||||||
{
|
|
||||||
//do the updates here
|
|
||||||
//for x
|
|
||||||
curr[0][0][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->ii[0][0][pos[0]][pos[1]][pos[2]];
|
|
||||||
curr[0][0][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->iv[0][0][pos[0]][pos[1]][pos[2]] * ( volt[0][2][pos[0]][pos[1]][pos[2]] + volt[1][2][pos[0]][pos[1]][pos[2]] - volt[0][2][pos[0]][pos[1]+1][pos[2]] - volt[1][2][pos[0]][pos[1]+1][pos[2]]);
|
|
||||||
|
|
||||||
curr[1][0][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->ii[1][0][pos[0]][pos[1]][pos[2]];
|
|
||||||
curr[1][0][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->iv[1][0][pos[0]][pos[1]][pos[2]] * ( volt[0][1][pos[0]][pos[1]][pos[2]+1] + volt[1][1][pos[0]][pos[1]][pos[2]+1] - volt[0][1][pos[0]][pos[1]][pos[2]] - volt[1][1][pos[0]][pos[1]][pos[2]]);
|
|
||||||
|
|
||||||
// cerr << volt[0][0][pos[0]][pos[1]][pos[2]] << " ";
|
|
||||||
// cerr << volt[0][1][pos[0]][pos[1]][pos[2]] << " ";
|
|
||||||
// cerr << volt[0][2][pos[0]][pos[1]][pos[2]] << endl;
|
|
||||||
|
|
||||||
//for y
|
|
||||||
curr[0][1][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->ii[0][1][pos[0]][pos[1]][pos[2]];
|
|
||||||
curr[0][1][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->iv[0][1][pos[0]][pos[1]][pos[2]] * ( volt[0][0][pos[0]][pos[1]][pos[2]] + volt[1][0][pos[0]][pos[1]][pos[2]] - volt[0][0][pos[0]][pos[1]][pos[2]+1] - volt[1][0][pos[0]][pos[1]][pos[2]+1]);
|
|
||||||
|
|
||||||
curr[1][1][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->ii[1][1][pos[0]][pos[1]][pos[2]];
|
|
||||||
curr[1][1][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->iv[1][1][pos[0]][pos[1]][pos[2]] * ( volt[0][2][pos[0]+1][pos[1]][pos[2]] + volt[1][2][pos[0]+1][pos[1]][pos[2]] - volt[0][2][pos[0]][pos[1]][pos[2]] - volt[1][2][pos[0]][pos[1]][pos[2]]);
|
|
||||||
|
|
||||||
//for z
|
|
||||||
curr[0][2][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->ii[0][2][pos[0]][pos[1]][pos[2]];
|
|
||||||
curr[0][2][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->iv[0][2][pos[0]][pos[1]][pos[2]] * ( volt[0][1][pos[0]][pos[1]][pos[2]] + volt[1][1][pos[0]][pos[1]][pos[2]] - volt[0][1][pos[0]+1][pos[1]][pos[2]] - volt[1][1][pos[0]+1][pos[1]][pos[2]]);
|
|
||||||
|
|
||||||
curr[1][2][pos[0]][pos[1]][pos[2]] *= m_Op_PML_SF->ii[1][2][pos[0]][pos[1]][pos[2]];
|
|
||||||
curr[1][2][pos[0]][pos[1]][pos[2]] += m_Op_PML_SF->iv[1][2][pos[0]][pos[1]][pos[2]] * ( volt[0][0][pos[0]][pos[1]+1][pos[2]] + volt[1][0][pos[0]][pos[1]+1][pos[2]] - volt[0][0][pos[0]][pos[1]][pos[2]] - volt[1][0][pos[0]][pos[1]][pos[2]]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
++pos[0];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/************************************************ Engine_Ext_PML_SF_Plane **************************************************************************/
|
|
||||||
Engine_Ext_PML_SF_Plane::Engine_Ext_PML_SF_Plane(Operator_Ext_PML_SF_Plane* op_ext) : Engine_Ext_PML_SF(op_ext)
|
|
||||||
{
|
|
||||||
m_Op_PML_SF_PL = op_ext;
|
|
||||||
}
|
|
||||||
|
|
||||||
Engine_Ext_PML_SF_Plane::~Engine_Ext_PML_SF_Plane()
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
void Engine_Ext_PML_SF_Plane::Apply2Voltages()
|
|
||||||
{
|
|
||||||
unsigned int pos[] = {0,0,0};
|
|
||||||
unsigned int pml_pos[] = {0,0,0};
|
|
||||||
unsigned int m_ny = m_Op_PML_SF_PL->m_ny;
|
|
||||||
unsigned int m_nyP = m_Op_PML_SF_PL->m_nyP;
|
|
||||||
unsigned int m_nyPP = m_Op_PML_SF_PL->m_nyPP;
|
|
||||||
|
|
||||||
pos[m_ny] = 0;
|
|
||||||
pml_pos[m_ny] = m_Op_PML_SF_PL->m_numLines[m_ny]-1;
|
|
||||||
// copy currents data from main engine to pml engine (lowest main line to highest pml line)
|
|
||||||
if (m_Op_PML_SF_PL->m_top==false)
|
|
||||||
{
|
|
||||||
for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1; ++pos[m_nyP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyP] = pos[m_nyP];
|
|
||||||
for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1; ++pos[m_nyPP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyPP] = pos[m_nyPP];
|
|
||||||
curr[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetCurr(0,pos);
|
|
||||||
curr[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
curr[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetCurr(1,pos);
|
|
||||||
curr[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
curr[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetCurr(2,pos);
|
|
||||||
curr[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// do the voltage updates
|
|
||||||
UpdateVoltages(0,m_Op_PML_SF->m_numLines[0]);
|
|
||||||
|
|
||||||
if (m_Op_PML_SF_PL->m_top==false)
|
|
||||||
{
|
|
||||||
// copy voltage data from pml engine to main engine (highest pml line to lowest main line)
|
|
||||||
pos[m_ny] = 0;
|
|
||||||
pml_pos[m_ny] = m_Op_PML_SF_PL->m_numLines[m_ny]-1;
|
|
||||||
for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]; ++pos[m_nyP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyP] = pos[m_nyP];
|
|
||||||
for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]; ++pos[m_nyPP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyPP] = pos[m_nyPP];
|
|
||||||
m_Eng->SetVolt(0,pos, volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] );
|
|
||||||
m_Eng->SetVolt(1,pos, volt[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] );
|
|
||||||
m_Eng->SetVolt(2,pos, volt[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] );
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void Engine_Ext_PML_SF_Plane::Apply2Current()
|
|
||||||
{
|
|
||||||
unsigned int pos[] = {0,0,0};
|
|
||||||
unsigned int pml_pos[] = {0,0,0};
|
|
||||||
unsigned int m_ny = m_Op_PML_SF_PL->m_ny;
|
|
||||||
unsigned int m_nyP = m_Op_PML_SF_PL->m_nyP;
|
|
||||||
unsigned int m_nyPP = m_Op_PML_SF_PL->m_nyPP;
|
|
||||||
|
|
||||||
pos[m_ny] = m_Op_PML_SF_PL->m_LineNr;
|
|
||||||
pml_pos[m_ny] = 0;
|
|
||||||
// copy voltage data from main engine to pml engine (highest main line to lowest pml line)
|
|
||||||
if (m_Op_PML_SF_PL->m_top)
|
|
||||||
{
|
|
||||||
for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]; ++pos[m_nyP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyP] = pos[m_nyP];
|
|
||||||
for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]; ++pos[m_nyPP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyPP] = pos[m_nyPP];
|
|
||||||
volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(0,pos);
|
|
||||||
volt[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
volt[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(1,pos);
|
|
||||||
volt[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
volt[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(2,pos);
|
|
||||||
volt[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
pos[m_ny] = 0;
|
|
||||||
pml_pos[m_ny] = m_Op_PML_SF_PL->m_numLines[m_ny]-1;
|
|
||||||
// copy (back again) voltage data from main engine to pml engine (lowest main line to highest pml line)
|
|
||||||
// this is necessary to catch the voltage excitation on the lowest main voltage line...
|
|
||||||
if (m_Op_PML_SF_PL->m_top==false)
|
|
||||||
{
|
|
||||||
for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1; ++pos[m_nyP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyP] = pos[m_nyP];
|
|
||||||
for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1; ++pos[m_nyPP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyPP] = pos[m_nyPP];
|
|
||||||
volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(0,pos);
|
|
||||||
volt[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
volt[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(1,pos);
|
|
||||||
volt[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
volt[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(2,pos);
|
|
||||||
volt[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] = 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
UpdateCurrents(0,m_Op_PML_SF->m_numLines[0]-1);
|
|
||||||
|
|
||||||
pos[m_ny] = m_Op_PML_SF_PL->m_LineNr;
|
|
||||||
pml_pos[m_ny] = 0;
|
|
||||||
// copy currents data from pml engine to main engine (lowest pml line to highest main line)
|
|
||||||
if (m_Op_PML_SF_PL->m_top)
|
|
||||||
{
|
|
||||||
for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1; ++pos[m_nyP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyP] = pos[m_nyP];
|
|
||||||
for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1; ++pos[m_nyPP])
|
|
||||||
{
|
|
||||||
pml_pos[m_nyPP] = pos[m_nyPP];
|
|
||||||
|
|
||||||
m_Eng->SetCurr(0,pos, curr[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] );
|
|
||||||
m_Eng->SetCurr(1,pos, curr[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] );
|
|
||||||
m_Eng->SetCurr(2,pos, curr[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] );
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
@ -1,64 +0,0 @@
|
||||||
/*
|
|
||||||
* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
|
||||||
*
|
|
||||||
* This program is free software: you can redistribute it and/or modify
|
|
||||||
* it under the terms of the GNU General Public License as published by
|
|
||||||
* the Free Software Foundation, either version 3 of the License, or
|
|
||||||
* (at your option) any later version.
|
|
||||||
*
|
|
||||||
* This program is distributed in the hope that it will be useful,
|
|
||||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
* GNU General Public License for more details.
|
|
||||||
*
|
|
||||||
* You should have received a copy of the GNU General Public License
|
|
||||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef ENGINE_EXT_PML_SF_H
|
|
||||||
#define ENGINE_EXT_PML_SF_H
|
|
||||||
|
|
||||||
#include "engine_extension.h"
|
|
||||||
#include "FDTD/engine.h"
|
|
||||||
#include "FDTD/operator.h"
|
|
||||||
|
|
||||||
class Operator_Ext_PML_SF;
|
|
||||||
class Operator_Ext_PML_SF_Plane;
|
|
||||||
|
|
||||||
//! Split field pml engine base class
|
|
||||||
class Engine_Ext_PML_SF : public Engine_Extension
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
virtual ~Engine_Ext_PML_SF();
|
|
||||||
|
|
||||||
inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return volt[0][n][x][y][z]+volt[1][n][x][y][z]; }
|
|
||||||
inline virtual FDTD_FLOAT GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return curr[0][n][x][y][z]+curr[1][n][x][y][z]; }
|
|
||||||
|
|
||||||
protected:
|
|
||||||
Engine_Ext_PML_SF(Operator_Ext_PML_SF* op_ext);
|
|
||||||
|
|
||||||
void UpdateVoltages(unsigned int startX, unsigned int numX);
|
|
||||||
void UpdateCurrents(unsigned int startX, unsigned int numX);
|
|
||||||
|
|
||||||
Operator_Ext_PML_SF* m_Op_PML_SF;
|
|
||||||
|
|
||||||
//split field voltages and currents
|
|
||||||
FDTD_FLOAT**** volt[2];
|
|
||||||
FDTD_FLOAT**** curr[2];
|
|
||||||
};
|
|
||||||
|
|
||||||
//! Split field pml plane engine class
|
|
||||||
class Engine_Ext_PML_SF_Plane : public Engine_Ext_PML_SF
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
Engine_Ext_PML_SF_Plane(Operator_Ext_PML_SF_Plane* op_ext);
|
|
||||||
virtual ~Engine_Ext_PML_SF_Plane();
|
|
||||||
|
|
||||||
virtual void Apply2Voltages();
|
|
||||||
virtual void Apply2Current();
|
|
||||||
|
|
||||||
protected:
|
|
||||||
Operator_Ext_PML_SF_Plane* m_Op_PML_SF_PL;
|
|
||||||
};
|
|
||||||
|
|
||||||
#endif // ENGINE_EXT_PML_SF_H
|
|
||||||
|
|
@ -1,430 +0,0 @@
|
||||||
/*
|
|
||||||
* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
|
||||||
*
|
|
||||||
* This program is free software: you can redistribute it and/or modify
|
|
||||||
* it under the terms of the GNU General Public License as published by
|
|
||||||
* the Free Software Foundation, either version 3 of the License, or
|
|
||||||
* (at your option) any later version.
|
|
||||||
*
|
|
||||||
* This program is distributed in the hope that it will be useful,
|
|
||||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
* GNU General Public License for more details.
|
|
||||||
*
|
|
||||||
* You should have received a copy of the GNU General Public License
|
|
||||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include "operator_ext_pml_sf.h"
|
|
||||||
#include "engine_ext_pml_sf.h"
|
|
||||||
#include "FDTD/operator_cylinder.h"
|
|
||||||
#include "tools/array_ops.h"
|
|
||||||
#include "fparser.hh"
|
|
||||||
|
|
||||||
bool Build_Split_Field_PML(Operator* op, int BC[6], int size[6], string gradFunc)
|
|
||||||
{
|
|
||||||
for (int n=0; n<6; ++n)
|
|
||||||
{
|
|
||||||
if (BC[n]==3) //split field PML
|
|
||||||
{
|
|
||||||
Operator_Ext_PML_SF_Plane* op_pml_sf = new Operator_Ext_PML_SF_Plane(op);
|
|
||||||
op_pml_sf->SetDirection(n/2,n%2);
|
|
||||||
if ((size[n]<4) || (size[n]>50))
|
|
||||||
{
|
|
||||||
cerr << "Build_Split_Field_PML: Warning, pml size invalid, skipping pml..." << endl;
|
|
||||||
delete op_pml_sf;
|
|
||||||
continue;
|
|
||||||
}
|
|
||||||
op_pml_sf->SetPMLLength(size[n]);
|
|
||||||
op_pml_sf->SetBoundaryCondition(BC);
|
|
||||||
|
|
||||||
if (!op_pml_sf->SetGradingFunction(gradFunc))
|
|
||||||
{
|
|
||||||
cerr << "Build_Split_Field_PML: Warning, pml grading function invalid, skipping pml..." << endl;
|
|
||||||
delete op_pml_sf;
|
|
||||||
continue;
|
|
||||||
}
|
|
||||||
|
|
||||||
cerr << "Build_Split_Field_PML:: Warning, currently only pml planes are implemented... edges and corner coming soon..." << endl;
|
|
||||||
op->AddExtension(op_pml_sf);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
/************************************************ Operator_Ext_PML_SF **************************************************************************/
|
|
||||||
Operator_Ext_PML_SF::Operator_Ext_PML_SF(Operator* op) : Operator_Extension(op)
|
|
||||||
{
|
|
||||||
m_SetupDone = false;
|
|
||||||
|
|
||||||
m_numLines[0]=0;
|
|
||||||
m_numLines[1]=0;
|
|
||||||
m_numLines[2]=0;
|
|
||||||
|
|
||||||
vv[0] = NULL;
|
|
||||||
vv[1] = NULL;
|
|
||||||
vi[0] = NULL;
|
|
||||||
vi[1] = NULL;
|
|
||||||
ii[0] = NULL;
|
|
||||||
ii[1] = NULL;
|
|
||||||
iv[0] = NULL;
|
|
||||||
iv[1] = NULL;
|
|
||||||
|
|
||||||
for (int n=0; n<6; ++n)
|
|
||||||
m_BC[n]=0;
|
|
||||||
|
|
||||||
m_GradingFunction = new FunctionParser();
|
|
||||||
//default grading function
|
|
||||||
SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z ");
|
|
||||||
}
|
|
||||||
|
|
||||||
Operator_Ext_PML_SF::~Operator_Ext_PML_SF()
|
|
||||||
{
|
|
||||||
delete m_GradingFunction;
|
|
||||||
m_GradingFunction = NULL;
|
|
||||||
DeleteOP();
|
|
||||||
}
|
|
||||||
|
|
||||||
void Operator_Ext_PML_SF::InitOP()
|
|
||||||
{
|
|
||||||
if (!m_SetupDone)
|
|
||||||
return;
|
|
||||||
vv[0] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
vv[1] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
|
|
||||||
vi[0] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
vi[1] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
|
|
||||||
ii[0] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
ii[1] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
|
|
||||||
iv[0] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
iv[1] = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
void Operator_Ext_PML_SF::DeleteOP()
|
|
||||||
{
|
|
||||||
if (!m_SetupDone)
|
|
||||||
return;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(vv[0],m_numLines);
|
|
||||||
vv[0] = NULL;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(vv[1],m_numLines);
|
|
||||||
vv[1] = NULL;
|
|
||||||
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(vi[0],m_numLines);
|
|
||||||
vi[0] = NULL;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(vi[1],m_numLines);
|
|
||||||
vi[1] = NULL;
|
|
||||||
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(ii[0],m_numLines);
|
|
||||||
ii[0] = NULL;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(ii[1],m_numLines);
|
|
||||||
ii[1] = NULL;
|
|
||||||
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(iv[0],m_numLines);
|
|
||||||
iv[0] = NULL;
|
|
||||||
Delete_N_3DArray<FDTD_FLOAT>(iv[1],m_numLines);
|
|
||||||
iv[1] = NULL;
|
|
||||||
}
|
|
||||||
|
|
||||||
bool Operator_Ext_PML_SF::SetGradingFunction(string func)
|
|
||||||
{
|
|
||||||
if (func.empty())
|
|
||||||
return true;
|
|
||||||
|
|
||||||
m_GradFunc = func;
|
|
||||||
int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N");
|
|
||||||
if (res < 0) return true;
|
|
||||||
|
|
||||||
cerr << "Operator_Ext_PML_SF::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl;
|
|
||||||
cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n";
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
|
|
||||||
bool Operator_Ext_PML_SF::BuildExtension()
|
|
||||||
{
|
|
||||||
if (!m_SetupDone)
|
|
||||||
{
|
|
||||||
cerr << "Operator_Ext_PML_SF::BuildExtension: Warning, Extension not initialized! Abort build!!" << endl;
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
double dT = m_Op->GetTimestep();
|
|
||||||
unsigned int pos[] = {0,0,0};
|
|
||||||
DeleteOP();
|
|
||||||
InitOP();
|
|
||||||
|
|
||||||
double inEC[4];
|
|
||||||
|
|
||||||
for (int n=0; n<3; ++n)
|
|
||||||
{
|
|
||||||
for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
|
|
||||||
{
|
|
||||||
for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
|
|
||||||
{
|
|
||||||
for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
|
|
||||||
{
|
|
||||||
Calc_ECPos(0,n,pos,inEC);
|
|
||||||
if (inEC[0]>0)
|
|
||||||
GetVV(0,n,pos[0],pos[1],pos[2]) = (1-dT*inEC[1]/2/inEC[0])/(1+dT*inEC[1]/2/inEC[0]);
|
|
||||||
if (inEC[2]>0)
|
|
||||||
GetII(0,n,pos[0],pos[1],pos[2]) = (1-dT*inEC[3]/2/inEC[2])/(1+dT*inEC[3]/2/inEC[2]);
|
|
||||||
|
|
||||||
if (inEC[0]>0)
|
|
||||||
GetVI(0,n,pos[0],pos[1],pos[2]) = (dT/inEC[0])/(1+dT*inEC[1]/2/inEC[0]);
|
|
||||||
if (inEC[2]>0)
|
|
||||||
GetIV(0,n,pos[0],pos[1],pos[2]) = (dT/inEC[2])/(1+dT*inEC[3]/2/inEC[2]);
|
|
||||||
|
|
||||||
// if (n==0)
|
|
||||||
// cerr << pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[1] << endl;
|
|
||||||
|
|
||||||
Calc_ECPos(1,n,pos,inEC);
|
|
||||||
if (inEC[0]>0)
|
|
||||||
GetVV(1,n,pos[0],pos[1],pos[2]) = (1-dT*inEC[1]/2/inEC[0])/(1+dT*inEC[1]/2/inEC[0]);
|
|
||||||
if (inEC[2]>0)
|
|
||||||
GetII(1,n,pos[0],pos[1],pos[2]) = (1-dT*inEC[3]/2/inEC[2])/(1+dT*inEC[3]/2/inEC[2]);
|
|
||||||
|
|
||||||
if (inEC[0]>0)
|
|
||||||
GetVI(1,n,pos[0],pos[1],pos[2]) = (dT/inEC[0])/(1+dT*inEC[1]/2/inEC[0]);
|
|
||||||
if (inEC[2]>0)
|
|
||||||
GetIV(1,n,pos[0],pos[1],pos[2]) = (dT/inEC[2])/(1+dT*inEC[3]/2/inEC[2]);
|
|
||||||
|
|
||||||
// if (n==0)
|
|
||||||
// cerr << pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[1] << endl;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
ApplyBC();
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
/************************************************ Operator_Ext_PML_SF_Plane **************************************************************************/
|
|
||||||
Operator_Ext_PML_SF_Plane::Operator_Ext_PML_SF_Plane(Operator* op) : Operator_Ext_PML_SF(op)
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
Operator_Ext_PML_SF_Plane::~Operator_Ext_PML_SF_Plane()
|
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
void Operator_Ext_PML_SF_Plane::SetDirection(int ny, bool top_ny)
|
|
||||||
{
|
|
||||||
if ((ny<0) || (ny>2))
|
|
||||||
return;
|
|
||||||
|
|
||||||
m_ny = ny;
|
|
||||||
m_nyP = (ny+1)%3;
|
|
||||||
m_nyPP = (ny+2)%3;
|
|
||||||
|
|
||||||
m_top = top_ny;
|
|
||||||
|
|
||||||
m_numLines[m_ny] = 8; //default width of the pml plane
|
|
||||||
m_numLines[m_nyP] = m_Op->GetNumberOfLines(m_nyP);
|
|
||||||
m_numLines[m_nyPP] = m_Op->GetNumberOfLines(m_nyPP);
|
|
||||||
|
|
||||||
unsigned int pos[] = {0,0,0};
|
|
||||||
m_LineNr = (unsigned int)((int)m_top * (int)(m_Op->GetNumberOfLines(m_ny)-1));
|
|
||||||
pos[m_ny] = m_LineNr;
|
|
||||||
|
|
||||||
m_pml_delta = m_Op->GetEdgeLength(m_ny,pos);
|
|
||||||
}
|
|
||||||
|
|
||||||
void Operator_Ext_PML_SF_Plane::SetPMLLength(int width)
|
|
||||||
{
|
|
||||||
if (m_ny<0)
|
|
||||||
{
|
|
||||||
cerr << "Operator_Ext_PML_SF_Plane::SetPMLLength: Warning, Direction not set! Use SetDirection first!!" << endl;
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (width<4)
|
|
||||||
{
|
|
||||||
cerr << "Operator_Ext_PML_SF_Plane::SetPMLLength: Warning: A pml width smaller than 4 lines is not allowed, skipping..." << endl;
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
if (width>50)
|
|
||||||
{
|
|
||||||
cerr << "Operator_Ext_PML_SF_Plane::SetPMLLength: Warning: A pml width greater than 20 lines is not allowed, skipping..." << endl;
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
m_SetupDone = true;
|
|
||||||
m_numLines[m_ny] = width;
|
|
||||||
|
|
||||||
m_pml_width = (width - 1.5) * m_pml_delta;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
double Operator_Ext_PML_SF_Plane::GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh) const
|
|
||||||
{
|
|
||||||
unsigned int l_pos[] = {pos[0],pos[1],pos[2]};
|
|
||||||
l_pos[m_ny] = m_LineNr;
|
|
||||||
|
|
||||||
return m_Op->GetEdgeArea(ny,l_pos,dualMesh);
|
|
||||||
}
|
|
||||||
|
|
||||||
double Operator_Ext_PML_SF_Plane::GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh) const
|
|
||||||
{
|
|
||||||
if (ny==m_ny)
|
|
||||||
return m_pml_delta;
|
|
||||||
|
|
||||||
unsigned int l_pos[] = {pos[0],pos[1],pos[2]};
|
|
||||||
l_pos[m_ny] = m_LineNr;
|
|
||||||
return m_Op->GetEdgeLength(ny,l_pos,dualMesh);
|
|
||||||
}
|
|
||||||
|
|
||||||
double Operator_Ext_PML_SF_Plane::GetKappaGraded(double depth, double Zm) const
|
|
||||||
{
|
|
||||||
if (depth<0)
|
|
||||||
return 0.0;
|
|
||||||
|
|
||||||
double vars[5] = {depth, m_pml_delta, m_pml_width, Zm, (double)m_numLines[m_ny]};
|
|
||||||
return m_GradingFunction->Eval(vars);
|
|
||||||
}
|
|
||||||
|
|
||||||
bool Operator_Ext_PML_SF_Plane::Calc_ECPos(int nP, int n, unsigned int* pos, double* inEC) const
|
|
||||||
{
|
|
||||||
unsigned int l_pos[] = {pos[0],pos[1],pos[2]};
|
|
||||||
l_pos[m_ny] = m_LineNr;
|
|
||||||
|
|
||||||
double inMat[4];
|
|
||||||
m_Op->Calc_EffMatPos(n,l_pos,inMat);
|
|
||||||
|
|
||||||
double Zm2 = inMat[2] / inMat[0]; // Zm^2 = mue/eps
|
|
||||||
double Zm = sqrt(Zm2); // Zm = sqrt(Zm^2) = sqrt(mue/eps)
|
|
||||||
double kappa = 0;
|
|
||||||
double sigma = 0;
|
|
||||||
double depth = 0;
|
|
||||||
if ( (n + nP + 1)%3 == m_ny )
|
|
||||||
{
|
|
||||||
if (m_top)
|
|
||||||
{
|
|
||||||
depth = pos[m_ny]*m_pml_delta - 0.5*m_pml_delta;
|
|
||||||
kappa = GetKappaGraded(depth, Zm);
|
|
||||||
sigma = GetKappaGraded(depth + 0.5*m_pml_delta, Zm) * Zm2;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
depth = m_pml_width - (pos[m_ny])*m_pml_delta;
|
|
||||||
kappa = GetKappaGraded(depth, Zm) ;
|
|
||||||
sigma = GetKappaGraded(depth-0.5*m_pml_delta, Zm) * Zm2;
|
|
||||||
}
|
|
||||||
if ((inMat[0]<=0) || (inMat[2]<=0)) //check if material properties are valid (necessary for cylindrical coords)
|
|
||||||
{
|
|
||||||
kappa = sigma = 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
double geomFactor = GetEdgeArea(n,pos) / GetEdgeLength(n,pos);
|
|
||||||
if (geomFactor<=0 || isnan(geomFactor) || isinf(geomFactor)) //check if geomFactor is positive, not zero and a valid number (necessary for cylindrical coords)
|
|
||||||
geomFactor = 0;
|
|
||||||
inEC[0] = inMat[0] * geomFactor;
|
|
||||||
inEC[1] = (inMat[1]+kappa) * geomFactor;
|
|
||||||
|
|
||||||
geomFactor = GetEdgeArea(n,pos,true) / GetEdgeLength(n,pos,true);
|
|
||||||
if (geomFactor<=0 || isnan(geomFactor) || isinf(geomFactor)) //check if geomFactor is positive, not zero and a valid number (necessary for cylindrical coords)
|
|
||||||
geomFactor = 0;
|
|
||||||
inEC[2] = inMat[2] * geomFactor;
|
|
||||||
inEC[3] = (inMat[3]+sigma) * geomFactor;
|
|
||||||
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
void Operator_Ext_PML_SF_Plane::ApplyBC()
|
|
||||||
{
|
|
||||||
bool PEC[6] = {1,1,1,1,1,1};
|
|
||||||
bool PMC[6] = {0,0,0,0,0,0};
|
|
||||||
|
|
||||||
if (m_top==false)
|
|
||||||
PEC[2*m_ny+1] = 0;
|
|
||||||
|
|
||||||
for (int n=0; n<6; ++n)
|
|
||||||
{
|
|
||||||
PMC[n] = (m_BC[n] == 1);
|
|
||||||
if (n/2 == m_ny)
|
|
||||||
PMC[n] = false;
|
|
||||||
}
|
|
||||||
|
|
||||||
//apply BC
|
|
||||||
unsigned int pos[3] = {0,0,0};
|
|
||||||
for (int n=0; n<3; ++n)
|
|
||||||
{
|
|
||||||
int nP = (n+1)%3;
|
|
||||||
int nPP = (n+2)%3;
|
|
||||||
for (pos[nP]=0; pos[nP]<m_numLines[nP]; ++pos[nP])
|
|
||||||
{
|
|
||||||
for (pos[nPP]=0; pos[nPP]<m_numLines[nPP]; ++pos[nPP])
|
|
||||||
{
|
|
||||||
for (int m=0; m<2; ++m)
|
|
||||||
{
|
|
||||||
pos[n]=0;
|
|
||||||
GetVV(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n];
|
|
||||||
GetVI(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n];
|
|
||||||
GetVV(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n];
|
|
||||||
GetVI(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n];
|
|
||||||
|
|
||||||
GetII(m,n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n];
|
|
||||||
GetIV(m,n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n];
|
|
||||||
GetII(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n];
|
|
||||||
GetIV(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n];
|
|
||||||
GetII(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n];
|
|
||||||
GetIV(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n];
|
|
||||||
|
|
||||||
pos[n]=m_numLines[n]-1;
|
|
||||||
GetVV(m,n,pos[0],pos[1],pos[2]) = 0; // these are outside the FDTD-domain as defined by the main disc
|
|
||||||
GetVI(m,n,pos[0],pos[1],pos[2]) = 0; // these are outside the FDTD-domain as defined by the main disc
|
|
||||||
|
|
||||||
GetVV(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n+1];
|
|
||||||
GetVI(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n+1];
|
|
||||||
GetVV(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n+1];
|
|
||||||
GetVI(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n+1];
|
|
||||||
|
|
||||||
pos[n]=m_numLines[n]-2;
|
|
||||||
GetII(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n+1];
|
|
||||||
GetIV(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n+1];
|
|
||||||
GetII(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n+1];
|
|
||||||
GetIV(m,nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PMC[2*n+1];
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
Engine_Extension* Operator_Ext_PML_SF_Plane::CreateEngineExtention()
|
|
||||||
{
|
|
||||||
Engine_Ext_PML_SF_Plane* eng_ext = new Engine_Ext_PML_SF_Plane(this);
|
|
||||||
return eng_ext;
|
|
||||||
}
|
|
||||||
|
|
||||||
bool Operator_Ext_PML_SF_Plane::IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const
|
|
||||||
{
|
|
||||||
UNUSED(closedAlpha);
|
|
||||||
UNUSED(R0_included);
|
|
||||||
|
|
||||||
if (m_ny==2)
|
|
||||||
{
|
|
||||||
Operator_Cylinder* op_cyl = dynamic_cast<Operator_Cylinder*>(m_Op);
|
|
||||||
if (op_cyl==NULL)
|
|
||||||
{
|
|
||||||
cerr << "Operator_Ext_PML_SF_Plane::IsCylinderCoordsSave(): Error!!! Sanity check failed!!! ==> Developer is not sane.... this should never have happend.. exit..." << endl;
|
|
||||||
exit(0);
|
|
||||||
}
|
|
||||||
if (op_cyl->GetClosedAlpha())
|
|
||||||
{
|
|
||||||
cerr << "Operator_Ext_PML_SF_Plane::IsCylinderCoordsSave(): Warning... this extension can not handle a closed alpha cylinder operator... " << endl;
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
|
|
||||||
void Operator_Ext_PML_SF_Plane::ShowStat(ostream &ostr) const
|
|
||||||
{
|
|
||||||
Operator_Extension::ShowStat(ostr);
|
|
||||||
string XYZ[3] = {"x","y","z"};
|
|
||||||
string top_bot[2] = {"bottom", "top"};
|
|
||||||
ostr << " Active direction\t: " << XYZ[m_ny] << " (" << top_bot[m_top] << ")" << endl;
|
|
||||||
ostr << " PML width (cells)\t: " << m_numLines[m_ny] << endl;
|
|
||||||
ostr << " Grading function\t: \"" << m_GradFunc << "\"" << endl;
|
|
||||||
}
|
|
||||||
|
|
@ -1,134 +0,0 @@
|
||||||
/*
|
|
||||||
* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
|
||||||
*
|
|
||||||
* This program is free software: you can redistribute it and/or modify
|
|
||||||
* it under the terms of the GNU General Public License as published by
|
|
||||||
* the Free Software Foundation, either version 3 of the License, or
|
|
||||||
* (at your option) any later version.
|
|
||||||
*
|
|
||||||
* This program is distributed in the hope that it will be useful,
|
|
||||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
* GNU General Public License for more details.
|
|
||||||
*
|
|
||||||
* You should have received a copy of the GNU General Public License
|
|
||||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef OPERATOR_EXT_PML_SF_H
|
|
||||||
#define OPERATOR_EXT_PML_SF_H
|
|
||||||
|
|
||||||
#include "operator_extension.h"
|
|
||||||
#include "FDTD/operator.h"
|
|
||||||
|
|
||||||
class FunctionParser;
|
|
||||||
|
|
||||||
//! Insert split field pml planes, edges and corner as necessary by the given boundary conditions
|
|
||||||
bool Build_Split_Field_PML(Operator* op, int BC[6], int size[6], string gradFunc);
|
|
||||||
|
|
||||||
//! This is the abstract operator extension for truncating the FDTD domain with a split field pml
|
|
||||||
class Operator_Ext_PML_SF : public Operator_Extension
|
|
||||||
{
|
|
||||||
friend class Engine_Ext_PML_SF;
|
|
||||||
friend class Engine_Ext_PML_SF_Plane;
|
|
||||||
public:
|
|
||||||
~Operator_Ext_PML_SF();
|
|
||||||
|
|
||||||
virtual void SetBoundaryCondition(int* BCs) {for (int n=0; n<6; ++n) m_BC[n]=BCs[n];}
|
|
||||||
|
|
||||||
inline virtual FDTD_FLOAT& GetVV(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[nP][n][x][y][z]; }
|
|
||||||
inline virtual FDTD_FLOAT& GetVI(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[nP][n][x][y][z]; }
|
|
||||||
inline virtual FDTD_FLOAT& GetII(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[nP][n][x][y][z]; }
|
|
||||||
inline virtual FDTD_FLOAT& GetIV(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[nP][n][x][y][z]; }
|
|
||||||
|
|
||||||
virtual double GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
|
|
||||||
virtual double GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
|
|
||||||
|
|
||||||
//! This will resturn the pml parameter grading
|
|
||||||
virtual double GetKappaGraded(double depth, double Zm) const {UNUSED(depth); UNUSED(Zm); return 0.0;}
|
|
||||||
|
|
||||||
//! Set the grading function for the pml
|
|
||||||
/*!
|
|
||||||
Define the pml grading grading function.
|
|
||||||
Predefined variables in this grading function are:
|
|
||||||
D = depth in the pml in meter
|
|
||||||
dl = mesh delta inside the pml in meter
|
|
||||||
W = width (length) of the pml in meter
|
|
||||||
N = number of cells for the pml
|
|
||||||
Z = wave impedance at the current depth and position
|
|
||||||
example: SetGradingFunction("-log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z");
|
|
||||||
|
|
||||||
An empty function string will be ignored.
|
|
||||||
*/
|
|
||||||
virtual bool SetGradingFunction(string func);
|
|
||||||
|
|
||||||
virtual bool BuildExtension();
|
|
||||||
|
|
||||||
virtual string GetExtensionName() const {return string("Split Field PML Extension");}
|
|
||||||
|
|
||||||
// virtual void ShowStat(ostream &ostr) const;
|
|
||||||
|
|
||||||
protected:
|
|
||||||
Operator_Ext_PML_SF(Operator* op);
|
|
||||||
|
|
||||||
virtual void ApplyBC() {};
|
|
||||||
|
|
||||||
virtual bool Calc_ECPos(int nP, int n, unsigned int* pos, double* inEC) const {UNUSED(n); UNUSED(nP); UNUSED(pos); UNUSED(inEC); return true;};
|
|
||||||
|
|
||||||
unsigned int m_numLines[3];
|
|
||||||
bool m_SetupDone;
|
|
||||||
|
|
||||||
int m_BC[6];
|
|
||||||
|
|
||||||
string m_GradFunc;
|
|
||||||
FunctionParser* m_GradingFunction;
|
|
||||||
|
|
||||||
void InitOP();
|
|
||||||
void DeleteOP();
|
|
||||||
//split field EC operator
|
|
||||||
//the first array-index is the splitted field part
|
|
||||||
FDTD_FLOAT**** vv[2]; //calc new voltage from old voltage
|
|
||||||
FDTD_FLOAT**** vi[2]; //calc new voltage from old current
|
|
||||||
FDTD_FLOAT**** ii[2]; //calc new current from old current
|
|
||||||
FDTD_FLOAT**** iv[2]; //calc new current from old voltage
|
|
||||||
};
|
|
||||||
|
|
||||||
//! This is an operator extension for truncating the FDTD domain with a split field pml plane
|
|
||||||
class Operator_Ext_PML_SF_Plane : public Operator_Ext_PML_SF
|
|
||||||
{
|
|
||||||
friend class Engine_Ext_PML_SF_Plane;
|
|
||||||
public:
|
|
||||||
Operator_Ext_PML_SF_Plane(Operator* op);
|
|
||||||
~Operator_Ext_PML_SF_Plane();
|
|
||||||
|
|
||||||
//! Define the direction of this PML plane: ny=0,1,2 -> x,y,z and if at bottom_ny -> e.g. x=0 or x=end
|
|
||||||
void SetDirection(int ny, bool top_ny);
|
|
||||||
void SetPMLLength(int width);
|
|
||||||
|
|
||||||
virtual double GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh = false) const;
|
|
||||||
virtual double GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh = false) const;
|
|
||||||
|
|
||||||
virtual double GetKappaGraded(double depth, double Zm) const;
|
|
||||||
|
|
||||||
virtual bool Calc_ECPos(int nP, int n, unsigned int* pos, double* inEC) const;
|
|
||||||
|
|
||||||
virtual Engine_Extension* CreateEngineExtention();
|
|
||||||
|
|
||||||
virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const;
|
|
||||||
|
|
||||||
virtual string GetExtensionName() const {return string("Split Field PML Plane Extension");}
|
|
||||||
|
|
||||||
virtual void ShowStat(ostream &ostr) const;
|
|
||||||
|
|
||||||
protected:
|
|
||||||
virtual void ApplyBC();
|
|
||||||
|
|
||||||
int m_ny,m_nyP,m_nyPP;
|
|
||||||
bool m_top;
|
|
||||||
unsigned int m_LineNr;
|
|
||||||
|
|
||||||
double m_pml_delta;
|
|
||||||
double m_pml_width;
|
|
||||||
};
|
|
||||||
|
|
||||||
#endif // OPERATOR_EXT_PML_SF_H
|
|
||||||
|
|
@ -141,8 +141,6 @@ SOURCES += FDTD/extensions/engine_extension.cpp \
|
||||||
FDTD/extensions/operator_ext_conductingsheet.cpp \
|
FDTD/extensions/operator_ext_conductingsheet.cpp \
|
||||||
FDTD/extensions/engine_ext_dispersive.cpp \
|
FDTD/extensions/engine_ext_dispersive.cpp \
|
||||||
FDTD/extensions/engine_ext_lorentzmaterial.cpp \
|
FDTD/extensions/engine_ext_lorentzmaterial.cpp \
|
||||||
FDTD/extensions/operator_ext_pml_sf.cpp \
|
|
||||||
FDTD/extensions/engine_ext_pml_sf.cpp \
|
|
||||||
FDTD/extensions/engine_ext_cylindermultigrid.cpp \
|
FDTD/extensions/engine_ext_cylindermultigrid.cpp \
|
||||||
FDTD/extensions/operator_ext_upml.cpp \
|
FDTD/extensions/operator_ext_upml.cpp \
|
||||||
FDTD/extensions/engine_ext_upml.cpp \
|
FDTD/extensions/engine_ext_upml.cpp \
|
||||||
|
|
@ -214,8 +212,6 @@ HEADERS += FDTD/extensions/operator_extension.h \
|
||||||
FDTD/extensions/cond_sheet_parameter.h \
|
FDTD/extensions/cond_sheet_parameter.h \
|
||||||
FDTD/extensions/engine_ext_dispersive.h \
|
FDTD/extensions/engine_ext_dispersive.h \
|
||||||
FDTD/extensions/engine_ext_lorentzmaterial.h \
|
FDTD/extensions/engine_ext_lorentzmaterial.h \
|
||||||
FDTD/extensions/operator_ext_pml_sf.h \
|
|
||||||
FDTD/extensions/engine_ext_pml_sf.h \
|
|
||||||
FDTD/extensions/engine_ext_cylindermultigrid.h \
|
FDTD/extensions/engine_ext_cylindermultigrid.h \
|
||||||
FDTD/extensions/operator_ext_upml.h \
|
FDTD/extensions/operator_ext_upml.h \
|
||||||
FDTD/extensions/engine_ext_upml.h \
|
FDTD/extensions/engine_ext_upml.h \
|
||||||
|
|
|
||||||
|
|
@ -28,7 +28,6 @@
|
||||||
#include "FDTD/extensions/operator_ext_excitation.h"
|
#include "FDTD/extensions/operator_ext_excitation.h"
|
||||||
#include "FDTD/extensions/operator_ext_tfsf.h"
|
#include "FDTD/extensions/operator_ext_tfsf.h"
|
||||||
#include "FDTD/extensions/operator_ext_mur_abc.h"
|
#include "FDTD/extensions/operator_ext_mur_abc.h"
|
||||||
#include "FDTD/extensions/operator_ext_pml_sf.h"
|
|
||||||
#include "FDTD/extensions/operator_ext_upml.h"
|
#include "FDTD/extensions/operator_ext_upml.h"
|
||||||
#include "FDTD/extensions/operator_ext_lorentzmaterial.h"
|
#include "FDTD/extensions/operator_ext_lorentzmaterial.h"
|
||||||
#include "FDTD/extensions/operator_ext_conductingsheet.h"
|
#include "FDTD/extensions/operator_ext_conductingsheet.h"
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue