diff --git a/FDTD/extensions/engine_ext_pml_sf.cpp b/FDTD/extensions/engine_ext_pml_sf.cpp deleted file mode 100644 index 28defdf..0000000 --- a/FDTD/extensions/engine_ext_pml_sf.cpp +++ /dev/null @@ -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 . -*/ - -#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(m_Op_PML_SF->m_numLines); - volt[1] = Create_N_3DArray(m_Op_PML_SF->m_numLines); - curr[0] = Create_N_3DArray(m_Op_PML_SF->m_numLines); - curr[1] = Create_N_3DArray(m_Op_PML_SF->m_numLines); -} - -Engine_Ext_PML_SF::~Engine_Ext_PML_SF() -{ - Delete_N_3DArray(volt[0],m_Op_PML_SF->m_numLines); - volt[0]=NULL; - Delete_N_3DArray(volt[1],m_Op_PML_SF->m_numLines); - volt[1]=NULL; - Delete_N_3DArray(curr[0],m_Op_PML_SF->m_numLines); - curr[0]=NULL; - Delete_N_3DArray(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; posXm_numLines[1]; ++pos[1]) - { - shift[1]=pos[1]; - for (pos[2]=0; pos[2]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; posXm_numLines[1]-1; ++pos[1]) - { - for (pos[2]=0; pos[2]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_numLines[m_nyP]-1; ++pos[m_nyP]) - { - pml_pos[m_nyP] = pos[m_nyP]; - for (pos[m_nyPP]=0; pos[m_nyPP]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_numLines[m_nyP]; ++pos[m_nyP]) - { - pml_pos[m_nyP] = pos[m_nyP]; - for (pos[m_nyPP]=0; pos[m_nyPP]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_numLines[m_nyP]; ++pos[m_nyP]) - { - pml_pos[m_nyP] = pos[m_nyP]; - for (pos[m_nyPP]=0; pos[m_nyPP]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_numLines[m_nyP]-1; ++pos[m_nyP]) - { - pml_pos[m_nyP] = pos[m_nyP]; - for (pos[m_nyPP]=0; pos[m_nyPP]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_numLines[m_nyP]-1; ++pos[m_nyP]) - { - pml_pos[m_nyP] = pos[m_nyP]; - for (pos[m_nyPP]=0; pos[m_nyPP]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]] ); - } - } - } -} diff --git a/FDTD/extensions/engine_ext_pml_sf.h b/FDTD/extensions/engine_ext_pml_sf.h deleted file mode 100644 index 8e910d5..0000000 --- a/FDTD/extensions/engine_ext_pml_sf.h +++ /dev/null @@ -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 . -*/ - -#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 diff --git a/FDTD/extensions/operator_ext_pml_sf.cpp b/FDTD/extensions/operator_ext_pml_sf.cpp deleted file mode 100644 index d067cce..0000000 --- a/FDTD/extensions/operator_ext_pml_sf.cpp +++ /dev/null @@ -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 . -*/ - -#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(m_numLines); - vv[1] = Create_N_3DArray(m_numLines); - - vi[0] = Create_N_3DArray(m_numLines); - vi[1] = Create_N_3DArray(m_numLines); - - ii[0] = Create_N_3DArray(m_numLines); - ii[1] = Create_N_3DArray(m_numLines); - - iv[0] = Create_N_3DArray(m_numLines); - iv[1] = Create_N_3DArray(m_numLines); -} - - -void Operator_Ext_PML_SF::DeleteOP() -{ - if (!m_SetupDone) - return; - Delete_N_3DArray(vv[0],m_numLines); - vv[0] = NULL; - Delete_N_3DArray(vv[1],m_numLines); - vv[1] = NULL; - - Delete_N_3DArray(vi[0],m_numLines); - vi[0] = NULL; - Delete_N_3DArray(vi[1],m_numLines); - vi[1] = NULL; - - Delete_N_3DArray(ii[0],m_numLines); - ii[0] = NULL; - Delete_N_3DArray(ii[1],m_numLines); - ii[1] = NULL; - - Delete_N_3DArray(iv[0],m_numLines); - iv[0] = NULL; - Delete_N_3DArray(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]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_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; -} diff --git a/FDTD/extensions/operator_ext_pml_sf.h b/FDTD/extensions/operator_ext_pml_sf.h deleted file mode 100644 index 41a9a2a..0000000 --- a/FDTD/extensions/operator_ext_pml_sf.h +++ /dev/null @@ -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 . -*/ - -#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 diff --git a/openEMS.pro b/openEMS.pro index e62581f..a391be9 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -141,8 +141,6 @@ SOURCES += FDTD/extensions/engine_extension.cpp \ FDTD/extensions/operator_ext_conductingsheet.cpp \ FDTD/extensions/engine_ext_dispersive.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/operator_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/engine_ext_dispersive.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/operator_ext_upml.h \ FDTD/extensions/engine_ext_upml.h \ diff --git a/openems.cpp b/openems.cpp index 9f1539b..63f9bb1 100644 --- a/openems.cpp +++ b/openems.cpp @@ -28,7 +28,6 @@ #include "FDTD/extensions/operator_ext_excitation.h" #include "FDTD/extensions/operator_ext_tfsf.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_lorentzmaterial.h" #include "FDTD/extensions/operator_ext_conductingsheet.h"