Merge branch with multigrid

Conflicts:
	FDTD/operator.h
	FDTD/operator_cylinder.cpp

FDTD/engine.cpp

@@ -48,7 +48,6 @@ Engine::~Engine()
 
 void Engine::Init()
 {
-	Reset();
 	numTS = 0;
 	volt = Create_N_3DArray<FDTD_FLOAT>(numLines);
 	curr = Create_N_3DArray<FDTD_FLOAT>(numLines);
@@ -73,6 +72,13 @@ void Engine::InitExtensions()
 	}
 }
 
+void Engine::ClearExtensions()
+{
+	for (size_t n=0;n<m_Eng_exts.size();++n)
+		delete m_Eng_exts.at(n);
+	m_Eng_exts.clear();
+}
+
 void Engine::Reset()
 {
 	Delete_N_3DArray(volt,numLines);
@@ -83,9 +89,7 @@ void Engine::Reset()
 	file_et.close();
 	file_ht.close();
 
-	for (size_t n=0;n<m_Eng_exts.size();++n)
+	ClearExtensions();
-		delete m_Eng_exts.at(n);
-	m_Eng_exts.clear();
 }
 
 void Engine::UpdateVoltages(unsigned int startX, unsigned int numX)
@@ -125,13 +129,19 @@ void Engine::UpdateVoltages(unsigned int startX, unsigned int numX)
 void Engine::ApplyVoltageExcite()
 {
 	int exc_pos;
+	unsigned int ny;
+	unsigned int pos[3];
 	//soft voltage excitation here (E-field excite)
 	for (unsigned int n=0;n<Op->Exc->Volt_Count;++n)
 	{
 		exc_pos = (int)numTS - (int)Op->Exc->Volt_delay[n];
 		exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
 //			if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl;
-		GetVolt(Op->Exc->Volt_dir[n],Op->Exc->Volt_index[0][n],Op->Exc->Volt_index[1][n],Op->Exc->Volt_index[2][n]) += Op->Exc->Volt_amp[n]*Op->Exc->Signal_volt[exc_pos];
+		ny = Op->Exc->Volt_dir[n];
+		pos[0]=Op->Exc->Volt_index[0][n];
+		pos[1]=Op->Exc->Volt_index[1][n];
+		pos[2]=Op->Exc->Volt_index[2][n];
+		SetVolt(ny,pos, GetVolt(ny,pos) + Op->Exc->Volt_amp[n]*Op->Exc->Signal_volt[exc_pos]);
 	}
 
 	// write the first excitation into the file "et"
@@ -172,13 +182,19 @@ void Engine::UpdateCurrents(unsigned int startX, unsigned int numX)
 void Engine::ApplyCurrentExcite()
 {
 	int exc_pos;
+	unsigned int ny;
+	unsigned int pos[3];
 	//soft current excitation here (H-field excite)
 	for (unsigned int n=0;n<Op->Exc->Curr_Count;++n)
 	{
 		exc_pos = (int)numTS - (int)Op->Exc->Curr_delay[n];
 		exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
 //			if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl;
-		GetCurr(Op->Exc->Curr_dir[n],Op->Exc->Curr_index[0][n],Op->Exc->Curr_index[1][n],Op->Exc->Curr_index[2][n]) += Op->Exc->Curr_amp[n]*Op->Exc->Signal_curr[exc_pos];
+		ny = Op->Exc->Curr_dir[n];
+		pos[0]=Op->Exc->Curr_index[0][n];
+		pos[1]=Op->Exc->Curr_index[1][n];
+		pos[2]=Op->Exc->Curr_index[2][n];
+		SetCurr(ny,pos, GetCurr(ny,pos) + Op->Exc->Curr_amp[n]*Op->Exc->Signal_curr[exc_pos]);
 	}
 
 	// write the first excitation into the file "ht"

FDTD/engine.h

@@ -47,10 +47,15 @@ public:
 	virtual unsigned int GetNumberOfTimesteps() {return numTS;};
 
 	//this access functions muss be overloaded by any new engine using a different storage model
-	inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return volt[n][x][y][z]; }
+	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& GetVolt( unsigned int n, unsigned int pos[3] ) { return volt[n][pos[0]][pos[1]][pos[2]]; }
+	inline virtual FDTD_FLOAT GetVolt( unsigned int n, const unsigned int pos[3] )							const { return volt[n][pos[0]][pos[1]][pos[2]]; }
-	inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return curr[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]; }
-	inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int pos[3] ) { return curr[n][pos[0]][pos[1]][pos[2]]; }
+	inline virtual FDTD_FLOAT GetCurr( unsigned int n, const unsigned int pos[3] )							const { return curr[n][pos[0]][pos[1]][pos[2]]; }
+
+	inline virtual void SetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value)	{ volt[n][x][y][z]=value; }
+	inline virtual void SetVolt( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value )						{ volt[n][pos[0]][pos[1]][pos[2]]=value; }
+	inline virtual void SetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value)	{ curr[n][x][y][z]=value; }
+	inline virtual void SetCurr( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value )						{ curr[n][pos[0]][pos[1]][pos[2]]=value; }
 
 	virtual void UpdateVoltages(unsigned int startX, unsigned int numX);
 	virtual void ApplyVoltageExcite();
@@ -75,6 +80,7 @@ protected:
 	unsigned int numTS;
 
 	virtual void InitExtensions();
+	virtual void ClearExtensions();
 	vector<Engine_Extension*> m_Eng_exts;
 
 	ofstream file_et, file_ht; //excite signal dump file

FDTD/engine_cylindermultigrid.cpp

@@ -0,0 +1,260 @@
+/*
+*	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_cylindermultigrid.h"
+#include "operator_cylindermultigrid.h"
+#include "engine_ext_cylindermultigrid.h"
+
+Engine_CylinderMultiGrid* Engine_CylinderMultiGrid::New(const Operator_CylinderMultiGrid* op, unsigned int numThreads)
+{
+	cout << "Create FDTD engine (cylindrical multi grid mesh using sse compression + multithreading)" << endl;
+	Engine_CylinderMultiGrid* e = new Engine_CylinderMultiGrid(op);
+	e->setNumThreads( numThreads );
+	e->Init();
+	return e;
+}
+
+Engine_CylinderMultiGrid::Engine_CylinderMultiGrid(const Operator_CylinderMultiGrid* op) : Engine_Multithread(op)
+{
+	Op_CMG = op;
+
+	m_WaitOnBase = new boost::barrier(2);
+	m_WaitOnChild = new boost::barrier(2);
+	m_WaitOnSync = new boost::barrier(2);
+
+	m_Eng_Ext_MG = new Engine_Ext_CylinderMultiGrid(NULL,true);
+	m_Eng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync);
+	m_Eng_Ext_MG->SetEngine(this);
+
+	Engine* eng = op->GetInnerOperator()->CreateEngine();
+	m_InnerEngine = dynamic_cast<Engine_Multithread*>(eng);
+
+	Engine_Ext_CylinderMultiGrid* m_InnerEng_Ext_MG = new Engine_Ext_CylinderMultiGrid(NULL,false);
+	m_InnerEng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync);
+
+	// if already has a base extension, switch places ... seems to be faster...
+	for (size_t n=0;n<m_InnerEngine->m_Eng_exts.size();++n)
+	{
+		Engine_Ext_CylinderMultiGrid* eng_mg = dynamic_cast<Engine_Ext_CylinderMultiGrid*>(m_InnerEngine->m_Eng_exts.at(n));
+		if (eng_mg)
+		{
+			m_InnerEngine->m_Eng_exts.at(n) = m_InnerEng_Ext_MG;
+			m_InnerEng_Ext_MG = eng_mg;
+			break;
+		}
+	}
+	m_InnerEngine->m_Eng_exts.push_back(m_InnerEng_Ext_MG);
+}
+
+Engine_CylinderMultiGrid::~Engine_CylinderMultiGrid()
+{
+	m_Thread_NumTS = 0;
+	m_startBarrier->wait();
+
+	m_IteratorThread_Group.join_all();
+
+	delete m_InnerEngine;
+	m_InnerEngine = NULL;
+
+	delete m_WaitOnBase;
+	m_WaitOnBase = NULL;
+	delete m_WaitOnChild;
+	m_WaitOnChild = NULL;
+	delete m_WaitOnSync;
+	m_WaitOnSync = NULL;
+
+	delete m_startBarrier;
+	m_startBarrier = NULL;
+	delete m_stopBarrier;
+	m_stopBarrier = NULL;
+
+}
+
+void Engine_CylinderMultiGrid::Init()
+{
+	Engine_Multithread::Init();
+
+	m_Eng_exts.push_back(m_Eng_Ext_MG);
+
+	m_startBarrier = new boost::barrier(3); //both engines + organizer
+	m_stopBarrier = new boost::barrier(3); //both engines + organizer
+
+	boost::thread *t = NULL;
+
+	t = new boost::thread( Engine_CylinderMultiGrid_Thread(this,m_startBarrier,m_stopBarrier,&m_Thread_NumTS, true) );
+	m_IteratorThread_Group.add_thread( t );
+
+	t = new boost::thread( Engine_CylinderMultiGrid_Thread(m_InnerEngine,m_startBarrier,m_stopBarrier,&m_Thread_NumTS, false) );
+	m_IteratorThread_Group.add_thread( t );
+}
+
+bool Engine_CylinderMultiGrid::IterateTS(unsigned int iterTS)
+{
+	m_Thread_NumTS = iterTS;
+
+	m_startBarrier->wait(); //start base and child iterations
+
+	m_stopBarrier->wait();  //tell base and child to wait for another start event...
+
+	//interpolate child data to base mesh...
+	for (unsigned int n=0;n<Op_CMG->m_Split_Pos-1;++n)
+		InterpolVoltChild2Base(n);
+	for (unsigned int n=0;n<Op_CMG->m_Split_Pos-2;++n)
+		InterpolCurrChild2Base(n);
+
+	return true;
+}
+
+void Engine_CylinderMultiGrid::InterpolVoltChild2Base(unsigned int rzPlane)
+{
+	//interpolate voltages from child engine to the base engine...
+	unsigned int pos[3];
+	unsigned int pos1_rz1, pos1_rz2;
+	unsigned int pos1_a1, pos1_a2;
+	pos[0] = rzPlane;
+	bool isOdd, isEven;
+	f4vector half, one_eighth, three_eighth;
+	for (int n=0;n<4;++n)
+	{
+		half.f[n]=0.5;
+		one_eighth.f[n]   = 1.0/8.0;
+		three_eighth.f[n] = 3.0/8.0;
+	}
+	for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+	{
+		isOdd = (pos[1]%2);
+		isEven = !isOdd * (pos[1]>0);
+
+		pos1_rz1 = pos[1]/2;
+		pos1_rz2 = pos[1]/2 + isOdd;
+
+		pos1_rz1 += pos[1]==1;
+		pos1_rz2 -= pos[1]==(numLines[1]-2);
+
+		pos1_a1 = pos[1]/2;
+		pos1_a2 = pos[1]/2 + isOdd - isEven;
+		pos1_a2 += pos[1]==(numLines[1]-1);
+		pos1_a2 -= pos[1]==(numLines[1]-2);
+
+		for (pos[2]=0; pos[2]<numVectors; ++pos[2])
+		{
+			//r - direction
+			f4_volt[0][pos[0]][pos[1]][pos[2]].v = half.v * ( m_InnerEngine->f4_volt[0][pos[0]][pos1_rz1][pos[2]].v + m_InnerEngine->f4_volt[0][pos[0]][pos1_rz2][pos[2]].v);
+
+			//z - direction
+			f4_volt[2][pos[0]][pos[1]][pos[2]].v = half.v * ( m_InnerEngine->f4_volt[2][pos[0]][pos1_rz1][pos[2]].v + m_InnerEngine->f4_volt[2][pos[0]][pos1_rz2][pos[2]].v);
+
+			//alpha - direction
+			f4_volt[1][pos[0]][pos[1]][pos[2]].v  = three_eighth.v * m_InnerEngine->f4_volt[1][pos[0]][pos1_a1][pos[2]].v;
+			f4_volt[1][pos[0]][pos[1]][pos[2]].v += one_eighth.v * m_InnerEngine->f4_volt[1][pos[0]][pos1_a2][pos[2]].v;
+		}
+	}
+
+	//r,z - interpolation correction...
+	pos[1]=1;
+	for (pos[2]=0; pos[2]<numVectors; ++pos[2])
+	{
+		f4_volt[0][pos[0]][pos[1]][pos[2]].v  += half.v * (m_InnerEngine->f4_volt[0][pos[0]][1][pos[2]].v - m_InnerEngine->f4_volt[0][pos[0]][2][pos[2]].v);
+		f4_volt[2][pos[0]][pos[1]][pos[2]].v  += half.v * (m_InnerEngine->f4_volt[2][pos[0]][1][pos[2]].v - m_InnerEngine->f4_volt[2][pos[0]][2][pos[2]].v);
+	}
+	pos[1]=numLines[1]-2;
+	for (pos[2]=0; pos[2]<numVectors; ++pos[2])
+	{
+		f4_volt[0][pos[0]][pos[1]][pos[2]].v  += half.v * (m_InnerEngine->f4_volt[0][pos[0]][pos[1]/2][pos[2]].v - m_InnerEngine->f4_volt[0][pos[0]][pos[1]/2-1][pos[2]].v);
+		f4_volt[2][pos[0]][pos[1]][pos[2]].v  += half.v * (m_InnerEngine->f4_volt[2][pos[0]][pos[1]/2][pos[2]].v - m_InnerEngine->f4_volt[2][pos[0]][pos[1]/2-1][pos[2]].v);
+	}
+}
+
+void Engine_CylinderMultiGrid::InterpolCurrChild2Base(unsigned int rzPlane)
+{
+	//interpolate voltages from child engine to the base engine...
+	unsigned int pos[3];
+	unsigned int pos1_rz1, pos1_rz2;
+	unsigned int pos1_a1, pos1_a2;
+	pos[0] = rzPlane;
+	bool isOdd, isEven;
+	f4vector quarter, one_fourth, three_fourth;
+	for (int n=0;n<4;++n)
+	{
+		quarter.f[n]=0.25;
+		one_fourth.f[n] = 1.0/4.0;
+		three_fourth.f[n] = 3.0/4.0;
+	}
+	for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+	{
+		isOdd = (pos[1]%2);
+		isEven = !isOdd * (pos[1]>0);
+
+		pos1_a1 = pos[1]/2;
+		pos1_a2 = pos[1]/2 + isOdd;
+
+		pos1_a1 += pos[1]==1;
+		pos1_a2 -= pos[1]==(numLines[1]-2);
+
+		pos1_rz1 = pos[1]/2;
+		pos1_rz2 = pos[1]/2 + isOdd - isEven;
+		pos1_rz2 += pos[1]==(numLines[1]-1);
+		pos1_rz2 -= pos[1]==(numLines[1]-2);
+
+		for (pos[2]=0; pos[2]<numVectors; ++pos[2])
+		{
+			//r - direction
+			f4_curr[0][pos[0]][pos[1]][pos[2]].v  = three_fourth.v * m_InnerEngine->f4_curr[0][pos[0]][pos1_rz1][pos[2]].v;
+			f4_curr[0][pos[0]][pos[1]][pos[2]].v += one_fourth.v * m_InnerEngine->f4_curr[0][pos[0]][pos1_rz2][pos[2]].v;
+
+			//z - direction
+			f4_curr[2][pos[0]][pos[1]][pos[2]].v  = three_fourth.v * m_InnerEngine->f4_curr[2][pos[0]][pos1_rz1][pos[2]].v;
+			f4_curr[2][pos[0]][pos[1]][pos[2]].v += one_fourth.v * m_InnerEngine->f4_curr[2][pos[0]][pos1_rz2][pos[2]].v;
+
+			//alpha - direction
+			f4_curr[1][pos[0]][pos[1]][pos[2]].v  = quarter.v * (m_InnerEngine->f4_curr[1][pos[0]][pos1_a1][pos[2]].v + m_InnerEngine->f4_curr[1][pos[0]][pos1_a2][pos[2]].v);
+		}
+
+	}
+	//alpha - interpolation correction...
+	pos[1]=1;
+	for (pos[2]=0; pos[2]<numVectors; ++pos[2])
+		f4_curr[1][pos[0]][pos[1]][pos[2]].v  += quarter.v * (m_InnerEngine->f4_curr[1][pos[0]][1][pos[2]].v - m_InnerEngine->f4_curr[1][pos[0]][2][pos[2]].v);
+	pos[1]=numLines[1]-2;
+	for (pos[2]=0; pos[2]<numVectors; ++pos[2])
+		f4_curr[1][pos[0]][pos[1]][pos[2]].v  += quarter.v * (m_InnerEngine->f4_curr[1][pos[0]][pos[1]/2][pos[2]].v - m_InnerEngine->f4_curr[1][pos[0]][pos[1]/2-1][pos[2]].v);
+}
+
+/****************************************************************************************/
+Engine_CylinderMultiGrid_Thread::Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS, bool isBase)
+{
+	m_startBarrier = start;
+	m_stopBarrier = stop;
+	m_Eng=engine;
+	m_isBase=isBase;
+	m_numTS = numTS;
+}
+
+void Engine_CylinderMultiGrid_Thread::operator()()
+{
+	m_startBarrier->wait(); //wait for Base engine to start the iterations...
+
+	while(*m_numTS>0)	//m_numTS==0 request to terminate this thread...
+	{
+		if (m_isBase)
+			m_Eng->Engine_Multithread::IterateTS(*m_numTS);
+		else
+			m_Eng->IterateTS(*m_numTS);
+		m_stopBarrier->wait();		//sync all workers after iterations are performed
+		m_startBarrier->wait();		//wait for Base engine to start the iterations again ...
+	}
+}

FDTD/engine_cylindermultigrid.h

@@ -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/>.
+*/
+
+#ifndef ENGINE_CYLINDERMULTIGRID_H
+#define ENGINE_CYLINDERMULTIGRID_H
+
+#include "engine_multithread.h"
+
+class Operator_CylinderMultiGrid;
+class Engine_CylinderMultiGrid_Thread;
+class Engine_Ext_CylinderMultiGrid;
+
+class Engine_CylinderMultiGrid : public Engine_Multithread
+{
+	friend class Engine_Ext_CylinderMultiGrid;
+public:
+    Engine_CylinderMultiGrid();
+
+	static Engine_CylinderMultiGrid* New(const Operator_CylinderMultiGrid* op, unsigned int numThreads = 0);
+	virtual ~Engine_CylinderMultiGrid();
+
+	virtual void InterpolVoltChild2Base(unsigned int rzPlane);
+	virtual void InterpolCurrChild2Base(unsigned int rzPlane);
+
+	virtual void Init();
+
+	//! Iterate \a iterTS number of timesteps
+	virtual bool IterateTS(unsigned int iterTS);
+
+protected:
+	Engine_CylinderMultiGrid(const Operator_CylinderMultiGrid* op);
+	const Operator_CylinderMultiGrid* Op_CMG;
+
+	Engine_Multithread* m_InnerEngine;
+
+	volatile unsigned int m_Thread_NumTS;
+	boost::thread_group m_IteratorThread_Group;
+	boost::barrier *m_startBarrier;
+	boost::barrier *m_stopBarrier;
+	Engine_CylinderMultiGrid_Thread* m_IteratorThread;
+	Engine_CylinderMultiGrid_Thread* m_InnerIteratorThread;
+	
+	//extension barrier
+	boost::barrier *m_WaitOnBase;
+	boost::barrier *m_WaitOnChild;
+	boost::barrier *m_WaitOnSync;
+
+	Engine_Ext_CylinderMultiGrid* m_Eng_Ext_MG;
+};
+
+
+class Engine_CylinderMultiGrid_Thread
+{
+public:
+	Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS, bool isBase);
+	void operator()();
+
+protected:
+	Engine_Multithread *m_Eng;
+	bool m_isBase;
+	boost::barrier *m_startBarrier;
+	boost::barrier *m_stopBarrier;
+	volatile unsigned int *m_numTS;
+};
+
+#endif // ENGINE_CYLINDERMULTIGRID_H

FDTD/engine_ext_cylinder.cpp

@@ -40,18 +40,18 @@ void Engine_Ext_Cylinder::Apply2Voltages()
 		pos[0] = 0;
 		for (pos[2]=0;pos[2]<numLines[2];++pos[2])
 		{
-			m_Eng->GetVolt(2,0,0,pos[2]) *= cyl_Op->vv_R0[pos[2]];
+			m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2])*cyl_Op->vv_R0[pos[2]]);
 			for (pos[1]=0;pos[1]<numLines[1]-1;++pos[1])
 			{
-					m_Eng->GetVolt(2,0,0,pos[2]) += cyl_Op->vi_R0[pos[2]] *  m_Eng->GetCurr(1,0,pos[1],pos[2]);
+					m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2]) + cyl_Op->vi_R0[pos[2]] *  m_Eng->GetCurr(1,0,pos[1],pos[2]) );
 			}
 		}
 		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
 		{
 			for (pos[2]=0;pos[2]<numLines[2];++pos[2])
 			{
-				m_Eng->GetVolt(1,0,pos[1],pos[2]) = 0; //no voltage in alpha-direction at r=0
+				m_Eng->SetVolt(1,0,pos[1],pos[2], 0); //no voltage in alpha-direction at r=0
-				m_Eng->GetVolt(2,0,pos[1],pos[2]) = m_Eng->GetVolt(2,0,0,pos[2]);
+				m_Eng->SetVolt(2,0,pos[1],pos[2], m_Eng->GetVolt(2,0,0,pos[2]) );
 			}
 		}
 	}
@@ -64,9 +64,9 @@ void Engine_Ext_Cylinder::Apply2Voltages()
 	{
 		for (pos[2]=0;pos[2]<numLines[2];++pos[2])
 		{
-			m_Eng->GetVolt(0,pos[0],0,pos[2]) = m_Eng->GetVolt(0,pos[0],last_A_Line,pos[2]);
+			m_Eng->SetVolt(0,pos[0],0,pos[2], m_Eng->GetVolt(0,pos[0],last_A_Line,pos[2]) );
-			m_Eng->GetVolt(1,pos[0],0,pos[2]) = m_Eng->GetVolt(1,pos[0],last_A_Line,pos[2]);
+			m_Eng->SetVolt(1,pos[0],0,pos[2], m_Eng->GetVolt(1,pos[0],last_A_Line,pos[2]) );
-			m_Eng->GetVolt(2,pos[0],0,pos[2]) = m_Eng->GetVolt(2,pos[0],last_A_Line,pos[2]);
+			m_Eng->SetVolt(2,pos[0],0,pos[2], m_Eng->GetVolt(2,pos[0],last_A_Line,pos[2]) );
 		}
 	}
 }
@@ -83,9 +83,9 @@ void Engine_Ext_Cylinder::Apply2Current()
 		unsigned int last_A_Line = numLines[1]-1;
 		for (pos[2]=0;pos[2]<numLines[2]-1;++pos[2])
 		{
-			m_Eng->GetCurr(0,pos[0],last_A_Line,pos[2]) = m_Eng->GetCurr(0,pos[0],0,pos[2]);
+			m_Eng->SetCurr(0,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(0,pos[0],0,pos[2]) );
-			m_Eng->GetCurr(1,pos[0],last_A_Line,pos[2]) = m_Eng->GetCurr(1,pos[0],0,pos[2]);
+			m_Eng->SetCurr(1,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(1,pos[0],0,pos[2]) );
-			m_Eng->GetCurr(2,pos[0],last_A_Line,pos[2]) = m_Eng->GetCurr(2,pos[0],0,pos[2]);
+			m_Eng->SetCurr(2,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(2,pos[0],0,pos[2]) );
 		}
 	}
 }

FDTD/engine_ext_cylindermultigrid.cpp

@@ -0,0 +1,158 @@
+/*
+*	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.h"
+#include "engine_ext_cylindermultigrid.h"
+#include "engine_cylindermultigrid.h"
+
+Engine_Ext_CylinderMultiGrid::Engine_Ext_CylinderMultiGrid(Operator_Extension* op_ext, bool isBase) : Engine_Extension(op_ext)
+{
+	m_IsBase = isBase;
+	m_Eng_MG = NULL;
+}
+
+Engine_Ext_CylinderMultiGrid::~Engine_Ext_CylinderMultiGrid()
+{
+}
+
+void Engine_Ext_CylinderMultiGrid::SetBarrier(boost::barrier* waitBase, boost::barrier* waitChild, boost::barrier* waitSync)
+{
+	m_WaitOnBase = waitBase;
+	m_WaitOnChild = waitChild;
+	m_WaitOnSync = waitSync;
+}
+
+void Engine_Ext_CylinderMultiGrid::SetEngine(Engine* eng)
+{
+	m_Eng_MG = dynamic_cast<Engine_CylinderMultiGrid*>(eng);
+	if (m_Eng_MG==NULL)
+	{
+		cerr << "Engine_Ext_CylinderMultiGrid::SetEngine(): Error" << endl;
+		exit(0);
+	}
+}
+
+void Engine_Ext_CylinderMultiGrid::DoPreVoltageUpdates()
+{
+	//cerr << "Engine_Ext_CylinderMultiGrid::DoPreVoltageUpdates() for " << m_IsBase << endl;
+	if (!m_IsBase)
+	{
+		//cerr << "child: volt wait on base " << endl;
+		m_WaitOnBase->wait(); //wait on base to finisch current sync and/or to finisch voltage updates, than start child voltage updates
+	}
+}
+
+void Engine_Ext_CylinderMultiGrid::DoPostVoltageUpdates()
+{
+
+}
+
+void Engine_Ext_CylinderMultiGrid::Apply2Voltages()
+{
+	if (m_IsBase)
+	{
+		m_WaitOnBase->wait();	//base voltage updates are done, tell child to start its voltage updates
+		m_WaitOnChild->wait();	//wait for child to finisch its updates
+		SyncVoltages();			//child is finisch, run sync and go to current updates next
+		m_WaitOnSync->wait();	//sync is done... move on and tell child to move on...
+	}
+	else
+	{
+		m_WaitOnChild->wait();	//child is finished voltage updates, will tell base to run sync
+		m_WaitOnSync->wait();	//wait for base to finisch sync before going to wait for current updates
+	}
+}
+
+void Engine_Ext_CylinderMultiGrid::SyncVoltages()
+{
+	if (m_Eng_MG==NULL)
+	{
+		cerr << "Engine_Ext_CylinderMultiGrid::SyncVoltages: Error engine is NULL" << endl;
+		return;
+	}
+
+	unsigned int* numLines = m_Eng_MG->numLines;
+
+	Engine_Multithread* m_InnerEng = m_Eng_MG->m_InnerEngine;
+
+	//interpolate voltages from base engine to child engine...
+	unsigned int pos[3];
+	pos[0] = m_Eng_MG->Op_CMG->GetSplitPos()-1;
+	unsigned int pos1_half = 0;
+	f4vector v_null; v_null.f[0] = 0;v_null.f[1] = 0;v_null.f[2] = 0;v_null.f[3] = 0;
+	for (pos[1]=0; pos[1]<numLines[1]-1; pos[1]+=2)
+	{
+		pos1_half = pos[1]/2;
+		for (pos[2]=0; pos[2]<m_Eng_MG->numVectors; ++pos[2])
+		{
+			//r - direczion
+			m_InnerEng->f4_volt[0][pos[0]][pos1_half][pos[2]].v = v_null.v;
+
+			//z - direction
+			m_InnerEng->f4_volt[2][pos[0]][pos1_half][pos[2]].v = m_Eng_MG->f4_volt[2][pos[0]][pos[1]][pos[2]].v;
+
+			//alpha - direction
+			m_InnerEng->f4_volt[1][pos[0]][pos1_half][pos[2]].v  = m_Eng_MG->f4_volt[1][pos[0]][pos[1]][pos[2]].v;
+			m_InnerEng->f4_volt[1][pos[0]][pos1_half][pos[2]].v += m_Eng_MG->f4_volt[1][pos[0]][pos[1]+1][pos[2]].v;
+		}
+	}
+}
+
+void Engine_Ext_CylinderMultiGrid::DoPreCurrentUpdates()
+{
+	//cerr << "Engine_Ext_CylinderMultiGrid::DoPreCurrentUpdates() for " << m_IsBase << endl;
+	if (!m_IsBase)
+	{
+		//cerr << "child: curr wait on base " << endl;
+		m_WaitOnBase->wait(); //wait on base to finisch voltage sync and current updates, than start child current updates
+	}
+}
+
+void Engine_Ext_CylinderMultiGrid::DoPostCurrentUpdates()
+{
+
+}
+
+void Engine_Ext_CylinderMultiGrid::Apply2Current()
+{
+	if (m_IsBase)
+	{
+		//cerr << "Base: curr wait on base done, wait on sync" << endl;
+		m_WaitOnBase->wait();	//base current updates are done, tell child to start its current updates
+		m_WaitOnChild->wait();	//wait for child to finisch its updates
+		SyncCurrents();			//child is finisch, run sync and go to voltage updates next
+		m_WaitOnSync->wait();	//sync is done... move on and tell child to move on...
+	}
+	else
+	{
+		m_WaitOnChild->wait();	//child is finished current updates, will tell base to run sync...
+		m_WaitOnSync->wait();	//wait for base to finisch sync before going to wait for next voltage updates
+		//cerr << "Child: curr done, wait on sync" << endl;
+	}
+}
+
+void Engine_Ext_CylinderMultiGrid::SyncCurrents()
+{
+	if (m_Eng_MG==NULL)
+	{
+		cerr << "Engine_Ext_CylinderMultiGrid::SyncCurrents: Error engine is NULL" << endl;
+		return;
+	}
+
+	m_Eng_MG->InterpolCurrChild2Base(m_Eng_MG->Op_CMG->GetSplitPos()-2);
+	return;
+}

FDTD/engine_ext_cylindermultigrid.h

@@ -0,0 +1,58 @@
+/*
+*	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_CYLINDERMULTIGRID_H
+#define ENGINE_EXT_CYLINDERMULTIGRID_H
+
+#include "engine_cylindermultigrid.h"
+#include "engine_extension.h"
+#include "operator_cylindermultigrid.h"
+
+class Operator_Ext_CylinderMultiGrid;
+
+class Engine_Ext_CylinderMultiGrid : public Engine_Extension
+{
+public:
+	Engine_Ext_CylinderMultiGrid(Operator_Extension* op_ext, bool isBase);
+	virtual ~Engine_Ext_CylinderMultiGrid();
+
+	void SetBarrier(boost::barrier* waitBase, boost::barrier* waitChild, boost::barrier* waitSync);
+
+	virtual void DoPreVoltageUpdates();
+	virtual void DoPostVoltageUpdates();
+	virtual void Apply2Voltages();
+
+	virtual void DoPreCurrentUpdates();
+	virtual void DoPostCurrentUpdates();
+	virtual void Apply2Current();
+
+	virtual void SetEngine(Engine* eng);
+
+protected:
+	void SyncVoltages();
+	void SyncCurrents();
+
+	Engine_CylinderMultiGrid* m_Eng_MG;
+
+	boost::barrier *m_WaitOnBase;
+	boost::barrier *m_WaitOnChild;
+	boost::barrier *m_WaitOnSync;
+
+	bool m_IsBase;
+};
+
+#endif // ENGINE_EXT_CYLINDERMULTIGRID_H

FDTD/engine_ext_dispersive.cpp

@@ -68,9 +68,9 @@ void Engine_Ext_Dispersive::Apply2Voltages()
 		{
 			for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i)
 			{
-				m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[0][i];
+				m_Eng->Engine::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
-				m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[1][i];
+				m_Eng->Engine::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
-				m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[2][i];
+				m_Eng->Engine::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]);
 			}
 			break;
 		}
@@ -79,18 +79,18 @@ void Engine_Ext_Dispersive::Apply2Voltages()
 			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
 			for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i)
 			{
-				eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[0][i];
+				eng_sse->Engine_sse::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
-				eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[1][i];
+				eng_sse->Engine_sse::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
-				eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[2][i];
+				eng_sse->Engine_sse::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]);
 			}
 			break;
 		}
 	default:
 		for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i)
 		{
-			m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[0][i];
+			m_Eng->SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
-			m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[1][i];
+			m_Eng->SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
-			m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[2][i];
+			m_Eng->SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]);
 		}
 		break;
 	}
@@ -109,9 +109,9 @@ void Engine_Ext_Dispersive::Apply2Current()
 		{
 			for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i)
 			{
-				m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[0][i];
+				m_Eng->Engine::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
-				m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[1][i];
+				m_Eng->Engine::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
-				m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[2][i];
+				m_Eng->Engine::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]);
 			}
 			break;
 		}
@@ -120,18 +120,18 @@ void Engine_Ext_Dispersive::Apply2Current()
 			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
 			for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i)
 			{
-				eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[0][i];
+				eng_sse->Engine_sse::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
-				eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[1][i];
+				eng_sse->Engine_sse::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
-				eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[2][i];
+				eng_sse->Engine_sse::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]);
 			}
 			break;
 		}
 	default:
 		for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i)
 		{
-			m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[0][i];
+			m_Eng->SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
-			m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[1][i];
+			m_Eng->SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
-			m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[2][i];
+			m_Eng->SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]);
 		}
 		break;
 	}

FDTD/engine_ext_mur_abc.cpp

@@ -194,8 +194,8 @@ void Engine_Ext_Mur_ABC::Apply2Voltages()
 			{
 				for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
 				{
-					m_Eng->Engine::GetVolt(m_nyP,pos) = m_volt_nyP[pos[m_nyP]][pos[m_nyPP]];
+					m_Eng->Engine::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]);
-					m_Eng->Engine::GetVolt(m_nyPP,pos) = m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]];
+					m_Eng->Engine::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]);
 				}
 			}
 			break;
@@ -208,8 +208,8 @@ void Engine_Ext_Mur_ABC::Apply2Voltages()
 			{
 				for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP])
 				{
-					eng_sse->Engine_sse::GetVolt(m_nyP,pos) = m_volt_nyP[pos[m_nyP]][pos[m_nyPP]];
+					eng_sse->Engine_sse::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]);
-					eng_sse->Engine_sse::GetVolt(m_nyPP,pos) = m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]];
+					eng_sse->Engine_sse::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]);
 				}
 			}
 			break;
@@ -220,8 +220,8 @@ void Engine_Ext_Mur_ABC::Apply2Voltages()
 		{
 			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->SetVolt(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]];
+				m_Eng->SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]);
 			}
 		}
 		break;

FDTD/engine_ext_pml_sf.cpp

@@ -181,9 +181,9 @@ void Engine_Ext_PML_SF_Plane::Apply2Voltages()
 			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->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->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->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->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->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]];
+				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]] );
 			}
 		}
 	}
@@ -254,9 +254,9 @@ void Engine_Ext_PML_SF_Plane::Apply2Current()
 			{
 				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->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->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->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->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]];
+				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]] );
 			}
 		}
 	}

FDTD/engine_multithread.cpp

@@ -44,8 +44,9 @@ Engine_Multithread* Engine_Multithread::New(const Operator_Multithread* op, unsi
 	return e;
 }
 
-Engine_Multithread::Engine_Multithread(const Operator_SSE_Compressed* op) : Engine_SSE_Compressed(op)
+Engine_Multithread::Engine_Multithread(const Operator_Multithread* op) : Engine_SSE_Compressed(op)
 {
+	m_Op_MT = op;
 	m_type = SSE;
 	m_barrier_VoltUpdate = 0;
 	m_barrier_VoltExcite = 0;
@@ -95,6 +96,10 @@ void Engine_Multithread::Init()
 	if (m_numThreads == 0)
 		m_numThreads = boost::thread::hardware_concurrency();
 
+	vector<unsigned int> m_Start_Lines;
+	vector<unsigned int> m_Stop_Lines;
+	m_Op_MT->CalcStartStopLines( m_numThreads, m_Start_Lines, m_Stop_Lines );
+
 	cout << "Multithreaded engine using " << m_numThreads << " threads. Utilization: (";
 	m_barrier_VoltUpdate = new boost::barrier(m_numThreads); // numThread workers
 	m_barrier_VoltExcite = new boost::barrier(m_numThreads+1); // numThread workers + 1 excitation thread
@@ -109,15 +114,14 @@ void Engine_Multithread::Init()
 	m_startBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
 	m_stopBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
 
-	unsigned int linesPerThread = round((float)numLines[0] / (float)m_numThreads);
 	for (unsigned int n=0; n<m_numThreads; n++)
 	{
-		unsigned int start = n * linesPerThread;
+		unsigned int start = m_Start_Lines.at(n);
-		unsigned int stop = (n+1) * linesPerThread - 1;
+		unsigned int stop = m_Stop_Lines.at(n);
 		unsigned int stop_h = stop;
-		if (n == m_numThreads-1) {
+		if (n == m_numThreads-1)
+		{
 			// last thread
-			stop = numLines[0]-1;
 			stop_h = stop-1;
 			cout << stop-start+1 << ")" << endl;
 		}
@@ -133,8 +137,9 @@ void Engine_Multithread::Init()
 
 void Engine_Multithread::Reset()
 {
-	if (!m_stopThreads) {
+	if (!m_stopThreads) // prevent multiple invocations
-		// prevent multiple invocations
+	{
+		ClearExtensions(); //prevent extensions from interfering with thread reset...
 
 		// stop the threads
 		//NS_Engine_Multithread::DBG().cout() << "stopping all threads" << endl;

FDTD/engine_multithread.h

@@ -78,6 +78,7 @@ class Engine_Multithread : public Engine_SSE_Compressed
 {
 	friend class NS_Engine_Multithread::thread;
 	friend class NS_Engine_Multithread::thread_e_excitation;
+	friend class Engine_CylinderMultiGrid;
 public:
 	static Engine_Multithread* New(const Operator_Multithread* op, unsigned int numThreads = 0);
 	virtual ~Engine_Multithread();
@@ -90,7 +91,8 @@ public:
 	virtual bool IterateTS(unsigned int iterTS);
 
 protected:
-	Engine_Multithread(const Operator_SSE_Compressed* op);
+	Engine_Multithread(const Operator_Multithread* op);
+	const Operator_Multithread* m_Op_MT;
 	boost::thread_group m_thread_group;
 	boost::barrier *m_startBarrier, *m_stopBarrier;
 	boost::barrier *m_barrier_VoltUpdate, *m_barrier_VoltExcite, *m_barrier_PreVolt, *m_barrier_PostVolt;

FDTD/engine_sse.h

@@ -33,10 +33,15 @@ public:
 	virtual unsigned int GetNumberOfTimesteps() {return numTS;};
 
 	//this access functions muss be overloaded by any new engine using a different storage model
-	inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return f4_volt[n][x][y][z%numVectors].f[z/numVectors]; }
+	inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z )	const { return f4_volt[n][x][y][z%numVectors].f[z/numVectors]; }
-	inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int pos[3] ) { return f4_volt[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; }
+	inline virtual FDTD_FLOAT GetVolt( unsigned int n, const unsigned int pos[3] )						const { return f4_volt[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; }
-	inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return f4_curr[n][x][y][z%numVectors].f[z/numVectors]; }
+	inline virtual FDTD_FLOAT GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z )	const { return f4_curr[n][x][y][z%numVectors].f[z/numVectors]; }
-	inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int pos[3] ) { return f4_curr[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; }
+	inline virtual FDTD_FLOAT GetCurr( unsigned int n, const unsigned int pos[3] )						const { return f4_curr[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; }
+
+	inline virtual void SetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value)	{ f4_volt[n][x][y][z%numVectors].f[z/numVectors]=value; }
+	inline virtual void SetVolt( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value )						{ f4_volt[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]=value; }
+	inline virtual void SetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value)	{ f4_curr[n][x][y][z%numVectors].f[z/numVectors]=value; }
+	inline virtual void SetCurr( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value )						{ f4_curr[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]=value; }
 
 protected:
 	Engine_sse(const Operator_sse* op);
@@ -47,6 +52,7 @@ protected:
 
 	unsigned int numVectors;
 
+public: //public access to the sse arrays for efficient extensions access... use careful...
 	f4vector**** f4_volt;
 	f4vector**** f4_curr;
 };

FDTD/operator.cpp

@@ -35,6 +35,8 @@ Operator::Operator()
 {
 	m_MeshType = ProcessFields::CARTESIAN_MESH;
 	Exc = 0;
+	dT = 0;
+	m_InvaildTimestep = false;
 }
 
 Operator::~Operator()
@@ -77,6 +79,7 @@ void Operator::Init()
 		m_BC[n]=0;
 
 	Exc = 0;
+	dT = 0;
 }
 
 void Operator::Reset()
@@ -302,7 +305,10 @@ void Operator::ShowStat() const
 	cout << "in " << GetDirName(1) << " direction\t\t: " << m_Nr_PEC[1] << endl;
 	cout << "in " << GetDirName(2) << " direction\t\t: " << m_Nr_PEC[2] << endl;
 	cout << "-----------------------------------" << endl;
-	cout << "Timestep (s)\t\t: " << dT << endl;
+	cout << "Timestep (s)\t\t: " << dT ;
+	if (opt_dT)
+		cout <<"\t(" << opt_dT << ")";
+	cout << endl;
 	cout << "Timestep method name\t: " << m_Used_TS_Name << endl;
 	cout << "Nyquist criteria (TS)\t: " << Exc->GetNyquistNum() << endl;
 	cout << "Nyquist criteria (s)\t: " << Exc->GetNyquistNum()*dT << endl;
@@ -458,7 +464,6 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo)
 {
 	if (geo==NULL) return false;
 
-	Reset();
 	CSX = geo;
 
 	CSRectGrid* grid=CSX->GetGrid();
@@ -527,6 +532,22 @@ int Operator::CalcECOperator()
 	if (Calc_EC()==0)
 		return -1;
 
+	m_InvaildTimestep = false;
+	opt_dT = 0;
+	if (dT>0)
+	{
+		double save_dT = dT;
+		CalcTimestep();
+		opt_dT = dT;
+		if (dT<save_dT)
+		{
+			cerr << "Operator::CalcECOperator: Warning, forced timestep: " << save_dT << "s is larger than calculated timestep: " << dT << "s! It is not recommended using this timestep!! " << endl;
+			m_InvaildTimestep = true;
+		}
+
+		dT = save_dT;
+	}
+	else
 		CalcTimestep();
 
 	InitOperator();
@@ -564,7 +585,8 @@ int Operator::CalcECOperator()
 	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;
+		if ((m_BC[n]==3) || (m_BC[n]==-1))
+			PEC[n] = false;
 	ApplyElectricBC(PEC);
 
 	InitExcitation();
@@ -600,8 +622,8 @@ void Operator::ApplyElectricBC(bool* dirs)
 				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]) = 0; // these are outside the FDTD-domain as defined by the main disc
+				GetVV(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; // 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
+				GetVI(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; // 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];

FDTD/operator.h

@@ -25,6 +25,7 @@
 
 class Operator_Extension;
 class Engine;
+class TiXmlElement;
 
 //! Abstract base-class for the FDTD-operator
 class Operator
@@ -50,6 +51,8 @@ public:
 	//! Calculate the effective/averaged material properties at the given position and direction ny.
 	virtual bool Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const;
 
+	virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) {return Exc->setupExcitation(Excite,maxTS);};
+
 	inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; }
 	inline virtual FDTD_FLOAT& GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; }
 
@@ -60,8 +63,11 @@ public:
 	virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
 	virtual void ApplyMagneticBC(bool* dirs);
 
+	//! Set a forced timestep to use by the operator
+	virtual void SetTimestep(double ts) {dT = ts;}
 	double GetTimestep() const {return dT;};
-	double GetNumberCells() const;
+	bool GetTimestepValid() const {return !m_InvaildTimestep;}
+	virtual double GetNumberCells() const;
 
 	//! Returns the number of lines as needed for post-processing etc. (for the engine, use GetOriginalNumLines())
 	virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];}
@@ -150,6 +156,8 @@ protected:
 	//Calc timestep only internal use
 	virtual double CalcTimestep();
 	double dT; //FDTD timestep!
+	double opt_dT;
+	bool m_InvaildTimestep;
 	string m_Used_TS_Name;
 
 	double CalcTimestep_Var1();

FDTD/operator_cylinder.cpp

@@ -30,31 +30,31 @@ Operator_Cylinder* Operator_Cylinder::New(unsigned int numThreads)
 	return op;
 }
 
-Operator_Cylinder::Operator_Cylinder() : __OP_CYLINDER_BASE_CLASS__()
+Operator_Cylinder::Operator_Cylinder() : Operator_Multithread()
 {
 	m_MeshType = ProcessFields::CYLINDRICAL_MESH;
 }
 
 Operator_Cylinder::~Operator_Cylinder()
 {
-	__OP_CYLINDER_BASE_CLASS__::Reset();
+	Operator_Multithread::Reset();
 }
 
 void Operator_Cylinder::Init()
 {
 	CC_closedAlpha = false;
 	CC_R0_included = false;
-	__OP_CYLINDER_BASE_CLASS__::Init();
+	Operator_Multithread::Init();
 }
 
 void Operator_Cylinder::Reset()
 {
-	__OP_CYLINDER_BASE_CLASS__::Reset();
+	Operator_Multithread::Reset();
 }
 
 void Operator_Cylinder::InitOperator()
 {
-	__OP_CYLINDER_BASE_CLASS__::InitOperator();
+	Operator_Multithread::InitOperator();
 }
 
 inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const
@@ -76,7 +76,7 @@ string Operator_Cylinder::GetDirName(int ny) const
 
 double Operator_Cylinder::GetMeshDelta(int n, const int* pos, bool dualMesh) const
 {
-	double delta = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(n,pos,dualMesh);
+	double delta = Operator_Multithread::GetMeshDelta(n,pos,dualMesh);
 	if (delta==0) return delta;
 	if (n==1)
 	{
@@ -103,7 +103,7 @@ double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) c
 {
 	if (ny==2)
 	{
-		double da = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(1,pos,dualMesh)/gridDelta;
+		double da = Operator_Multithread::GetMeshDelta(1,pos,dualMesh)/gridDelta;
 		double r1,r2;
 		if (!dualMesh)
 		{
@@ -127,12 +127,12 @@ double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) c
 		else
 			return fabs(MainOp->GetIndexDelta(1,pos[1]) * MainOp->GetIndexDelta(2,pos[2]) * GetDiscLine(0,pos[0],false) * gridDelta * gridDelta);
 	}
-	return __OP_CYLINDER_BASE_CLASS__::GetNodeArea(ny,pos,dualMesh);
+	return Operator_Multithread::GetNodeArea(ny,pos,dualMesh);
 }
 
 double Operator_Cylinder::GetEdgeLength(int ny, const int pos[3], bool dualMesh) const
 {
-	double length = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(ny,pos,dualMesh);
+	double length = Operator_Multithread::GetMeshDelta(ny,pos,dualMesh);
 	if (ny!=1)
 		return length;
 	return length * GetDiscLine(0,pos[0],dualMesh);
@@ -148,7 +148,7 @@ double Operator_Cylinder::GetEdgeArea(int ny, const int pos[3], bool dualMesh) c
 
 bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
 {
-	if (__OP_CYLINDER_BASE_CLASS__::SetGeometryCSX(geo)==false) return false;
+	if (Operator_Multithread::SetGeometryCSX(geo)==false) return false;
 
 	double minmaxA = fabs(discLines[1][numLines[1]-1]-discLines[1][0]);
 	if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width
@@ -209,7 +209,7 @@ void Operator_Cylinder::ApplyElectricBC(bool* dirs)
 			}
 		}
 	}
-	__OP_CYLINDER_BASE_CLASS__::ApplyElectricBC(dirs);
+	Operator_Multithread::ApplyElectricBC(dirs);
 }
 
 void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
@@ -223,7 +223,7 @@ void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
 	{
 		dirs[0]=0;  //no PMC in r_min directions...
 	}
-	__OP_CYLINDER_BASE_CLASS__::ApplyMagneticBC(dirs);
+	Operator_Multithread::ApplyMagneticBC(dirs);
 }
 
 void Operator_Cylinder::AddExtension(Operator_Extension* op_ext)

FDTD/operator_cylinder.h

@@ -18,9 +18,6 @@
 #ifndef OPERATOR_CYLINDER_H
 #define OPERATOR_CYLINDER_H
 
-//! define the base class for the cylindrical coordinate FDTD
-#define __OP_CYLINDER_BASE_CLASS__ Operator_Multithread
-
 #include "operator_multithread.h"
 
 //! This class creates an operator for a cylindrical FDTD.
@@ -28,21 +25,15 @@
 This class creates an operator for a cylindrical FDTD. No special engine is necessary,
 all special cases e.g. a closed alpha mesh or an included r=0 case is treated by an operator/engine extension \sa operator_ext_cylinder.
 */
-class Operator_Cylinder : public __OP_CYLINDER_BASE_CLASS__
+class Operator_Cylinder : public Operator_Multithread
 {
+	friend class Operator_CylinderMultiGrid;
 public:
 	static Operator_Cylinder* New(unsigned int numThreads = 0);
 	virtual ~Operator_Cylinder();
 
 	virtual bool SetGeometryCSX(ContinuousStructure* geo);
 
-
-//	virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC) const;
-//
-//
-//	//! Calculate the effective/averaged material properties at the given position and direction ny.
-//	virtual bool Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const;
-
 	virtual void ApplyElectricBC(bool* dirs);
 	virtual void ApplyMagneticBC(bool* dirs);
 

FDTD/operator_cylindermultigrid.cpp

@@ -0,0 +1,217 @@
+/*
+*	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_cylindermultigrid.h"
+#include "engine_cylindermultigrid.h"
+#include "operator_ext_cylinder.h"
+
+Operator_CylinderMultiGrid::Operator_CylinderMultiGrid(vector<double> Split_Radii) : Operator_Cylinder()
+{
+	m_Split_Radii = Split_Radii;
+	m_Split_Rad = m_Split_Radii.back();
+	m_Split_Radii.pop_back();
+}
+
+Operator_CylinderMultiGrid::~Operator_CylinderMultiGrid()
+{
+}
+
+Operator_CylinderMultiGrid* Operator_CylinderMultiGrid::New(vector<double> Split_Radii, unsigned int numThreads)
+{
+	cout << "Create cylindrical multi grid FDTD operator " << endl;
+	Operator_CylinderMultiGrid* op = new Operator_CylinderMultiGrid(Split_Radii);
+	op->setNumThreads(numThreads);
+	op->Init();
+
+	return op;
+}
+
+Engine* Operator_CylinderMultiGrid::CreateEngine() const
+{
+	Engine_CylinderMultiGrid* eng = Engine_CylinderMultiGrid::New(this,m_numThreads);
+	return eng;
+}
+
+double Operator_CylinderMultiGrid::GetNumberCells() const
+{
+	if (numLines)
+		return (numLines[0]-m_Split_Pos)*(numLines[1])*(numLines[2]) + m_InnerOp->GetNumberCells();
+	return 0;
+}
+
+bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
+{
+	if (Operator_Cylinder::SetGeometryCSX(geo)==false)
+		return false;
+
+	if (numLines[1]%2 != 1)
+	{
+		cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, number of line in alpha direction must be odd... found: " << numLines[1] << endl;
+		exit(0);
+	}
+
+	//check if mesh is homogenous in alpha-direction
+	double diff=discLines[1][1]-discLines[1][0];
+	for (unsigned int n=2;n<numLines[1];++n)
+	{
+		if ( fabs((discLines[1][n]-discLines[1][n-1]) - diff)/diff > 1e-10)
+		{
+			cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, mesh has to be homogenous in alpha direction for multi grid engine, violation found at: " << n << endl;
+			exit(0);
+		}
+	}
+
+	m_Split_Pos = 0;
+	for (unsigned int n=0;n<numLines[0];++n)
+	{
+		if (m_Split_Rad < discLines[0][n])
+		{
+			m_Split_Pos = n;
+			cout << "Operator_CylinderMultiGrid::SetGeometryCSX: Found mesh split position @" << m_Split_Pos << endl;
+			m_Split_Rad = discLines[0][n];
+			break;
+		}
+	}
+	if ((m_Split_Pos<4) || (m_Split_Pos>numLines[0]-4))
+	{
+		cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, split invalid..." << endl;
+		return false;
+	}
+
+	CSRectGrid* grid = geo->GetGrid();
+
+	grid->ClearLines(0);
+	grid->ClearLines(1);
+	for (unsigned int n=0; n<m_Split_Pos ;++n)
+		grid->AddDiscLine(0,discLines[0][n]);
+	for (unsigned int n=0; n<numLines[1]; n+=2)
+		grid->AddDiscLine(1,discLines[1][n]);
+
+	if (m_InnerOp->SetGeometryCSX(CSX)==false)
+		return false;
+
+	//restore grid to original mesh
+	grid->ClearLines(0);
+	grid->ClearLines(1);
+	for (unsigned int n=0; n<numLines[0]; ++n)
+		grid->AddDiscLine(0,discLines[0][n]);
+	for (unsigned int n=0; n<numLines[1]; ++n)
+		grid->AddDiscLine(1,discLines[1][n]);
+
+	return true;
+}
+
+void Operator_CylinderMultiGrid::Init()
+{
+	Operator_Cylinder::Init();
+
+	if (m_Split_Radii.empty())
+		m_InnerOp = Operator_Cylinder::New(m_numThreads);
+	else
+		m_InnerOp = Operator_CylinderMultiGrid::New(m_Split_Radii,m_numThreads);
+}
+
+void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const
+{
+	if (numLines[0]<numThreads) //in case more threads requested as lines in r-direction, reduce number of worker threads
+		numThreads = numLines[0];
+
+	unsigned int linesPerThread = round((float)(numLines[0] - m_Split_Pos + 1) / (float)numThreads);
+	if ((numThreads-1) * linesPerThread >= (numLines[0] - m_Split_Pos + 1))
+		--numThreads;
+
+	start.resize(numThreads);
+	stop.resize(numThreads);
+
+	for (unsigned int n=0; n<numThreads; n++)
+	{
+		start.at(n) = n * linesPerThread + m_Split_Pos - 1;
+		stop.at(n) = (n+1) * linesPerThread - 1 + m_Split_Pos - 1;
+		if (n == numThreads-1) // last thread
+			stop.at(n) = numLines[0]-1;
+	}
+}
+
+
+int Operator_CylinderMultiGrid::CalcECOperator()
+{
+	int retCode=0;
+	if (dT)
+		m_InnerOp->SetTimestep(dT);
+
+	//calc inner child first
+	m_InnerOp->CalcECOperator();
+
+	dT = m_InnerOp->GetTimestep();
+
+	retCode = Operator_Cylinder::CalcECOperator();
+	if (GetTimestepValid()==false)
+	{
+		cerr << "Operator_CylinderMultiGrid::CalcECOperator(): Warning, timestep invalid... resetting..." << endl;
+		dT = opt_dT;
+		m_InnerOp->SetTimestep(dT);
+		m_InnerOp->CalcECOperator();
+		return Operator_Cylinder::CalcECOperator();
+	}
+
+	return retCode;
+}
+
+bool Operator_CylinderMultiGrid::SetupExcitation(TiXmlElement* Excite, unsigned int maxTS)
+{
+	if (!m_InnerOp->SetupExcitation(Excite,maxTS))
+		return false;
+	return Exc->setupExcitation(Excite,maxTS);
+}
+
+void Operator_CylinderMultiGrid::Reset()
+{
+	Operator_Cylinder::Reset();
+	m_InnerOp->Reset();
+}
+
+void Operator_CylinderMultiGrid::SetBoundaryCondition(int* BCs)
+{
+	Operator_Cylinder::SetBoundaryCondition(BCs);
+	BCs[1] = 0; //always PEC in +r-direction
+	m_InnerOp->SetBoundaryCondition(BCs);
+}
+
+void Operator_CylinderMultiGrid::AddExtension(Operator_Extension* op_ext)
+{
+	if (dynamic_cast<Operator_Ext_Cylinder*>(op_ext))
+	{
+		Operator_Cylinder::AddExtension(op_ext);
+		return;
+	}
+	Operator_Extension* child_Ext = op_ext->Clone(m_InnerOp);
+	if (child_Ext==NULL)
+	{
+		cerr << "Operator_CylinderMultiGrid::AddExtension: Warning, extension: " << op_ext->GetExtensionName() << " can not be cloned for the child operator. Skipping Extension... " << endl;
+		return;
+	}
+	Operator_Cylinder::AddExtension(op_ext);
+	//give the copy to child
+	m_InnerOp->AddExtension(child_Ext);
+}
+
+void Operator_CylinderMultiGrid::ShowStat() const
+{
+	m_InnerOp->ShowStat();
+	m_InnerOp->ShowExtStat();
+	Operator_Cylinder::ShowStat();
+}

FDTD/operator_cylindermultigrid.h

@@ -0,0 +1,71 @@
+/*
+*	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_CYLINDERMULTIGRID_H
+#define OPERATOR_CYLINDERMULTIGRID_H
+
+#include "operator_cylinder.h"
+
+//! This is a cylindrical FDTD operator using a simple multi-grid approach.
+/*!
+  This cylindrical multi-grid operator itself is not calculating any real operator, instead it is hosting two separate "child" operators of type "Operator_Cylinder".
+  This operator class (or the corresponding engine) will perform the interpolation and connection between these two child-operator/engines.
+  One of the child operators itself may be another multi-grid operator to allow for a cascaded multi-grid approach.
+  */
+class Operator_CylinderMultiGrid : public Operator_Cylinder
+{
+	friend class Engine_CylinderMultiGrid;
+public:
+	static Operator_CylinderMultiGrid* New(vector<double> Split_Radii, unsigned int numThreads = 0);
+	virtual ~Operator_CylinderMultiGrid();
+
+	virtual double GetNumberCells() const;
+
+	virtual Engine* CreateEngine() const;
+
+	virtual bool SetGeometryCSX(ContinuousStructure* geo);
+
+	virtual unsigned int GetSplitPos() const {return m_Split_Pos;}
+
+	virtual int CalcECOperator();
+
+	virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS);
+
+	virtual void SetBoundaryCondition(int* BCs);
+
+	virtual void AddExtension(Operator_Extension* op_ext);
+
+	Operator_Cylinder* GetInnerOperator() const {return m_InnerOp;}
+
+	virtual void ShowStat() const;
+
+protected:
+	Operator_CylinderMultiGrid(vector<double> Split_Radii);
+	virtual void Init();
+//	virtual void InitOperator();
+	virtual void Reset();
+
+	double m_Split_Rad;
+	vector<double> m_Split_Radii;
+	unsigned int m_Split_Pos;
+
+	Operator_Cylinder* m_InnerOp;
+
+	virtual void CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const;
+};
+
+#endif // OPERATOR_CYLINDERMULTIGRID_H

FDTD/operator_ext_cylinder.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_ext_cylinder.h"
 #include "operator_cylinder.h"
 #include "engine_ext_cylinder.h"

FDTD/operator_ext_cylinder.h

@@ -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/>.
+*/
+
 #ifndef OPERATOR_EXT_CYLINDER_H
 #define OPERATOR_EXT_CYLINDER_H
 

FDTD/operator_ext_mur_abc.cpp

@@ -21,6 +21,33 @@
 #include "tools/array_ops.h"
 
 Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op)
+{
+	Initialize();
+}
+
+Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC()
+{
+	Delete2DArray(m_Mur_Coeff_nyP,m_numLines);
+	m_Mur_Coeff_nyP = NULL;
+	Delete2DArray(m_Mur_Coeff_nyPP,m_numLines);
+	m_Mur_Coeff_nyPP = NULL;
+}
+
+Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op, Operator_Ext_Mur_ABC* op_ext) : Operator_Extension(op, op_ext)
+{
+	Initialize();
+	m_v_phase = op_ext->m_v_phase;
+	SetDirection(op_ext->m_ny,op_ext->m_top);
+}
+
+Operator_Extension* Operator_Ext_Mur_ABC::Clone(Operator* op)
+{
+	if (dynamic_cast<Operator_Ext_Mur_ABC*>(this)==NULL)
+		return NULL;
+	return new Operator_Ext_Mur_ABC(op, this);
+}
+
+void Operator_Ext_Mur_ABC::Initialize()
 {
 	m_ny = -1;
 	m_nyP = -1;
@@ -37,14 +64,6 @@ Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op
 	m_numLines[1]=0;
 }
 
-Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC()
-{
-	Delete2DArray(m_Mur_Coeff_nyP,m_numLines);
-	m_Mur_Coeff_nyP = NULL;
-	Delete2DArray(m_Mur_Coeff_nyPP,m_numLines);
-	m_Mur_Coeff_nyPP = NULL;
-}
-
 void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny)
 {
 	if ((ny<0) || (ny>2))
@@ -54,6 +73,7 @@ void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny)
 	Delete2DArray(m_Mur_Coeff_nyPP,m_numLines);
 
 	m_ny = ny;
+	m_top = top_ny;
 	m_nyP = (ny+1)%3;
 	m_nyPP = (ny+2)%3;
 	if (!top_ny)

FDTD/operator_ext_mur_abc.h

@@ -28,6 +28,9 @@ public:
 	Operator_Ext_Mur_ABC(Operator* op);
 	~Operator_Ext_Mur_ABC();
 
+	Operator_Ext_Mur_ABC(Operator* op, Operator_Ext_Mur_ABC* op_ext);
+	virtual Operator_Extension* Clone(Operator* op);
+
 	//! Define the direction of this ABC: \a 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);
 
@@ -45,8 +48,10 @@ public:
 	virtual void ShowStat(ostream &ostr) const;
 
 protected:
+	void Initialize();
 	int m_ny;
 	int m_nyP,m_nyPP;
+	bool m_top;
 	unsigned int m_LineNr;
 	int m_LineNr_Shift;
 

FDTD/operator_extension.cpp

@@ -24,7 +24,20 @@ Operator_Extension::Operator_Extension(Operator* op)
 	m_Op = op;
 }
 
+Operator_Extension::Operator_Extension(Operator* op, Operator_Extension* op_ext)
+{
+	UNUSED(op_ext);
+	m_Op = op;
+}
+
 void Operator_Extension::ShowStat(ostream &ostr) const
 {
 	ostr << "--- " << GetExtensionName() << " ---" << endl;
 }
+
+Operator_Extension* Operator_Extension::Clone(Operator* op)
+{
+	if (dynamic_cast<Operator_Extension*>(this)==NULL)
+		return NULL;
+	return new Operator_Extension(op, this);
+}

FDTD/operator_extension.h

@@ -35,6 +35,13 @@ class Operator_Extension
 {
 	friend class Engine_Extension;
 public:
+	//! Create a clone of this extension, will return NULL if this is impossible
+	/*!
+		Create a clone of this extension, will return NULL if this is impossible (e.g. derived extension has no clone method and copy-constructor)...
+		BuildExtension has to be called separatly!
+	*/
+	virtual Operator_Extension* Clone(Operator* op);
+
 	virtual bool BuildExtension() {return true;}
 
 	virtual Engine_Extension* CreateEngineExtention() {return 0;}
@@ -47,6 +54,8 @@ public:
 
 protected:
 	Operator_Extension(Operator* op);
+	//! Copy constructor, returns NULL if extension cannot be copied...
+	Operator_Extension(Operator* op, Operator_Extension* op_ext);
 	Operator* m_Op;
 };
 

FDTD/operator_multithread.cpp

@@ -74,14 +74,35 @@ void Operator_Multithread::Reset()
 	delete m_CalcPEC_Stop;m_CalcPEC_Stop=NULL;
 }
 
+void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const
+{
+	if (numLines[0]<numThreads) //in case more threads requested as lines in x-direction, reduce number of worker threads
+		numThreads = numLines[0];
+
+	unsigned int linesPerThread = round((float)numLines[0] / (float)numThreads);
+	if ((numThreads-1) * linesPerThread >= numLines[0])
+		--numThreads;
+
+	start.resize(numThreads);
+	stop.resize(numThreads);
+
+	for (unsigned int n=0; n<numThreads; n++)
+	{
+		start.at(n) = n * linesPerThread;
+		stop.at(n) = (n+1) * linesPerThread - 1;
+		if (n == numThreads-1) // last thread
+			stop.at(n) = numLines[0]-1;
+	}
+}
+
 int Operator_Multithread::CalcECOperator()
 {
 	if (m_numThreads == 0)
 		m_numThreads = boost::thread::hardware_concurrency();
 
-	unsigned int linesPerThread = round((float)numLines[0] / (float)m_numThreads);
+	vector<unsigned int> m_Start_Lines;
-	if ((m_numThreads-1) * linesPerThread >= numLines[0])
+	vector<unsigned int> m_Stop_Lines;
-		--m_numThreads;
+	CalcStartStopLines( m_numThreads, m_Start_Lines, m_Stop_Lines );
 
 	cout << "Multithreaded operator using " << m_numThreads << " threads." << std::endl;
 
@@ -94,12 +115,7 @@ int Operator_Multithread::CalcECOperator()
 
 	for (unsigned int n=0; n<m_numThreads; n++)
 	{
-		unsigned int start = n * linesPerThread;
+		boost::thread *t = new boost::thread( Operator_Thread(this,m_Start_Lines.at(n),m_Stop_Lines.at(n),n) );
-		unsigned int stop = (n+1) * linesPerThread - 1;
-		if (n == m_numThreads-1) // last thread
-			stop = numLines[0]-1;
-
-		boost::thread *t = new boost::thread( Operator_Thread(this,start,stop,n) );
 		m_thread_group.add_thread( t );
 	}
 

FDTD/operator_multithread.h

@@ -24,6 +24,7 @@
 
 class Operator_Multithread : public Operator_SSE_Compressed
 {
+	friend class Engine_Multithread;
 	friend class Operator_Thread;
 public:
 	//! Create a new operator
@@ -55,6 +56,14 @@ protected:
 
 	boost::thread_group m_thread_group;
 	unsigned int m_numThreads; // number of worker threads
+
+	//! Calculate the start/stop lines for the multithreading operator and engine.
+	/*!
+		It depends on the number of threads and number of lines to simulate.
+		This method is also used by the multithreading engine!
+		This method may also reduce the usable number of thread in case of too few lines or otherwise bad utilization.
+	*/
+	virtual void CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const;
 };
 
 class Operator_Thread

matlab/examples/Coax_Cylindrical_MG.m

@@ -0,0 +1,155 @@
+close all
+clear
+clc
+
+%example for an cylindrical mesh, modeling a coaxial cable
+% this example is using a multi-grid approach
+
+
+%% setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Settings = [];
+Settings.LogFile = 'openEMS.log';
+
+physical_constants
+
+f0 = 0.5e9;
+epsR = 1;       %material filling
+
+length = 1000;
+port_dist = length/2;
+rad_i  = 10;    %inner radius
+rad_a = 200;    %outer radius
+partial = 0.5;  %e.g. 0.5 means only one half of a coax, should be <1 or change boundary cond.
+max_mesh = 10 / sqrt(epsR);
+max_alpha = max_mesh;
+N_alpha = ceil(rad_a * 2*pi * partial / max_alpha);
+
+%make it even...
+N_alpha = N_alpha + mod(N_alpha,2);
+%make sure it is multiple of 4, needed for 2 multi-grid steps
+N_alpha = ceil((N_alpha)/4) *4 + 1;
+
+openEMS_opts = '';
+% openEMS_opts = [openEMS_opts ' --disable-dumps'];
+% openEMS_opts = [openEMS_opts ' --debug-material'];
+% openEMS_opts = [openEMS_opts ' --numThreads=1'];
+
+def_refSimu = 0; % do a reference simulation without the multi-grid
+
+%% setup done %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+if (def_refSimu>0)
+    Sim_Path = 'tmp_ref';
+else
+    Sim_Path = 'tmp';
+end
+Sim_CSX = 'coax.xml';
+
+if (exist(Sim_Path,'dir'))
+    rmdir(Sim_Path,'s');
+end
+mkdir(Sim_Path);
+
+%setup FDTD parameter
+if (def_refSimu>0)
+    FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',5 );
+else
+    FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',5 ,'MultiGrid','60,120');
+end
+FDTD = SetGaussExcite(FDTD,f0,f0);
+BC = [0 0 1 1 2 2];
+FDTD = SetBoundaryCond(FDTD,BC);
+
+mesh_res = [max_mesh 2*pi*partial/N_alpha max_mesh];
+
+%setup CSXCAD geometry
+CSX = InitCSX();
+mesh.x = SmoothMeshLines([rad_i rad_a],mesh_res(1));
+mesh.y = linspace(-pi*partial,pi*partial,N_alpha);
+mesh.z = SmoothMeshLines([0 port_dist length],mesh_res(3));
+CSX = DefineRectGrid(CSX, 1e-3,mesh);
+
+start = [rad_i mesh.y(1)   mesh.z(3)];
+stop  = [rad_a mesh.y(end) mesh.z(3)];
+
+CSX = AddExcitation(CSX,'excite',0,[1 0 0]);
+weight{1} = '1/rho';
+weight{2} = 0;
+weight{3} = 0;
+CSX = SetExcitationWeight(CSX, 'excite', weight );
+CSX = AddBox(CSX,'excite',0 ,start,stop);
+ 
+
+start = [mesh.x(1)   mesh.y(1)   mesh.z(1)];
+stop  = [mesh.x(end) mesh.y(end) mesh.z(end)];
+CSX = AddMaterial(CSX,'material');
+CSX = SetMaterialProperty(CSX,'material','Epsilon',epsR);
+CSX = AddBox(CSX,'material',0 ,start,stop);
+
+%dump
+CSX = AddDump(CSX,'Et_rz_','DumpMode',0);
+start = [mesh.x(1)   0 mesh.z(1)];
+stop  = [mesh.x(end) 0 mesh.z(end)];
+CSX = AddBox(CSX,'Et_rz_',0 , start,stop);
+
+CSX = AddDump(CSX,'Ht_rz_','DumpType',1,'DumpMode',0);
+CSX = AddBox(CSX,'Ht_rz_',0 , start,stop);
+
+CSX = AddDump(CSX,'Et_','DumpType',0,'DumpMode',0);
+start = [mesh.x(1)   mesh.y(1)   length/2];
+stop  = [mesh.x(end) mesh.y(end) length/2];
+CSX = AddBox(CSX,'Et_',0,start,stop);
+
+CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0);
+start = [mesh.x(1)   mesh.y(1)   length/2];
+stop  = [mesh.x(end) mesh.y(end) length/2];
+CSX = AddBox(CSX,'Ht_',0,start,stop);
+
+% voltage calc (take a voltage average to be at the same spot as the
+% current calculation)
+CSX = AddProbe(CSX,'ut1_1',0);
+start = [ rad_i 0 port_dist ];stop = [ rad_a 0 port_dist ];
+CSX = AddBox(CSX,'ut1_1', 0 ,start,stop);
+CSX = AddProbe(CSX,'ut1_2',0);
+start = [ rad_i 0 port_dist+mesh_res(3) ];stop = [ rad_a 0 port_dist+mesh_res(3) ];
+CSX = AddBox(CSX,'ut1_2', 0 ,start,stop);
+
+% current calc
+CSX = AddProbe(CSX,'it1',1);
+mid = 75;
+start = [ 0 mesh.y(1) port_dist+mesh_res(3)/2 ];stop = [ mid mesh.y(end) port_dist+mesh_res(3)/2 ];
+CSX = AddBox(CSX,'it1', 0 ,start,stop);
+
+%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
+
+RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts, Settings)
+
+%%
+close all
+freq = linspace(0,2*f0,201);
+UI = ReadUI({'ut1_1','ut1_2','it1'},Sim_Path,freq);
+u_f = (UI.FD{1}.val + UI.FD{2}.val)/2;          %averaging voltages to fit current
+i_f = UI.FD{3}.val / partial;
+
+% plot(UI.TD{1}.t,UI.TD{1}.val);
+% grid on;
+% 
+% figure
+% plot(UI.TD{3}.t,UI.TD{3}.val);
+% grid on;
+
+%plot Z_L compare
+figure
+ZL = Z0/2/pi/sqrt(epsR)*log(rad_a/rad_i); %analytic line-impedance of a coax
+plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g','Linewidth',3);
+hold on;
+grid on;
+Z = u_f./i_f;
+plot(UI.FD{1}.f,real(Z),'k--','Linewidth',2);
+plot(UI.FD{1}.f,imag(Z),'r-','Linewidth',2);
+xlim([0 2*f0]);
+legend('Z_L - analytic','\Re\{Z\} - FDTD','\Im\{Z\} - FDTD','Location','Best');
+
+

openEMS.pro

@@ -81,7 +81,10 @@ SOURCES += main.cpp \
     FDTD/engine_ext_lorentzmaterial.cpp \
     FDTD/operator_ext_pml_sf.cpp \
     FDTD/engine_ext_pml_sf.cpp \
-    FDTD/processmodematch.cpp
+    FDTD/processmodematch.cpp \
+    FDTD/operator_cylindermultigrid.cpp \
+    FDTD/engine_cylindermultigrid.cpp \
+    FDTD/engine_ext_cylindermultigrid.cpp
 HEADERS += tools/ErrorMsg.h \
     tools/AdrOp.h \
     tools/constants.h \
@@ -119,7 +122,10 @@ HEADERS += tools/ErrorMsg.h \
     FDTD/engine_ext_lorentzmaterial.h \
     FDTD/operator_ext_pml_sf.h \
     FDTD/engine_ext_pml_sf.h \
-    FDTD/processmodematch.h
+    FDTD/processmodematch.h \
+    FDTD/operator_cylindermultigrid.h \
+    FDTD/engine_cylindermultigrid.h \
+    FDTD/engine_ext_cylindermultigrid.h
 QMAKE_CXXFLAGS_RELEASE = -O3 \
     -g \
 	-march=native

openems.cpp

@@ -19,6 +19,7 @@
 #include <iomanip>
 #include "tools/array_ops.h"
 #include "FDTD/operator_cylinder.h"
+#include "FDTD/operator_cylindermultigrid.h"
 #include "FDTD/engine_multithread.h"
 #include "FDTD/operator_multithread.h"
 #include "FDTD/operator_ext_mur_abc.h"
@@ -293,7 +294,6 @@ int openEMS::SetupFDTD(const char* file)
 	double maxTime=0;
 	FDTD_Opts->QueryDoubleAttribute("MaxTime",&maxTime);
 
-
 	TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond");
 	if (BC==NULL)
 	{
@@ -319,7 +319,15 @@ int openEMS::SetupFDTD(const char* file)
 	//*************** setup operator ************//
 	if (CylinderCoords)
 	{
+		const char* radii = FDTD_Opts->Attribute("MultiGrid");
+		if (radii)
+		{
+			string rad(radii);
+			FDTD_Op = Operator_CylinderMultiGrid::New(SplitString2Double(rad,','),m_engine_numThreads);
+		}
+		else
 			FDTD_Op = Operator_Cylinder::New(m_engine_numThreads);
+		CSX.SetCoordInputType(1); //tell CSX to use cylinder-coords
 	}
 	else if (m_engine == EngineType_SSE)
 	{
@@ -345,13 +353,18 @@ int openEMS::SetupFDTD(const char* file)
 	if (CSX.GetQtyPropertyType(CSProperties::LORENTZMATERIAL)>0)
 		FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op));
 
+	double timestep=0;
+	FDTD_Opts->QueryDoubleAttribute("TimeStep",&timestep);
+	if (timestep)
+		FDTD_Op->SetTimestep(timestep);
+
 	FDTD_Op->CalcECOperator();
 	
 	unsigned int maxTime_TS = (unsigned int)(maxTime/FDTD_Op->GetTimestep());
 	if ((maxTime_TS>0) && (maxTime_TS<NrTS))
 		NrTS = maxTime_TS;
 
-	if (!FDTD_Op->Exc->setupExcitation( FDTD_Opts->FirstChildElement("Excitation"), NrTS ))
+	if (!FDTD_Op->SetupExcitation( FDTD_Opts->FirstChildElement("Excitation"), NrTS ))
 		exit(2);
 
 	if (DebugMat)