Moved special cylinder operations into op extentions

FDTD/engine.cpp

@@ -34,7 +34,8 @@ Engine::Engine(const Operator* op)
 	Op = op;
 	for (int n=0;n<3;++n)
 	{
-		numLines[n] = Op->GetNumberOfLines(n);
+//		numLines[n] = Op->GetNumberOfLines(n);
+		numLines[n] = Op->numLines[n];
 	}
 
 	for (size_t n=0;n<Op->GetNumberOfExtentions();++n)

FDTD/engine_cylinder.cpp

@@ -1,136 +0,0 @@
-/*
-*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
-*
-*	This program is free software: you can redistribute it and/or modify
-*	it under the terms of the GNU General Public License as published by
-*	the Free Software Foundation, either version 3 of the License, or
-*	(at your option) any later version.
-*
-*	This program is distributed in the hope that it will be useful,
-*	but WITHOUT ANY WARRANTY; without even the implied warranty of
-*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-*	GNU General Public License for more details.
-*
-*	You should have received a copy of the GNU General Public License
-*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include "engine_cylinder.h"
-#include "engine_extension.h"
-#include "operator_extension.h"
-
-Engine_Cylinder* Engine_Cylinder::New(const Operator_Cylinder* op)
-{
-	Engine_Cylinder* e = new Engine_Cylinder(op);
-	e->Init();
-	return e;
-}
-
-Engine_Cylinder::Engine_Cylinder(const Operator_Cylinder* op) : Engine(op)
-{
-	cyl_Op = op;
-	if (cyl_Op->GetClosedAlpha())
-	{
-		++numLines[1]; //necessary for dobled voltage and current line in alpha-dir, operator will return one smaller for correct post-processing
-	}
-}
-
-inline void Engine_Cylinder::CloseAlphaVoltages()
-{
-	unsigned int pos[3];
-	// copy voltages from last alpha-plane to first
-	unsigned int last_A_Line = numLines[1]-1;
-	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
-	{
-		for (pos[2]=0;pos[2]<numLines[2];++pos[2])
-		{
-			volt[0][pos[0]][0][pos[2]] = volt[0][pos[0]][last_A_Line][pos[2]];
-			volt[1][pos[0]][0][pos[2]] = volt[1][pos[0]][last_A_Line][pos[2]];
-			volt[2][pos[0]][0][pos[2]] = volt[2][pos[0]][last_A_Line][pos[2]];
-		}
-	}
-}
-
-void Engine_Cylinder::R0IncludeVoltages()
-{
-	unsigned int pos[3];
-	pos[0] = 0;
-	for (pos[2]=0;pos[2]<numLines[2];++pos[2])
-	{
-		volt[2][0][0][pos[2]] *= cyl_Op->vv_R0[pos[2]];
-		for (pos[1]=0;pos[1]<numLines[1]-cyl_Op->GetClosedAlpha();++pos[1])
-		{
-				volt[2][0][0][pos[2]] += cyl_Op->vi_R0[pos[2]] *  curr[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])
-		{
-			volt[1][0][pos[1]][pos[2]] = 0; //no voltage in alpha-direction at r=0
-			volt[2][0][pos[1]][pos[2]] = volt[2][0][0][pos[2]];
-		}
-	}
-}
-
-inline void Engine_Cylinder::CloseAlphaCurrents()
-{
-	unsigned int pos[3];
-	// copy currents from first alpha-plane to last
-	for (pos[0]=0;pos[0]<numLines[0]-1;++pos[0])
-	{
-		unsigned int last_A_Line = numLines[1]-1;
-		for (pos[2]=0;pos[2]<numLines[2]-1;++pos[2])
-		{
-			curr[0][pos[0]][last_A_Line][pos[2]] = curr[0][pos[0]][0][pos[2]];
-			curr[1][pos[0]][last_A_Line][pos[2]] = curr[1][pos[0]][0][pos[2]];
-			curr[2][pos[0]][last_A_Line][pos[2]] = curr[2][pos[0]][0][pos[2]];
-		}
-	}
-}
-
-bool Engine_Cylinder::IterateTS(unsigned int iterTS)
-{
-	if (cyl_Op->GetClosedAlpha()==false)
-		return Engine::IterateTS(iterTS);
-
-	for (unsigned int iter=0;iter<iterTS;++iter)
-	{
-		//voltage updates with extensions
-		for (size_t n=0;n<m_Eng_exts.size();++n)
-			m_Eng_exts.at(n)->DoPreVoltageUpdates();
-
-		UpdateVoltages();
-
-		for (size_t n=0;n<m_Eng_exts.size();++n)
-			m_Eng_exts.at(n)->DoPostVoltageUpdates();
-		for (size_t n=0;n<m_Eng_exts.size();++n)
-			m_Eng_exts.at(n)->Apply2Voltages();
-
-		if (cyl_Op->GetR0Included())
-			R0IncludeVoltages();
-
-		ApplyVoltageExcite();
-
-		CloseAlphaVoltages();
-
-		//current updates with extensions
-		for (size_t n=0;n<m_Eng_exts.size();++n)
-			m_Eng_exts.at(n)->DoPreCurrentUpdates();
-
-		UpdateCurrents();
-
-		for (size_t n=0;n<m_Eng_exts.size();++n)
-			m_Eng_exts.at(n)->DoPostCurrentUpdates();
-		for (size_t n=0;n<m_Eng_exts.size();++n)
-			m_Eng_exts.at(n)->Apply2Current();
-
-		ApplyCurrentExcite();
-
-		CloseAlphaCurrents();
-
-		++numTS;
-	}
-
-	return true;
-}

FDTD/engine_ext_cylinder.cpp

@@ -0,0 +1,89 @@
+/*
+*	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_cylinder.h"
+#include "operator_ext_cylinder.h"
+
+Engine_Ext_Cylinder::Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext) : Engine_Extension(op_ext)
+{
+	CC_closedAlpha = op_ext->CC_closedAlpha;
+	CC_R0_included = op_ext->CC_R0_included;
+
+	for (int n=0;n<3;++n)
+		numLines[n] = op_ext->m_Op->GetNumberOfLines(n);
+}
+
+void Engine_Ext_Cylinder::Apply2Voltages()
+{
+	if (CC_closedAlpha==false) return;
+
+	if (CC_R0_included)
+	{
+		unsigned int pos[3];
+		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]];
+			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]);
+			}
+		}
+		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->GetVolt(2,0,pos[1],pos[2]) = m_Eng->GetVolt(2,0,0,pos[2]);
+			}
+		}
+	}
+
+	//close alpha
+	unsigned int pos[3];
+	// copy voltages from last alpha-plane to first
+	unsigned int last_A_Line = numLines[1]-1;
+	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+	{
+		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->GetVolt(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]);
+		}
+	}
+}
+
+void Engine_Ext_Cylinder::Apply2Current()
+{
+	if (CC_closedAlpha==false) return;
+
+	//close alpha
+	unsigned int pos[3];
+	// copy currents from first alpha-plane to last
+	for (pos[0]=0;pos[0]<numLines[0]-1;++pos[0])
+	{
+		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->GetCurr(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]);
+		}
+	}
+}

FDTD/engine_ext_cylinder.h

@@ -19,25 +19,27 @@
 #define ENGINE_CYLINDER_H
 
 #include "engine.h"
+#include "engine_extension.h"
 #include "operator_cylinder.h"
 
-class Engine_Cylinder : public Engine
+class Operator_Ext_Cylinder;
+
+class Engine_Ext_Cylinder : public Engine_Extension
 {
 public:
-	static Engine_Cylinder* New(const Operator_Cylinder* op);
+	Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext);
 
-	//!Iterate a number of timesteps
+	virtual void Apply2Voltages();
-	virtual bool IterateTS(unsigned int iterTS);
+
+	virtual void Apply2Current();
 
 protected:
-	Engine_Cylinder(const Operator_Cylinder* op);
+	Operator_Ext_Cylinder* cyl_Op;
 
-	virtual inline void CloseAlphaVoltages();
+	unsigned int numLines[3];
-	virtual inline void CloseAlphaCurrents();
 
-	virtual inline void R0IncludeVoltages();
+	bool CC_closedAlpha;
-
+	bool CC_R0_included;
-	const Operator_Cylinder* cyl_Op;
 };
 
 #endif // ENGINE_CYLINDER_H

FDTD/operator.h

@@ -28,6 +28,7 @@ class Operator_Extension;
 //! Abstract base-class for the FDTD-operator
 class Operator
 {
+	friend class Engine;
 public:
 	//! Create a new operator
 	static Operator* New();

FDTD/operator_cylinder.cpp

@@ -17,6 +17,7 @@
 
 #include "operator_cylinder.h"
 #include "operator_extension.h"
+#include "operator_ext_cylinder.h"
 
 Operator_Cylinder* Operator_Cylinder::New()
 {
@@ -38,26 +39,20 @@ void Operator_Cylinder::Init()
 {
 	CC_closedAlpha = false;
 	CC_R0_included = false;
-	vv_R0 = NULL;
-	vi_R0 = NULL;
 	Operator::Init();
 }
 
 void Operator_Cylinder::Reset()
 {
 	Operator::Reset();
-	delete[] vv_R0;vv_R0=NULL;
-	delete[] vi_R0;vi_R0=NULL;
 }
 
 void Operator_Cylinder::InitOperator()
 {
-	if (CC_R0_included)
-	{
-		vv_R0 = new FDTD_FLOAT[numLines[2]];
-		vi_R0 = new FDTD_FLOAT[numLines[2]];
-	}
 	Operator::InitOperator();
+
+	if (CC_closedAlpha || CC_R0_included)
+		this->AddExtension(new Operator_Ext_Cylinder(this));
 }
 
 inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const
@@ -123,43 +118,12 @@ bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
 	else if (discLines[0][0]==0.0)
 	{
 		cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl;
-		CC_R0_included=true;
+		CC_R0_included= true;  //also needed for correct ec-calculation
 	}
 
 	return true;
 }
 
-int Operator_Cylinder::CalcECOperator()
-{
-	int val = Operator::CalcECOperator();
-	if (val)
-		return val;
-
-	//if r=0 is not included -> obviously no special treatment for r=0
-	//if alpha direction is not closed, PEC-BC at r=0 necessary and already set...
-	if ((CC_R0_included==false) || (CC_closedAlpha==false))
-		return val;
-
-	unsigned int pos[3];
-	double inEC[4];
-	pos[0]=0;
-	for (pos[2]=0;pos[2]<numLines[2];++pos[2])
-	{
-		double C=0;
-		double G=0;
-		for (pos[1]=0;pos[1]<numLines[1]-1;++pos[1])
-		{
-			Calc_ECPos(2,pos,inEC);
-			C+=inEC[0]*0.5;
-			G+=inEC[1]*0.5;
-		}
-		vv[2][0][0][pos[2]] = 1;
-		vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C);
-		vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C);
-	}
-	return 0;
-}
-
 void Operator_Cylinder::ApplyElectricBC(bool* dirs)
 {
 	if (dirs==NULL) return;

FDTD/operator_cylinder.h

@@ -28,7 +28,7 @@ public:
 	
 	virtual bool SetGeometryCSX(ContinuousStructure* geo);
 
-	virtual int CalcECOperator();
+//	virtual int CalcECOperator();
 
 	virtual void ApplyElectricBC(bool* dirs);
 	virtual void ApplyMagneticBC(bool* dirs);
@@ -46,6 +46,9 @@ public:
 
 	virtual void AddExtension(Operator_Extension* op_ext);
 
+	virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
+	virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
+
 protected:
 	Operator_Cylinder();
 	virtual void Init();
@@ -55,14 +58,6 @@ protected:
 	bool CC_closedAlpha;
 	bool CC_R0_included;
 
-	virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
-	virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
-
-public:
-	//special EC operator for R0
-	FDTD_FLOAT* vv_R0; //calc new voltage from old voltage
-	FDTD_FLOAT* vi_R0; //calc new voltage from old current
-
 };
 
 #endif // OPERATOR_CYLINDER_H

FDTD/operator_ext_cylinder.cpp

@@ -0,0 +1,60 @@
+#include "operator_ext_cylinder.h"
+#include "operator_cylinder.h"
+#include "engine_ext_cylinder.h"
+
+Operator_Ext_Cylinder::Operator_Ext_Cylinder(Operator_Cylinder* op) : Operator_Extension(op)
+{
+	m_Op_Cyl = op;
+
+	CC_R0_included=m_Op_Cyl->GetR0Included();
+	CC_closedAlpha=m_Op_Cyl->GetClosedAlpha();
+
+	vv_R0 = NULL;
+	vi_R0 = NULL;
+}
+
+Operator_Ext_Cylinder::~Operator_Ext_Cylinder()
+{
+	delete[] vv_R0;vv_R0=NULL;
+	delete[] vi_R0;vi_R0=NULL;
+}
+
+bool Operator_Ext_Cylinder::BuildExtension()
+{
+	delete[] vv_R0;vv_R0=NULL;
+	delete[] vi_R0;vi_R0=NULL;
+
+	//if r=0 is not included -> obviously no special treatment for r=0
+	//if alpha direction is not closed, PEC-BC at r=0 necessary and already set...
+	if ((CC_R0_included==false) || (CC_closedAlpha==false))
+		return true;
+
+	vv_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2)];
+	vi_R0 = new FDTD_FLOAT[m_Op->GetNumberOfLines(2)];
+
+	unsigned int pos[3];
+	double inEC[4];
+	double dT = m_Op->GetTimestep();
+	pos[0]=0;
+	for (pos[2]=0;pos[2]<m_Op->GetNumberOfLines(2);++pos[2])
+	{
+		double C=0;
+		double G=0;
+		for (pos[1]=0;pos[1]<m_Op->GetNumberOfLines(1)-1;++pos[1])
+		{
+			m_Op_Cyl->Calc_ECPos(2,pos,inEC);
+			C+=inEC[0]*0.5;
+			G+=inEC[1]*0.5;
+		}
+		m_Op->GetVV(2,0,0,pos[2]) = 1;
+		vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C);
+		vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C);
+	}
+	return true;
+}
+
+Engine_Extension* Operator_Ext_Cylinder::CreateEngineExtention()
+{
+	Engine_Ext_Cylinder* eng_ext = new Engine_Ext_Cylinder(this);
+	return eng_ext;
+}

FDTD/operator_ext_cylinder.h

@@ -0,0 +1,36 @@
+#ifndef OPERATOR_EXT_CYLINDER_H
+#define OPERATOR_EXT_CYLINDER_H
+
+#include "operator_extension.h"
+#include "operator.h"
+
+class Operator_Cylinder;
+
+class Operator_Ext_Cylinder : public Operator_Extension
+{
+	friend class Engine_Ext_Cylinder;
+public:
+	Operator_Ext_Cylinder(Operator_Cylinder* op);
+	~Operator_Ext_Cylinder();
+
+	virtual bool BuildExtension();
+
+	virtual Engine_Extension* CreateEngineExtention();
+
+	virtual bool IsCylinderCoordsSave() {return true;}
+
+	virtual std::string GetExtensionName() {return std::string("Extension for the Cylinder-Coords Operator");}
+
+protected:
+	Operator_Cylinder* m_Op_Cyl;
+
+	bool CC_closedAlpha;
+	bool CC_R0_included;
+
+	//special EC operator for R0
+	FDTD_FLOAT* vv_R0; //calc new voltage from old voltage
+	FDTD_FLOAT* vi_R0; //calc new voltage from old current
+
+};
+
+#endif // OPERATOR_EXT_CYLINDER_H

openEMS.pro

@@ -44,14 +44,15 @@ SOURCES += main.cpp \
     openems.cpp \
     FDTD/engine_multithread.cpp \
     FDTD/operator_cylinder.cpp \
-    FDTD/engine_cylinder.cpp \
     FDTD/engine_sse.cpp \
     FDTD/operator_sse.cpp \
     FDTD/operator_extension.cpp \
     FDTD/engine_extension.cpp \
     FDTD/engine_ext_mur_abc.cpp \
     FDTD/operator_ext_mur_abc.cpp \
-    FDTD/excitation.cpp
+    FDTD/excitation.cpp \
+    FDTD/operator_ext_cylinder.cpp \
+    FDTD/engine_ext_cylinder.cpp
 HEADERS += tools/ErrorMsg.h \
     tools/AdrOp.h \
     tools/constants.h \
@@ -68,12 +69,13 @@ HEADERS += tools/ErrorMsg.h \
     FDTD/operator_cylinder.h \
     FDTD/engine_sse.h \
     FDTD/operator_sse.h \
-    FDTD/engine_cylinder.h \
     FDTD/operator_extension.h \
     FDTD/engine_extension.h \
     FDTD/engine_ext_mur_abc.h \
     FDTD/operator_ext_mur_abc.h \
-    FDTD/excitation.h
+    FDTD/excitation.h \
+    FDTD/operator_ext_cylinder.h \
+    FDTD/engine_ext_cylinder.h
 QMAKE_CXXFLAGS_RELEASE = -O3 \
     -g \
 	-march=native

openems.cpp

@@ -19,7 +19,7 @@
 #include <iomanip>
 #include "tools/array_ops.h"
 #include "FDTD/engine.h"
-#include "FDTD/engine_cylinder.h"
+#include "FDTD/operator_cylinder.h"
 #include "FDTD/engine_multithread.h"
 #include "FDTD/engine_sse.h"
 #include "FDTD/operator_ext_mur_abc.h"
@@ -255,7 +255,17 @@ int openEMS::SetupFDTD(const char* file)
 	if (CylinderCoords)
 	{
 		cerr << "openEMS: creating cylinder coordinate FDTD engine..." << endl;
-		FDTD_Eng = Engine_Cylinder::New((Operator_Cylinder*)FDTD_Op);
+		switch (m_engine) {
+//		case EngineType_Multithreaded:
+//			FDTD_Eng = Engine_Multithread::New(FDTD_Op,m_engine_numThreads);
+//			break;
+//		case EngineType_SSE:
+//			FDTD_Eng = Engine_sse::New(dynamic_cast<Operator_sse*>(FDTD_Op));
+//			break;
+		default:
+			FDTD_Eng = Engine::New(FDTD_Op);
+			break;
+		}
 	}
 	else
 	{