NEW: first implementation of Cylindrical Coord FDTD

todo:
  - allow closed alpha mesh
  - allow r=0 case
  - adept field processing
  - ...

FDTD/engine.cpp

@@ -20,7 +20,7 @@
 
 //! \brief construct an Engine instance
 //! it's the responsibility of the caller to free the returned pointer
-Engine* Engine::createEngine(const Operator* op)
+Engine* Engine::New(const Operator* op)
 {
 	Engine* e = new Engine(op);
 	e->Init();
@@ -52,14 +52,11 @@ void Engine::Reset()
 	curr=NULL;
 }
 
-bool Engine::IterateTS(unsigned int iterTS)
+inline void Engine::UpdateVoltages()
 {
 	unsigned int pos[3];
-	int exc_pos;
 	bool shift[3];
 
-	for (unsigned int iter=0;iter<iterTS;++iter)
-	{
 	//voltage updates
 	for (pos[0]=0;pos[0]<Op->numLines[0];++pos[0])
 	{
@@ -79,13 +76,17 @@ bool Engine::IterateTS(unsigned int iterTS)
 				volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]];
 				volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-shift[0]][pos[1]][pos[2]]);
 
-					//for x
+				//for z
 				volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]];
 				volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-shift[1]][pos[2]]);
 			}
 		}
 	}
+}
 
+inline void Engine::ApplyVoltageExcite()
+{
+	int exc_pos;
 	//soft voltage excitation here (E-field excite)
 	for (unsigned int n=0;n<Op->E_Exc_Count;++n)
 	{
@@ -94,8 +95,11 @@ bool Engine::IterateTS(unsigned int iterTS)
 //			if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl;
 		volt[Op->E_Exc_dir[n]][Op->E_Exc_index[0][n]][Op->E_Exc_index[1][n]][Op->E_Exc_index[2][n]] += Op->E_Exc_amp[n]*Op->ExciteSignal[exc_pos];
 	}
+}
 
-		//current updates
+inline void Engine::UpdateCurrents()
+{
+	unsigned int pos[3];
 	for (pos[0]=0;pos[0]<Op->numLines[0]-1;++pos[0])
 	{
 		for (pos[1]=0;pos[1]<Op->numLines[1]-1;++pos[1])
@@ -111,14 +115,27 @@ bool Engine::IterateTS(unsigned int iterTS)
 				curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]];
 				curr[1][pos[0]][pos[1]][pos[2]] += Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]+1] - volt[2][pos[0]][pos[1]][pos[2]] + volt[2][pos[0]+1][pos[1]][pos[2]]);
 
-					//for x
+				//for z
 				curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]];
 				curr[2][pos[0]][pos[1]][pos[2]] += Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]+1][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]] + volt[0][pos[0]][pos[1]+1][pos[2]]);
 			}
 		}
 	}
+}
 
-		//soft current excitation here (H-field excite)
+inline void Engine::ApplyCurrentExcite()
+{
+	int exc_pos;
+}
+
+bool Engine::IterateTS(unsigned int iterTS)
+{
+	for (unsigned int iter=0;iter<iterTS;++iter)
+	{
+		UpdateVoltages();
+		ApplyVoltageExcite();
+		UpdateCurrents();
+		ApplyCurrentExcite();
 		++numTS;
 	}
 	return true;

FDTD/engine.h

@@ -23,7 +23,7 @@
 class Engine
 {
 public:
-	static Engine* createEngine(const Operator* op);
+	static Engine* New(const Operator* op);
 	virtual ~Engine();
 
 	virtual void Init();
@@ -41,6 +41,11 @@ protected:
 	Engine(const Operator* op);
 	const Operator* Op;
 
+	virtual inline void UpdateVoltages();
+	virtual inline void ApplyVoltageExcite();
+	virtual inline void UpdateCurrents();
+	virtual inline void ApplyCurrentExcite();
+
 	FDTD_FLOAT**** volt;
 	FDTD_FLOAT**** curr;
 	unsigned int numTS;

FDTD/engine_cylinder.cpp

@@ -0,0 +1,44 @@
+/*
+*	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"
+
+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)
+{
+}
+
+Engine_Cylinder::~Engine_Cylinder()
+{
+}
+
+void Engine_Cylinder::Init()
+{
+	cerr << "Engine_Cylinder::Init()" << endl;
+	Engine::Init();
+}
+
+void Engine_Cylinder::Reset()
+{
+	Engine::Reset();
+}

FDTD/engine_cylinder.h

@@ -0,0 +1,37 @@
+/*
+*	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_CYLINDER_H
+#define ENGINE_CYLINDER_H
+
+#include "engine.h"
+#include "operator_cylinder.h"
+
+class Engine_Cylinder : public Engine
+{
+public:
+	static Engine_Cylinder* New(const Operator_Cylinder* op);
+
+	virtual void Init();
+	virtual void Reset();
+
+protected:
+	Engine_Cylinder(const Operator_Cylinder* op);
+	~Engine_Cylinder();
+};
+
+#endif // ENGINE_CYLINDER_H

FDTD/engine_multithread.cpp

@@ -34,7 +34,7 @@
 
 //! \brief construct an Engine_Multithread instance
 //! it's the responsibility of the caller to free the returned pointer
-Engine_Multithread* Engine_Multithread::createEngine(const Operator* op, unsigned int numThreads)
+Engine_Multithread* Engine_Multithread::New(const Operator* op, unsigned int numThreads)
 {
 	Engine_Multithread* e = new Engine_Multithread(op);
 	e->setNumThreads( numThreads );

FDTD/engine_multithread.h

@@ -73,7 +73,7 @@ class Engine_Multithread : public Engine
 	friend class NS_Engine_Multithread::thread;
 	friend class NS_Engine_Multithread::thread_e_excitation;
 public:
-	static Engine_Multithread* createEngine(const Operator* op, unsigned int numThreads = 0);
+	static Engine_Multithread* New(const Operator* op, unsigned int numThreads = 0);
 	virtual ~Engine_Multithread();
 
 	virtual void setNumThreads( unsigned int numThreads );

FDTD/operator.cpp

@@ -21,14 +21,20 @@
 #include "tools/array_ops.h"
 #include "fparser.hh"
 
+Operator* Operator::New()
+{
+	Operator* op = new Operator();
+	op->Init();
+	return op;
+}
+
 Operator::Operator()
 {
-	Operator::Init();
 }
 
 Operator::~Operator()
 {
-	Operator::Reset();
+	Reset();
 }
 
 void Operator::Init()
@@ -86,7 +92,7 @@ void Operator::Reset()
 		delete[] EC_R[n];
 	}
 
-	Operator::Init();
+	Init();
 }
 
 unsigned int Operator::CalcNyquistNum(double fmax)
@@ -449,6 +455,16 @@ void Operator::InitOperator()
 	ii = Create_N_3DArray(numLines);
 }
 
+inline void Operator::Calc_ECOperatorPos(int n, unsigned int* pos)
+{
+	unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]);
+	vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
+	vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
+
+	ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
+	iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
+}
+
 int Operator::CalcECOperator()
 {
 	if (Calc_EC()==0)
@@ -469,12 +485,7 @@ int Operator::CalcECOperator()
 			{
 				for (pos[2]=0;pos[2]<numLines[2];++pos[2])
 				{
-					i = MainOp->SetPos(pos[0],pos[1],pos[2]);
+					Calc_ECOperatorPos(n,pos);
-					vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
-					vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
-
-					ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
-					iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
 				}
 			}
 		}

FDTD/operator.h

@@ -28,7 +28,8 @@
 class Operator
 {
 public:
-    Operator();
+	//! Create a new operator
+	static Operator* New();
 	virtual ~Operator();
 
 	virtual bool SetGeometryCSX(ContinuousStructure* geo);
@@ -69,6 +70,9 @@ public:
 	bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false);
 
 protected:
+	//! use New() for creating a new Operator
+	Operator();
+
 	virtual void Init();
 	virtual void InitOperator();
 
@@ -92,10 +96,13 @@ protected:
 	double dT; //FDTD timestep!
 	unsigned int m_nyquistTS;
 
+	//! Calc operator at certain pos
+	virtual inline void Calc_ECOperatorPos(int n, unsigned int* pos);
+
 	//EC elements, internal only!
-	bool Calc_EC();
+	virtual bool Calc_EC();
-	bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
+	virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
-	bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
+	virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
 	double* EC_C[3];
 	double* EC_G[3];
 	double* EC_L[3];

FDTD/operator_cylinder.cpp

@@ -0,0 +1,376 @@
+/*
+*	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_cylinder.h"
+
+Operator_Cylinder* Operator_Cylinder::New()
+{
+	Operator_Cylinder* op = new Operator_Cylinder();
+	op->Init();
+	return op;
+}
+
+Operator_Cylinder::Operator_Cylinder()
+{
+}
+
+Operator_Cylinder::~Operator_Cylinder()
+{
+	Operator::Reset();
+}
+
+void Operator_Cylinder::Init()
+{
+	CC_closedAlpha = false;
+	CC_R0_included = false;
+	Operator::Init();
+}
+
+void Operator_Cylinder::Reset()
+{
+	Operator::Reset();
+}
+
+bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
+{
+	if (Operator::SetGeometryCSX(geo)==false) return false;
+
+	double minmaxA = fabs(discLines[1][numLines[1]-1]-discLines[1][0]);
+//			cerr << minmaxA -2*PI << " < " << (2*PI)/10/numLines[1] << endl;
+	if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width
+	{
+		CC_closedAlpha = true;
+		--numLines[1];
+		cout << "Operator_Cylinder::SetGeometryCSX: Alpha is a full 2*PI => closed Cylinder..." << endl;
+		cerr << "Operator_Cylinder::SetGeometryCSX: closed cylinder not yet implemented... exit... " << endl; exit(1);
+	}
+	else if (minmaxA>2*PI)
+		{cerr << "Operator_Cylinder::SetGeometryCSX: Alpha Max-Min must not be larger than 2*PI!!!" << endl; Reset(); return false;}
+	else CC_closedAlpha=false;
+
+	if (discLines[0][0]<0)
+		{cerr << "Operator_Cylinder::SetGeometryCSX: r<0 not allowed in Cylinder Coordinates!!!" << endl; Reset(); return false;}
+	else if (discLines[0][0]==0.0)
+	{
+		cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl;
+		cerr << "Operator_Cylinder::SetGeometryCSX: r=0 included not yet implemented... exit... " << endl; exit(1);
+		CC_R0_included=true;
+	}
+
+	return true;
+}
+
+inline void Operator_Cylinder::Calc_ECOperatorPos(int n, unsigned int* pos)
+{
+	unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]);
+	if (EC_C[n][i]>0)
+	{
+		vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
+		vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
+	}
+	else
+	{
+		vv[n][pos[0]][pos[1]][pos[2]] = 0;
+		vi[n][pos[0]][pos[1]][pos[2]] = 0;
+	}
+
+	if (EC_L[n][i]>0)
+	{
+		ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
+		iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
+	}
+	else
+	{
+		ii[n][pos[0]][pos[1]][pos[2]] = 0;
+		iv[n][pos[0]][pos[1]][pos[2]] = 0;
+	}
+}
+
+void Operator_Cylinder::ApplyElectricBC(bool* dirs)
+{
+	if (dirs==NULL) return;
+	if (CC_closedAlpha)
+	{
+		dirs[2]=0;dirs[3]=0; //no PEC in alpha directions...
+	}
+	if (CC_R0_included)
+	{
+		dirs[2]=0;  //no PEC in r_min directions...
+	}
+
+	Operator::ApplyElectricBC(dirs);
+}
+
+void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
+{
+	if (dirs==NULL) return;
+	if (CC_closedAlpha)
+	{
+		dirs[2]=0;dirs[3]=0; //no PMC in alpha directions...
+	}
+	if (CC_R0_included)
+	{
+		dirs[2]=0;  //no PMC in r_min directions...
+	}
+	Operator::ApplyMagneticBC(dirs);
+}
+
+bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
+{
+	double coord[3];
+	double shiftCoord[3];
+	int nP = (n+1)%3;
+	int nPP = (n+2)%3;
+	coord[0] = discLines[0][pos[0]];
+	coord[1] = discLines[1][pos[1]];
+	coord[2] = discLines[2][pos[2]];
+	unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+	double delta=MainOp->GetIndexDelta(n,pos[n]);
+	double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
+	double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
+	double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
+	double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
+	double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
+	double geom_factor,A_n;
+
+	//******************************* epsilon,kappa averaging *****************************//
+	//shift up-right
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]+deltaP*0.25;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
+	CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	switch (n)
+	{
+	case 0:
+		geom_factor = fabs(deltaPP*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		break;
+	case 1:
+		geom_factor = fabs(deltaP*deltaPP/(delta*coord[0]))*0.25;
+		break;
+	case 2:
+		geom_factor = fabs(deltaPP/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0))*0.25;
+		break;
+	}
+	geom_factor*=gridDelta;
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] = mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__;
+		inEC[1] = mat->GetKappaWeighted(n,shiftCoord)*geom_factor;
+	}
+	else
+	{
+		inEC[0] = 1*geom_factor*__EPS0__;
+		inEC[1] = 0;
+	}
+
+	//shift up-left
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	switch (n)
+	{
+	case 0:
+		geom_factor = fabs(deltaPP*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		break;
+	case 1:
+		geom_factor = fabs(deltaP_M*deltaPP/(delta*coord[0]))*0.25;
+		break;
+	case 2:
+		geom_factor = fabs(deltaPP/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0))*0.25;
+		break;
+	}
+	geom_factor*=gridDelta;
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__;
+		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*geom_factor;
+	}
+	else
+	{
+		inEC[0] += 1*geom_factor*__EPS0__;
+		inEC[1] += 0;
+	}
+
+	//shift down-right
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]+deltaP*0.25;
+	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	switch (n)
+	{
+	case 0:
+		geom_factor = fabs(deltaPP_M*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		break;
+	case 1:
+		geom_factor = fabs(deltaP*deltaPP_M/(delta*coord[0]))*0.25;
+		break;
+	case 2:
+		geom_factor = fabs(deltaPP_M/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0))*0.25;
+		break;
+	}
+	geom_factor*=gridDelta;
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__;
+		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*geom_factor;
+	}
+	else
+	{
+		inEC[0] += 1*geom_factor*__EPS0__;
+		inEC[1] += 0;
+	}
+
+	//shift down-left
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
+	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	switch (n)
+	{
+	case 0:
+		geom_factor = fabs(deltaPP_M*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		break;
+	case 1:
+		geom_factor = fabs(deltaP_M*deltaPP_M/(delta*coord[0]))*0.25;
+		break;
+	case 2:
+		geom_factor = fabs(deltaPP_M/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0))*0.25;
+		break;
+	}
+	geom_factor*=gridDelta;
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__;
+		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*geom_factor;
+	}
+	else
+	{
+		inEC[0] += 1*geom_factor*__EPS0__;
+		inEC[1] += 0;
+	}
+
+	if (CC_R0_included && (n>0) && (coord[0]==0))
+	{
+		inEC[0]=0;
+		inEC[1]=0;
+	}
+
+//	if ((n==2) && (pos[1]==0) && (pos[2]==0))
+//	cerr << n << " -> " <<  pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[0] << endl;
+
+	//******************************* mu,sigma averaging *****************************//
+	//shift down
+	shiftCoord[n] = coord[n]-delta_M*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*0.5;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	double delta_n = fabs(delta_M);
+	if (n==1)
+	{
+		delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius
+	}
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[2] = delta_n / mat->GetMueWeighted(n,shiftCoord);
+		if (mat->GetSigma(n))
+			inEC[3] = delta_n / mat->GetSigmaWeighted(n,shiftCoord);
+		else
+			inEC[3] = 0;
+	}
+	else
+	{
+		inEC[2] = delta_n;
+		inEC[3] = 0;
+	}
+	//shift up
+	shiftCoord[n] = coord[n]+delta*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*0.5;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	delta_n = fabs(delta);
+	if (n==1)
+	{
+		delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius
+	}
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[2] += mat->GetMue(n)*delta_n;
+		if (mat->GetSigmaWeighted(n,shiftCoord))
+			inEC[3] += delta_n/mat->GetSigmaWeighted(n,shiftCoord);
+		else
+			inEC[3] = 0;
+	}
+	else
+	{
+		inEC[2] += 1*delta_n;
+		inEC[3] = 0;
+	}
+
+	A_n = fabs(deltaP*deltaPP);
+	if (n==0) //z-direction n==0 -> r; nP==1 -> alpha; nPP==2 -> z
+	{
+		A_n = A_n * coord[0];
+	}
+	if (n==2) //z-direction n==2 -> z; nP==0 -> r; nPP==1 -> alpha
+	{
+		A_n = fabs(deltaPP) * (pow(coord[0]+fabs(deltaP),2.0) - pow(coord[0],2.0))*0.5;
+	}
+
+	inEC[2] = gridDelta * A_n * 2 * __MUE0__ / inEC[2];
+	if (inEC[3]) inEC[3]=gridDelta * A_n * 2 / inEC[3];
+
+//	if ((n==0) && (pos[1]==0) && (pos[2]==0))
+//		cerr <<  inEC[2]/(coord[0]) <<  endl;
+//	cerr << n << " -> " <<  pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[2] << endl;
+
+	return true;
+}
+
+bool Operator_Cylinder::Calc_EffMatPos(int n, unsigned int* pos, double* inMat)
+{
+	return false; //fixme
+
+//	int nP = (n+1)%3;
+//	int nPP = (n+2)%3;
+//
+//	unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+//	double delta=MainOp->GetIndexDelta(n,pos[n]);
+//	double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
+//	double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
+//
+//	double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
+//	double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
+//	double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
+//
+//	this->Calc_ECPos(n,pos,inMat);
+//
+//	inMat[0] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
+//	inMat[1] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
+//
+//	inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
+//	inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
+//
+//	return true;
+}
+

FDTD/operator_cylinder.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_CYLINDER_H
+#define OPERATOR_CYLINDER_H
+
+#include "operator.h"
+
+class Operator_Cylinder : public Operator
+{
+public:
+	static Operator_Cylinder* New();
+	virtual ~Operator_Cylinder();
+	
+	virtual bool SetGeometryCSX(ContinuousStructure* geo);
+
+	virtual void ApplyElectricBC(bool* dirs);
+	virtual void ApplyMagneticBC(bool* dirs);
+
+	virtual void Reset();
+
+protected:
+	Operator_Cylinder();
+	virtual void Init();
+
+	virtual inline void Calc_ECOperatorPos(int n, unsigned int* pos);
+
+	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);
+};
+
+#endif // OPERATOR_CYLINDER_H

FDTD/processfields.cpp

@@ -81,7 +81,7 @@ void ProcessFields::InitProcess()
 			H5::DataSet dataset = group->createDataSet( names[n].c_str(), datatype, dataspace );
 			//convert to float...
 			float* array = new float[NrLines[n]];
-			for (int i=0;i<NrLines[n];++i)
+			for (unsigned int i=0;i<NrLines[n];++i)
 				array[i] = Lines[n][i];
 			//write to dataset
 			dataset.write( array, H5::PredType::NATIVE_FLOAT );

matlab/InitCylindricalFDTD.m

@@ -0,0 +1,16 @@
+function FDTD = InitCylindricalFDTD(NrTS, endCrit, varargin)
+% function FDTD = InitCylindricalFDTD(NrTS, endCrit, varargin)
+%
+% see also InitFDTD
+%
+% e.g FDTD = InitCylindricalFDTD(5e5,1e-6,'OverSampling',10)
+%
+% openEMS matlab interface
+% -----------------------
+% author: Thorsten Liebig
+
+FDTD = InitFDTD(NrTS, endCrit, varargin);
+
+FDTD.ATTRIBUTE.CylinderCoords=1;
+
+

openEMS.pro

@@ -1,7 +1,8 @@
 # -------------------------------------------------
 # Project created by QtCreator 2010-02-26T22:34:51
 # -------------------------------------------------
-QT -= gui core
+QT -= gui \
+    core
 TARGET = openEMS
 CONFIG += console
 CONFIG -= app_bundle
@@ -17,11 +18,11 @@ LIBS += -L../CSXCAD \
     -L../tinyxml \
     -ltinyxml \
     -lboost_thread \
-	-lhdf5 -lhdf5_cpp
+    -lhdf5 \
+    -lhdf5_cpp
 QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../CSXCAD\'
 QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../fparser\'
 QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../tinyxml\'
-
 SOURCES += main.cpp \
     tools/ErrorMsg.cpp \
     tools/AdrOp.cpp \
@@ -35,7 +36,9 @@ SOURCES += main.cpp \
     FDTD/processcurrent.cpp \
     examples/FDTD_examples.cpp \
     openems.cpp \
-    FDTD/engine_multithread.cpp
+    FDTD/engine_multithread.cpp \
+    FDTD/operator_cylinder.cpp \
+    FDTD/engine_cylinder.cpp
 HEADERS += tools/ErrorMsg.h \
     tools/AdrOp.h \
     tools/constants.h \
@@ -49,7 +52,12 @@ HEADERS += tools/ErrorMsg.h \
     FDTD/processcurrent.h \
     examples/FDTD_examples.h \
     openems.h \
-    FDTD/engine_multithread.h
+    FDTD/engine_multithread.h \
-
+    FDTD/operator_cylinder.h \
-QMAKE_CXXFLAGS_RELEASE = -O2 -g -march=native
+    FDTD/engine_cylinder.h
-QMAKE_CXXFLAGS_DEBUG = -O0 -g -march=native
+QMAKE_CXXFLAGS_RELEASE = -O2 \
+    -g \
+    -march=native
+QMAKE_CXXFLAGS_DEBUG = -O0 \
+    -g \
+    -march=native

openems.cpp

@@ -18,8 +18,8 @@
 #include "openems.h"
 #include <iomanip>
 #include "tools/array_ops.h"
-#include "FDTD/operator.h"
 #include "FDTD/engine.h"
+#include "FDTD/engine_cylinder.h"
 #include "FDTD/engine_multithread.h"
 #include "FDTD/processvoltage.h"
 #include "FDTD/processcurrent.h"
@@ -43,6 +43,7 @@ openEMS::openEMS()
 	FDTD_Op=NULL;
 	FDTD_Eng=NULL;
 	PA=NULL;
+	CylinderCoords = false;
 	Enable_Dumps = true;
 	DebugMat = false;
 	DebugOp = false;
@@ -190,6 +191,13 @@ int openEMS::SetupFDTD(const char* file)
 	else
 		NrTS = help;
 
+	FDTD_Opts->QueryIntAttribute("CylinderCoords",&help);
+	if (help==1)
+	{
+		cout << "Using a cylinder coordinate FDTD..." << endl;
+		CylinderCoords = true;
+	}
+
 	FDTD_Opts->QueryDoubleAttribute("endCriteria",&endCrit);
 	if (endCrit==0)
 		endCrit=1e-6;
@@ -226,8 +234,16 @@ int openEMS::SetupFDTD(const char* file)
 
 	//*************** setup operator ************//
 	cout << "Create Operator..." << endl;
-	FDTD_Op = new Operator();
+	if (CylinderCoords)
-	if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1);
+	{
+		FDTD_Op = Operator_Cylinder::New();
+	}
+	else
+	{
+		FDTD_Op = Operator::New();
+	}
+
+	if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(2);
 
 	FDTD_Op->CalcECOperator();
 
@@ -251,14 +267,22 @@ int openEMS::SetupFDTD(const char* file)
 	cout << "Creation time for operator: " << difftime(OpDoneTime,startTime) << " s" << endl;
 
 	//create FDTD engine
+	if (CylinderCoords)
+	{
+		cerr << "openEMS: creating cylinder coordinate FDTD engine..." << endl;
+		FDTD_Eng = Engine_Cylinder::New((Operator_Cylinder*)FDTD_Op);
+	}
+	else
+	{
 		switch (m_engine) {
 		case EngineType_Multithreaded:
-		FDTD_Eng = Engine_Multithread::createEngine(FDTD_Op,m_engine_numThreads);
+			FDTD_Eng = Engine_Multithread::New(FDTD_Op,m_engine_numThreads);
 			break;
 		default:
-		FDTD_Eng = Engine::createEngine(FDTD_Op);
+			FDTD_Eng = Engine::New(FDTD_Op);
 			break;
 		}
+	}
 
 	time_t currTime = time(NULL);
 

openems.h

@@ -45,6 +45,9 @@ public:
 
 protected:
 	void SetupExcitation(TiXmlElement* Excite);
+
+	bool CylinderCoords;
+
 	//! Number of Timesteps
 	unsigned int NrTS;
 	bool Enable_Dumps;