diff --git a/FDTD/engine_ext_pml_sf.cpp b/FDTD/engine_ext_pml_sf.cpp
new file mode 100644
index 0000000..0a68de6
--- /dev/null
+++ b/FDTD/engine_ext_pml_sf.cpp
@@ -0,0 +1,241 @@
+/*
+* 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 "engine.h"
+#include "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->GetVolt(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->GetVolt(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->GetVolt(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;
+ }
+ }
+ }
+
+ 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->GetCurr(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->GetCurr(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->GetCurr(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/engine_ext_pml_sf.h b/FDTD/engine_ext_pml_sf.h
new file mode 100644
index 0000000..615beab
--- /dev/null
+++ b/FDTD/engine_ext_pml_sf.h
@@ -0,0 +1,64 @@
+/*
+* 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 "engine.h"
+#include "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/operator.cpp b/FDTD/operator.cpp
index 1c219d9..8e2d2a1 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -562,8 +562,11 @@ int Operator::CalcECOperator()
delete[] EC_R[n];EC_R[n]=NULL;
}
- //Always apply PEC to all boundary's
+ //Apply PEC to all boundary's
bool PEC[6]={1,1,1,1,1,1};
+ //exception for pml boundaries
+ for (int n=0;n<6;++n)
+ PEC[n] = m_BC[n]!=3;
ApplyElectricBC(PEC);
InitExcitation();
@@ -597,9 +600,11 @@ void Operator::ApplyElectricBC(bool* dirs)
GetVI(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
GetVV(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
GetVI(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
+
pos[n]=numLines[n]-1;
- GetVV(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
- GetVI(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
+ GetVV(n,pos[0],pos[1],pos[2]) = 0; // these are outside the FDTD-domain as defined by the main disc
+ GetVI(n,pos[0],pos[1],pos[2]) = 0; // these are outside the FDTD-domain as defined by the main disc
+
GetVV(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
GetVI(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
GetVV(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
diff --git a/FDTD/operator.h b/FDTD/operator.h
index eb056e7..463a168 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -31,6 +31,7 @@ class Operator
{
friend class Engine;
friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice
+ friend class Operator_Ext_PML_SF_Plane;
public:
//! Create a new operator
static Operator* New();
diff --git a/FDTD/operator_ext_pml_sf.cpp b/FDTD/operator_ext_pml_sf.cpp
new file mode 100644
index 0000000..7fbda01
--- /dev/null
+++ b/FDTD/operator_ext_pml_sf.cpp
@@ -0,0 +1,354 @@
+/*
+* 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 "tools/array_ops.h"
+
+bool Build_Split_Field_PML(Operator* op, int BC[6])
+{
+ cerr << "Build_Split_Field_PML:: Warning, currently only pml planes are implemented... edges and corner coming soon..." << endl;
+ 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);
+ op_pml_sf->SetPMLLength(8);
+ op_pml_sf->SetBoundaryCondition(BC);
+ 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;
+}
+
+Operator_Ext_PML_SF::~Operator_Ext_PML_SF()
+{
+ 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::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]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->GetMeshDelta(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::GetNodeArea(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;
+ if (ny==m_ny)
+ return m_Op->GetMeshDelta(m_nyP,l_pos,dualMesh)*m_Op->GetMeshDelta(m_nyPP,l_pos,dualMesh);
+ if (ny==m_nyP)
+ return m_Op->GetMeshDelta(m_ny,l_pos,dualMesh)*m_Op->GetMeshDelta(m_nyPP,l_pos,dualMesh);
+ if (ny==m_nyPP)
+ return m_Op->GetMeshDelta(m_nyP,l_pos,dualMesh)*m_Op->GetMeshDelta(m_ny,l_pos,dualMesh);
+ return 0.0;
+}
+
+double Operator_Ext_PML_SF_Plane::GetNodeLength(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->GetMeshDelta(ny,l_pos,dualMesh);
+}
+
+double Operator_Ext_PML_SF_Plane::GetKappaGraded(double depth) const
+{
+ if (depth<0)
+ return 0.0;
+
+ //todo: use fparser to allow arbitrary, user-defined profiles and parameter
+ double g = 2.5;
+ double R0 = 1e-6;
+ double kappa0 = -log(R0)*log(g)/(2*m_pml_delta * pow(g,m_pml_width/m_pml_delta) -1);
+ return pow(g,depth/m_pml_delta)*kappa0;
+}
+
+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);
+
+ UNUSED(nP);
+
+ double Zm = sqrt(inMat[2] / inMat[0]); // Zm = 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;
+ }
+ else
+ {
+ depth = m_pml_width - (pos[m_ny])*m_pml_delta;
+ kappa = GetKappaGraded(depth) / Zm ;
+ sigma = GetKappaGraded(depth-0.5*m_pml_delta) * Zm;
+ }
+// if ((pos[0]==0) && (pos[1]==0))
+// cerr << n << "\t" << pos[m_ny] << "\t"<< depth << "\t" << kappa << "\t" << sigma << endl;
+ }
+
+ inEC[0] = inMat[0] * GetNodeArea(n,pos) / GetNodeLength(n,pos);
+ inEC[1] = (inMat[1]+kappa) * GetNodeArea(n,pos) / GetNodeLength(n,pos);
+
+ inEC[2] = inMat[2] * GetNodeArea(n,pos,true) / GetNodeLength(n,pos,true);
+ inEC[3] = (inMat[3]+sigma) * GetNodeArea(n,pos,true) / GetNodeLength(n,pos,true);
+
+ 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].
+*/
+
+#ifndef OPERATOR_EXT_PML_SF_H
+#define OPERATOR_EXT_PML_SF_H
+
+#include "operator.h"
+#include "operator_extension.h"
+
+//! 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]);
+
+//! 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 GetNodeArea(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;}
+ virtual double GetNodeLength(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) const {UNUSED(depth);return 0.0;}
+
+ 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];
+
+ 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 GetNodeArea(int ny, unsigned int pos[3], bool dualMesh = false) const;
+ virtual double GetNodeLength(int ny, unsigned int pos[3], bool dualMesh = false) const;
+
+ virtual double GetKappaGraded(double depth) const;
+
+ virtual bool Calc_ECPos(int nP, int n, unsigned int* pos, double* inEC) const;
+
+ virtual Engine_Extension* CreateEngineExtention();
+
+ virtual bool IsCylinderCoordsSave() const {return false;}
+
+ 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/matlab/SetBoundaryCond.m b/matlab/SetBoundaryCond.m
index d26e834..1d69ca0 100644
--- a/matlab/SetBoundaryCond.m
+++ b/matlab/SetBoundaryCond.m
@@ -2,7 +2,7 @@ function FDTD = SetBoundaryCond(FDTD,BC)
% FDTD = SetBoundaryCond(FDTD,BC)
%
% BC = [xmin xmax ymin ymax zmin zmax];
-% ?min/?max: 0=PEC 1=PMC 2=MUR-ABC
+% ?min/?max: 0=PEC 1=PMC 2=MUR-ABC 3=PML-ABC
%
% openEMS matlab interface
% -----------------------
diff --git a/openEMS.pro b/openEMS.pro
index 5d1a2de..b451aed 100644
--- a/openEMS.pro
+++ b/openEMS.pro
@@ -77,7 +77,9 @@ SOURCES += main.cpp \
FDTD/operator_ext_dispersive.cpp \
FDTD/operator_ext_lorentzmaterial.cpp \
FDTD/engine_ext_dispersive.cpp \
- FDTD/engine_ext_lorentzmaterial.cpp
+ FDTD/engine_ext_lorentzmaterial.cpp \
+ FDTD/operator_ext_pml_sf.cpp \
+ FDTD/engine_ext_pml_sf.cpp
HEADERS += tools/ErrorMsg.h \
tools/AdrOp.h \
tools/constants.h \
@@ -111,7 +113,9 @@ HEADERS += tools/ErrorMsg.h \
FDTD/operator_ext_dispersive.h \
FDTD/operator_ext_lorentzmaterial.h \
FDTD/engine_ext_dispersive.h \
- FDTD/engine_ext_lorentzmaterial.h
+ FDTD/engine_ext_lorentzmaterial.h \
+ FDTD/operator_ext_pml_sf.h \
+ FDTD/engine_ext_pml_sf.h
QMAKE_CXXFLAGS_RELEASE = -O3 \
-g \
-march=native
diff --git a/openems.cpp b/openems.cpp
index a612bed..97cdd29 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -22,6 +22,7 @@
#include "FDTD/engine_multithread.h"
#include "FDTD/operator_multithread.h"
#include "FDTD/operator_ext_mur_abc.h"
+#include "FDTD/operator_ext_pml_sf.h"
#include "FDTD/processvoltage.h"
#include "FDTD/processcurrent.h"
#include "FDTD/process_efield.h"
@@ -266,13 +267,14 @@ int openEMS::SetupFDTD(const char* file)
//Mur-ABC, defined as extension to the operator
for (int n=0;n<6;++n)
{
- if (bounds[n]==2)
+ if (bounds[n]==2) //Mur-ABC
{
Operator_Ext_Mur_ABC* op_ext_mur = new Operator_Ext_Mur_ABC(FDTD_Op);
op_ext_mur->SetDirection(n/2,n%2);
FDTD_Op->AddExtension(op_ext_mur);
}
}
+ Build_Split_Field_PML(FDTD_Op,bounds);
if (CSX.GetQtyPropertyType(CSProperties::LORENTZMATERIAL)>0)
FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op));