extention updates & new extention: Mur's absorbing boundary condition

FDTD/engine.cpp

@@ -17,6 +17,7 @@
 
 #include "engine.h"
 #include "engine_extension.h"
+#include "operator_extension.h"
 #include "tools/array_ops.h"
 
 //! \brief construct an Engine instance
@@ -35,10 +36,25 @@ Engine::Engine(const Operator* op)
 	{
 		numLines[n] = Op->GetNumberOfLines(n);
 	}
+
+	for (size_t n=0;n<Op->GetNumberOfExtentions();++n)
+	{
+		Operator_Extension* op_ext = Op->GetExtension(n);
+		Engine_Extension* eng_ext = op_ext->CreateEngineExtention();
+		if (eng_ext)
+		{
+			eng_ext->SetEngine(this);
+			m_Eng_exts.push_back(eng_ext);
+		}
+	}
 }
 
 Engine::~Engine()
 {
+	for (size_t n=0;n<m_Eng_exts.size();++n)
+		delete m_Eng_exts.at(n);
+	m_Eng_exts.clear();
+
 	this->Reset();
 }
 
@@ -143,11 +159,6 @@ void Engine::ApplyCurrentExcite()
 {
 }
 
-void Engine::AddExtension(Engine_Extension* eng_ext)
-{
-	m_Eng_exts.push_back(eng_ext);
-}
-
 bool Engine::IterateTS(unsigned int iterTS)
 {
 	for (unsigned int iter=0;iter<iterTS;++iter)

FDTD/engine.h

@@ -40,7 +40,8 @@ public:
 	inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return volt[n][x][y][z]; }
 	inline virtual FDTD_FLOAT GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return curr[n][x][y][z]; }
 
-	virtual void AddExtension(Engine_Extension* eng_ext);
+	inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int pos[] ) { return volt[n][pos[0]][pos[1]][pos[2]]; }
+	inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int pos[] ) { return curr[n][pos[0]][pos[1]][pos[2]]; }
 
 protected:
 	Engine(const Operator* op);

FDTD/engine_ext_mur_abc.cpp

@@ -0,0 +1,105 @@
+/*
+*	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_mur_abc.h"
+#include "operator_ext_mur_abc.h"
+#include "engine.h"
+#include "tools/array_ops.h"
+
+Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext) : Engine_Extension(op_ext)
+{
+	m_Op_mur = op_ext;
+	m_numLines[0] = m_Op_mur->m_numLines[0];
+	m_numLines[1] = m_Op_mur->m_numLines[1];
+	m_ny = m_Op_mur->m_ny;
+	m_nyP = m_Op_mur->m_nyP;
+	m_nyPP = m_Op_mur->m_nyPP;
+	m_LineNr = m_Op_mur->m_LineNr;
+	m_LineNr_Shift = m_Op_mur->m_LineNr_Shift;
+
+	m_Mur_Coeff = m_Op_mur->m_Mur_Coeff;
+
+	m_volt_nyP = Create2DArray(m_numLines);
+	m_volt_nyPP = Create2DArray(m_numLines);
+}
+
+Engine_Ext_Mur_ABC::~Engine_Ext_Mur_ABC()
+{
+	Delete2DArray(m_volt_nyP,m_numLines);
+	m_volt_nyP = NULL;
+	Delete2DArray(m_volt_nyPP,m_numLines);
+	m_volt_nyPP = NULL;
+}
+
+void Engine_Ext_Mur_ABC::DoPreVoltageUpdates()
+{
+	if (m_Eng==NULL) return;
+	if (m_Mur_Coeff==0) return;
+	unsigned int pos[] = {0,0,0};
+	unsigned int pos_shift[] = {0,0,0};
+	pos[m_ny] = m_LineNr;
+	pos_shift[m_ny] = m_LineNr_Shift;
+
+	for (pos[m_nyP]=0;pos[m_nyP]<m_numLines[0];++pos[m_nyP])
+	{
+		pos_shift[m_nyP] = pos[m_nyP];
+		for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
+		{
+			pos_shift[m_nyPP] = pos[m_nyPP];
+			m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyP,pos_shift) - m_Mur_Coeff * m_Eng->GetVolt(m_nyP,pos);
+			m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyPP,pos_shift) - m_Mur_Coeff * m_Eng->GetVolt(m_nyPP,pos);
+		}
+	}
+}
+
+void Engine_Ext_Mur_ABC::DoPostVoltageUpdates()
+{
+	if (m_Eng==NULL) return;
+	if (m_Mur_Coeff==0) return;
+	unsigned int pos[] = {0,0,0};
+	unsigned int pos_shift[] = {0,0,0};
+	pos[m_ny] = m_LineNr;
+	pos_shift[m_ny] = m_LineNr_Shift;
+
+	for (pos[m_nyP]=0;pos[m_nyP]<m_numLines[0];++pos[m_nyP])
+	{
+		pos_shift[m_nyP] = pos[m_nyP];
+		for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
+		{
+			pos_shift[m_nyPP] = pos[m_nyPP];
+			m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Mur_Coeff * m_Eng->GetVolt(m_nyP,pos_shift);
+			m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Mur_Coeff * m_Eng->GetVolt(m_nyPP,pos_shift);
+		}
+	}
+}
+
+void Engine_Ext_Mur_ABC::Apply2Voltages()
+{
+	if (m_Eng==NULL) return;
+	if (m_Mur_Coeff==0) return;
+	unsigned int pos[] = {0,0,0};
+	pos[m_ny] = m_LineNr;
+
+	for (pos[m_nyP]=0;pos[m_nyP]<m_numLines[0];++pos[m_nyP])
+	{
+		for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
+		{
+			m_Eng->GetVolt(m_nyP,pos) = m_volt_nyP[pos[m_nyP]][pos[m_nyPP]];
+			m_Eng->GetVolt(m_nyPP,pos) = m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]];
+		}
+	}
+}

FDTD/engine_ext_mur_abc.h

@@ -0,0 +1,50 @@
+/*
+*	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_MUR_ABC_H
+#define ENGINE_EXT_MUR_ABC_H
+
+#include "engine_extension.h"
+#include "operator.h"
+
+class Operator_Ext_Mur_ABC;
+
+class Engine_Ext_Mur_ABC : public Engine_Extension
+{
+public:
+	Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext);
+	virtual ~Engine_Ext_Mur_ABC();
+
+	virtual void DoPreVoltageUpdates();
+	virtual void DoPostVoltageUpdates();
+	virtual void Apply2Voltages();
+
+protected:
+	Operator_Ext_Mur_ABC* m_Op_mur;
+
+	int m_ny;
+	int m_nyP,m_nyPP;
+	unsigned int m_LineNr;
+	int m_LineNr_Shift;
+	unsigned int m_numLines[2];
+
+	FDTD_FLOAT m_Mur_Coeff;
+	FDTD_FLOAT** m_volt_nyP; //n+1 direction
+	FDTD_FLOAT** m_volt_nyPP; //n+2 direction
+};
+
+#endif // ENGINE_EXT_MUR_ABC_H

FDTD/engine_extension.cpp

@@ -1,9 +1,26 @@
+/*
+*	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_extension.h"
 
 #include "engine.h"
 
-Engine_Extension::Engine_Extension(Operator_Extension* op_ext, Engine* eng)
+Engine_Extension::Engine_Extension(Operator_Extension* op_ext)
 {
 	m_Op_ext = op_ext;
-	m_Eng = eng;
+	m_Eng = NULL;
 }

FDTD/engine_extension.h

@@ -29,18 +29,20 @@ public:
 	virtual void DoPreVoltageUpdates() {}
 	//! This methode will be called __after__ the main engine does the usual voltage updates. This methode may __not__ change the engine voltages!!!
 	virtual void DoPostVoltageUpdates() {}
-	//! This methode will be called __after__ all updates to the voltages and extensions and may add its results to the engine voltages, but may __not__ rely on the current value of the engine voltages!!!
+	//! This methode will be called __after__ all updates to the voltages and extensions and may add/set its results to the engine voltages, but may __not__ rely on the current value of the engine voltages!!!
 	virtual void Apply2Voltages() {}
 
 	//! This methode will be called __before__ the main engine does the usual current updates. This methode may __not__ change the engine current!!!
 	virtual void DoPreCurrentUpdates() {}
 	//! This methode will be called __after__ the main engine does the usual current updates. This methode may __not__ change the engine current!!!
 	virtual void DoPostCurrentUpdates() {}
-	//! This methode will be called __after__ all updates to the current and extensions and may add its results to the engine current, but may __not__ rely on the current value of the engine current!!!
+	//! This methode will be called __after__ all updates to the current and extensions and may add/set its results to the engine current, but may __not__ rely on the current value of the engine current!!!
 	virtual void Apply2Current() {}
 
+	//! Set the Engine to this extention. This will usually done automatically by Engine::AddExtension
+	virtual void SetEngine(Engine* eng) {m_Eng=eng;}
 protected:
-	Engine_Extension(Operator_Extension* op_ext, Engine* eng);
+	Engine_Extension(Operator_Extension* op_ext);
 
 	Operator_Extension* m_Op_ext;
 	Engine* m_Eng;

FDTD/operator.cpp

@@ -35,6 +35,9 @@ Operator::Operator()
 
 Operator::~Operator()
 {
+	for (size_t n=0;n<m_Op_exts.size();++n)
+		delete m_Op_exts.at(n);
+	m_Op_exts.clear();
 	Reset();
 }
 

FDTD/operator.h

@@ -83,6 +83,8 @@ public:
 	virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL);
 
 	virtual void AddExtension(Operator_Extension* op_ext);
+	virtual size_t GetNumberOfExtentions() const {return m_Op_exts.size();}
+	virtual Operator_Extension* GetExtension(size_t index) const {return m_Op_exts.at(index);}
 
 protected:
 	//! use New() for creating a new Operator

FDTD/operator_ext_mur_abc.cpp

@@ -0,0 +1,80 @@
+/*
+*	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_mur_abc.h"
+#include "engine_ext_mur_abc.h"
+
+#include "tools/array_ops.h"
+
+Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op)
+{
+	m_Mur_Coeff = NULL;
+	m_ny = -1;
+	m_nyP = -1;
+	m_nyPP = -1;
+	m_LineNr = 0;
+	m_LineNr_Shift = 0;
+	m_Mur_Coeff = 0;
+}
+
+Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC()
+{
+}
+
+void Operator_Ext_Mur_ABC::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;
+	if (!top_ny)
+	{
+		m_LineNr = 0;
+		m_LineNr_Shift = 1;
+	}
+	else
+	{
+		m_LineNr = m_Op->GetNumberOfLines(m_ny)-1;
+		m_LineNr_Shift = m_Op->GetNumberOfLines(m_ny) - 2;
+	}
+
+	m_numLines[0] = m_Op->GetNumberOfLines(m_nyP);
+	m_numLines[1] = m_Op->GetNumberOfLines(m_nyPP);
+}
+
+bool Operator_Ext_Mur_ABC::BuildExtension()
+{
+	if (m_ny<0)
+	{
+		cerr << "Operator_Ext_Mur_ABC::BuildExtension: Warning, Extension not initialized! Use SetDirection!! Abort build!!" << endl;
+		return false;
+	}
+	double dT = m_Op->GetTimestep();
+	unsigned int pos[] = {0,0,0};
+	pos[m_ny] = m_LineNr;
+	double delta = fabs(m_Op->GetMeshDelta(m_ny,pos));
+	m_Mur_Coeff = (__C0__ * dT - delta) / (__C0__ * dT + delta);
+//	cerr << "Operator_Ext_Mur_ABC::BuildExtension(): " << m_Mur_Coeff << endl;
+	return true;
+}
+
+Engine_Extension* Operator_Ext_Mur_ABC::CreateEngineExtention()
+{
+	Engine_Ext_Mur_ABC* eng_ext = new Engine_Ext_Mur_ABC(this);
+	return eng_ext;
+}

FDTD/operator_ext_mur_abc.h

@@ -0,0 +1,49 @@
+/*
+*	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_MUR_ABC_H
+#define OPERATOR_EXT_MUR_ABC_H
+
+#include "operator.h"
+#include "operator_extension.h"
+
+class Operator_Ext_Mur_ABC : public Operator_Extension
+{
+	friend class Engine_Ext_Mur_ABC;
+public:
+	Operator_Ext_Mur_ABC(Operator* op);
+	~Operator_Ext_Mur_ABC();
+
+	//! Define the direction of this ABC: 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);
+
+	virtual bool BuildExtension();
+
+	virtual Engine_Extension* CreateEngineExtention();
+
+protected:
+	int m_ny;
+	int m_nyP,m_nyPP;
+	unsigned int m_LineNr;
+	int m_LineNr_Shift;
+
+	unsigned int m_numLines[2];
+
+	FDTD_FLOAT m_Mur_Coeff;
+};
+
+#endif // OPERATOR_EXT_MUR_ABC_H

FDTD/operator_extension.cpp

@@ -1,3 +1,20 @@
+/*
+*	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_extension.h"
 
 #include "operator.h"

FDTD/operator_extension.h

@@ -19,6 +19,7 @@
 #define OPERATOR_EXTENSION_H
 
 class Operator;
+class Engine_Extension;
 
 //! Abstract base-class for all operator extensions
 class Operator_Extension
@@ -26,6 +27,8 @@ class Operator_Extension
 public:
 	virtual bool BuildExtension() {return true;}
 
+	virtual Engine_Extension* CreateEngineExtention() {return 0;}
+
 protected:
 	Operator_Extension(Operator* op);
 	Operator* m_Op;

openEMS.pro

@@ -40,7 +40,9 @@ SOURCES += main.cpp \
     FDTD/operator_cylinder.cpp \
     FDTD/engine_cylinder.cpp \
     FDTD/operator_extension.cpp \
-    FDTD/engine_extension.cpp
+    FDTD/engine_extension.cpp \
+    FDTD/engine_ext_mur_abc.cpp \
+    FDTD/operator_ext_mur_abc.cpp
 HEADERS += tools/ErrorMsg.h \
     tools/AdrOp.h \
     tools/constants.h \
@@ -58,7 +60,9 @@ HEADERS += tools/ErrorMsg.h \
     FDTD/operator_cylinder.h \
     FDTD/engine_cylinder.h \
     FDTD/operator_extension.h \
-    FDTD/engine_extension.h
+    FDTD/engine_extension.h \
+    FDTD/engine_ext_mur_abc.h \
+    FDTD/operator_ext_mur_abc.h
 QMAKE_CXXFLAGS_RELEASE = -O2 \
     -g \
     -march=native

openems.cpp

@@ -21,6 +21,7 @@
 #include "FDTD/engine.h"
 #include "FDTD/engine_cylinder.h"
 #include "FDTD/engine_multithread.h"
+#include "FDTD/operator_ext_mur_abc.h"
 #include "FDTD/processvoltage.h"
 #include "FDTD/processcurrent.h"
 #include "FDTD/processfields_td.h"
@@ -254,6 +255,18 @@ int openEMS::SetupFDTD(const char* file)
 
 	if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(2);
 
+	/**************************** create all operator/engine extensions here !!!! **********************************/
+	//Mur-ABC
+	for (int n=0;n<6;++n)
+	{
+		if (bounds[n]==2)
+		{
+			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);
+		}
+	}
+
 	FDTD_Op->CalcECOperator();
 
 	SetupExcitation(FDTD_Opts->FirstChildElement("Excitation"));

tools/array_ops.cpp

@@ -18,6 +18,33 @@
 #include "array_ops.h"
 #include <ostream>
 
+FDTD_FLOAT** Create2DArray(const unsigned int* numLines)
+{
+	FDTD_FLOAT** array=NULL;
+	unsigned int pos[3];
+	array = new FDTD_FLOAT*[numLines[0]];
+	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+	{
+		array[pos[0]] = new FDTD_FLOAT[numLines[1]];
+		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+		{
+			array[pos[0]][pos[1]] = 0;
+		}
+	}
+	return array;
+}
+
+void Delete2DArray(FDTD_FLOAT** array, const unsigned int* numLines)
+{
+	if (array==NULL) return;
+	unsigned int pos[3];
+	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+	{
+		delete[] array[pos[0]];
+	}
+	delete[] array;
+}
+
 FDTD_FLOAT*** Create3DArray(const unsigned int* numLines)
 {
 	FDTD_FLOAT*** array=NULL;

tools/array_ops.h

@@ -20,6 +20,9 @@
 
 #include "../FDTD/operator.h"
 
+FDTD_FLOAT** Create2DArray(const unsigned int* numLines);
+void Delete2DArray(FDTD_FLOAT** array, const unsigned int* numLines);
+
 FDTD_FLOAT*** Create3DArray(const unsigned int* numLines);
 void Delete3DArray(FDTD_FLOAT*** array, const unsigned int* numLines);