Operator & Operator_Cylinder: changes to material averaging methods

Operator:
- new method to setup the mesh: SetupCSXGrid()
- Most methods handling material and operator calculations are now protected.
- New method for accessing the material distribution.

Operator_Cylinder:
- overloaded SetupCSXGrid() handling cylindrical specialties
- This commit adds an additional line in alpha-direction
in case of a closed cylinder simulation. Thereby the material averaging
will be accurate at the alpha-interface

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

FDTD/extensions/engine_ext_cylinder.cpp

@@ -61,14 +61,13 @@ void Engine_Ext_Cylinder::Apply2Voltages()
 
 	//close alpha
 	unsigned int pos[3];
-	// copy voltages from last alpha-plane to first
+	// copy tangential voltages from last alpha-plane to first
-	unsigned int last_A_Line = numLines[1]-1;
+	unsigned int last_A_Line = numLines[1]-2;
 	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 	{
 		for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
 		{
 			m_Eng->SetVolt(0,pos[0],0,pos[2], m_Eng->GetVolt(0,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->SetVolt(2,pos[0],0,pos[2], m_Eng->GetVolt(2,pos[0],last_A_Line,pos[2]) );
 		}
 	}
@@ -80,14 +79,13 @@ void Engine_Ext_Cylinder::Apply2Current()
 
 	//close alpha
 	unsigned int pos[3];
-	// copy currents from first alpha-plane to last
+	// copy tangential 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;
+		unsigned int last_A_Line = numLines[1]-2;
 		for (pos[2]=0; pos[2]<numLines[2]-1; ++pos[2])
 		{
 			m_Eng->SetCurr(0,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(0,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->SetCurr(2,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(2,pos[0],0,pos[2]) );
 		}
 	}

FDTD/operator.cpp

@@ -56,6 +56,8 @@ Engine* Operator::CreateEngine() const
 
 void Operator::Init()
 {
+	CSX = NULL;
+
 	Operator_Base::Init();
 
 	vv=NULL;
@@ -177,6 +179,15 @@ double Operator::GetEdgeLength(int n, const unsigned int* pos, bool dualMesh) co
 	}
 }
 
+double Operator::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const
+{
+	if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) )
+		return 0.0;
+
+	unsigned int uiPos[]={pos[0],pos[1],pos[2]};
+	return GetNodeWidth(ny, uiPos, dualMesh);
+}
+
 double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const
 {
 	int nyP = (ny+1)%3;
@@ -184,6 +195,15 @@ double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) c
 	return GetNodeWidth(nyP,pos,dualMesh) * GetNodeWidth(nyPP,pos,dualMesh);
 }
 
+double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
+{
+	if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) )
+		return 0.0;
+
+	unsigned int uiPos[]={pos[0],pos[1],pos[2]};
+	return GetNodeArea(ny, uiPos, dualMesh);
+}
+
 bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool lower, bool* inside) const
 {
 	bool ok=true;
@@ -577,6 +597,35 @@ void Operator::DumpMaterial2File(string filename)
 	cout << " done!" << endl;
 }
 
+ bool Operator::SetupCSXGrid(CSRectGrid* grid)
+ {
+	 for (int n=0; n<3; ++n)
+	 {
+		 discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
+		 if (numLines[n]<3)
+		 {
+			 cerr << "CartOperator::SetupCSXGrid: you need at least 3 disc-lines in every direction (3D!)!!!" << endl;
+			 Reset();
+			 return false;
+		 }
+	 }
+	 MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
+	 MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
+	 if (grid->GetDeltaUnit()<=0)
+	 {
+		 cerr << "CartOperator::SetupCSXGrid: grid delta unit must not be <=0 !!!" << endl;
+		 Reset();
+		 return false;
+	 }
+	 else gridDelta=grid->GetDeltaUnit();
+	 MainOp->SetGridDelta(1);
+	 MainOp->AddCellAdrOp();
+
+	 //delete the grid clone...
+	 delete grid;
+	 return true;
+ }
+
 bool Operator::SetGeometryCSX(ContinuousStructure* geo)
 {
 	if (geo==NULL) return false;
@@ -584,29 +633,8 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo)
 	CSX = geo;
 
 	CSRectGrid* grid=CSX->GetGrid();
-	for (int n=0; n<3; ++n)
+
-	{
+	return SetupCSXGrid(CSRectGrid::Clone(grid));
-		discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
-		if (n==1)
-			if (numLines[n]<3)
-			{
-				cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl;
-				Reset();
-				return false;
-			}
-	}
-	MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
-	MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
-	if (grid->GetDeltaUnit()<=0)
-	{
-		cerr << "CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!!" << endl;
-		Reset();
-		return false;
-	}
-	else gridDelta=grid->GetDeltaUnit();
-	MainOp->SetGridDelta(1);
-	MainOp->AddCellAdrOp();
-	return true;
 }
 
 void Operator::InitOperator()
@@ -982,6 +1010,44 @@ double Operator::GetRawDiscDelta(int ny, const int pos) const
 	return (discLines[ny][pos+1] - discLines[ny][pos]);
 }
 
+double Operator::GetMaterial(int ny, const double* coords, int MatType, bool markAsUsed) const
+{
+	CSProperties* prop = CSX->GetPropertyByCoordPriority(coords,CSProperties::MATERIAL,markAsUsed);
+	CSPropMaterial* mat = dynamic_cast<CSPropMaterial*>(prop);
+	if (mat)
+	{
+		switch (MatType)
+		{
+		case 0:
+			return mat->GetEpsilonWeighted(ny,coords);
+		case 1:
+			return mat->GetKappaWeighted(ny,coords);
+		case 2:
+			return mat->GetMueWeighted(ny,coords);
+		case 3:
+			return mat->GetSigmaWeighted(ny,coords);
+		default:
+			cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
+			return 0;
+		}
+	}
+
+	switch (MatType)
+	{
+	case 0:
+		return 1;
+	case 1:
+		return 0;
+	case 2:
+		return 1;
+	case 3:
+		return 0;
+	default:
+		cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
+		return 0;
+	}
+}
+
 bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const
 {
 	int n=ny;
@@ -1008,23 +1074,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	shiftCoord[n] = coord[n]+delta*0.5;
 	shiftCoord[nP] = coord[nP]+deltaP*0.25;
 	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
-	A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true);
+	A_n = GetNodeArea(ny,loc_pos,true);
-//	{
+	EffMat[0] = GetMaterial(n, shiftCoord, 0)*A_n;
-//		cerr << ny << " " << pos[0] << " " <<  pos[1] << " " <<  pos[2] << ": " << A_n << endl;
+	EffMat[1] = GetMaterial(n, shiftCoord, 1)*A_n;
-//		exit(0);
-//	}
-	CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true);
-	if (prop)
-	{
-		CSPropMaterial* mat = prop->ToMaterial();
-		EffMat[0] = mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		EffMat[1] = mat->GetKappaWeighted(n,shiftCoord)*A_n;
-	}
-	else
-	{
-		EffMat[0] = 1*A_n;
-		EffMat[1] = 0;
-	}
 	area+=A_n;
 
 	//shift up-left
@@ -1033,20 +1085,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
 
 	--loc_pos[nP];
-	A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true);
+	A_n = GetNodeArea(ny,loc_pos,true);
-//	cerr << A_n << endl;
+	EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n;
-	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true);
+	EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n;
-	if (prop)
-	{
-		CSPropMaterial* mat = prop->ToMaterial();
-		EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n;
-	}
-	else
-	{
-		EffMat[0] += 1*A_n;
-		EffMat[1] += 0;
-	}
 	area+=A_n;
 
 	//shift down-right
@@ -1055,19 +1096,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
 	++loc_pos[nP];
 	--loc_pos[nPP];
-	A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true);
+	A_n = GetNodeArea(ny,loc_pos,true);
-	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true);
+	EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n;
-	if (prop)
+	EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n;
-	{
-		CSPropMaterial* mat = prop->ToMaterial();
-		EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n;
-	}
-	else
-	{
-		EffMat[0] += 1*A_n;
-		EffMat[1] += 0;
-	}
 	area+=A_n;
 
 	//shift down-left
@@ -1075,19 +1106,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
 	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
 	--loc_pos[nP];
-	A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true);
+	A_n = GetNodeArea(ny,loc_pos,true);
-	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true);
+	EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n;
-	if (prop)
+	EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n;
-	{
-		CSPropMaterial* mat = prop->ToMaterial();
-		EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n;
-	}
-	else
-	{
-		EffMat[0] += 1*A_n;
-		EffMat[1] += 0;
-	}
 	area+=A_n;
 
 	EffMat[0]*=__EPS0__/area;
@@ -1098,27 +1119,19 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	loc_pos[1]=pos[1];
 	loc_pos[2]=pos[2];
 	double length=0;
+
 	//shift down
 	shiftCoord[n] = coord[n]-delta_M*0.25;
 	shiftCoord[nP] = coord[nP]+deltaP*0.5;
 	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
 	--loc_pos[n];
-	double delta_ny = GetNodeWidth(n,(unsigned int*)loc_pos,true);
+	double delta_ny = GetNodeWidth(n,loc_pos,true);
-	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true);
+	EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2);
-	if (prop)
+	double sigma = GetMaterial(n, shiftCoord, 3);
-	{
+	if (sigma)
-		CSPropMaterial* mat = prop->ToMaterial();
+		EffMat[3] = delta_ny / sigma;
-		EffMat[2] = delta_ny / mat->GetMueWeighted(n,shiftCoord);
-		if (mat->GetSigmaWeighted(n,shiftCoord))
-			EffMat[3] = delta_ny / mat->GetSigmaWeighted(n,shiftCoord);
-		else
-			EffMat[3] = 0;
-	}
 	else
-	{
-		EffMat[2] = delta_ny;
 		EffMat[3] = 0;
-	}
 	length=delta_ny;
 
 	//shift up
@@ -1126,22 +1139,13 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
 	shiftCoord[nP] = coord[nP]+deltaP*0.5;
 	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
 	++loc_pos[n];
-	delta_ny = GetNodeWidth(n,(unsigned int*)loc_pos,true);
+	delta_ny = GetNodeWidth(n,loc_pos,true);
-	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true);
+	EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2);
-	if (prop)
+	sigma = GetMaterial(n, shiftCoord, 3);
-	{
+	if (sigma)
-		CSPropMaterial* mat = prop->ToMaterial();
+		EffMat[3] += delta_ny / sigma;
-		EffMat[2] += delta_ny / mat->GetMueWeighted(n,shiftCoord);
-		if (mat->GetSigmaWeighted(n,shiftCoord))
-			EffMat[3] += delta_ny/mat->GetSigmaWeighted(n,shiftCoord);
-		else
-			EffMat[3] = 0;
-	}
 	else
-	{
-		EffMat[2] += 1*delta_ny;
 		EffMat[3] = 0;
-	}
 	length+=delta_ny;
 
 	EffMat[2] = length * __MUE0__ / EffMat[2];

FDTD/operator.h

@@ -35,6 +35,8 @@ class Operator : public Operator_Base
 	friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice
 	friend class Operator_Ext_PML_SF_Plane;
 	friend class Operator_Ext_Excitation;
+	friend class Operator_Ext_UPML;
+	friend class Operator_Ext_Cylinder;
 public:
 	enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4};
 
@@ -48,20 +50,6 @@ public:
 
 	virtual int CalcECOperator( DebugFlags debugFlags = None );
 
-	//! Calculate the FDTD equivalent circuit parameter for the given position and direction ny. \sa Calc_EffMat_Pos
-	virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC) const;
-
-	//! Get the FDTD raw disc delta, needed by Calc_EffMatPos() \sa Calc_EffMatPos
-	/*!
-	  Get the raw disc delta for a given position and direction.
-	  The result will be positive if a disc delta inside the simulation domain is requested.
-	  The result will be the negative value of the first or last disc delta respectivly if the position is outside the field domain.
-	  */
-	virtual double GetRawDiscDelta(int ny, const int pos) 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 bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) {return Exc->setupExcitation(Excite,maxTS);};
 
 	virtual void DumpExciationSignals();
@@ -112,9 +100,13 @@ public:
 
 	//! Get the node width for a given direction \a n and a given mesh position \a pos
 	virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeLength(ny,pos,!dualMesh);}
+	//! Get the node width for a given direction \a n and a given mesh position \a pos
+	virtual double GetNodeWidth(int ny, const int pos[3], bool dualMesh = false) const;
 
 	//! Get the node area for a given direction \a n and a given mesh position \a pos
 	virtual double GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const;
+	//! Get the node area for a given direction \a n and a given mesh position \a pos
+	virtual double GetNodeArea(int ny, const int pos[3], bool dualMesh = false) const;
 
 	//! Get the length of an FDTD edge (unit is meter).
 	virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const;
@@ -147,6 +139,8 @@ protected:
 	virtual void InitDataStorage();
 	virtual void InitExcitation();
 
+	virtual bool SetupCSXGrid(CSRectGrid* grid);
+
 	struct Grid_Path
 	{
 		vector<unsigned int> posPath[3];
@@ -174,6 +168,23 @@ protected:
 	double CalcTimestep_Var1();
 	double CalcTimestep_Var3();
 
+	//! Calculate the FDTD equivalent circuit parameter for the given position and direction ny. \sa Calc_EffMat_Pos
+	virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC) const;
+
+	//! Get the FDTD raw disc delta, needed by Calc_EffMatPos() \sa Calc_EffMatPos
+	/*!
+	  Get the raw disc delta for a given position and direction.
+	  The result will be positive if a disc delta inside the simulation domain is requested.
+	  The result will be the negative value of the first or last disc delta respectivly if the position is outside the field domain.
+	  */
+	virtual double GetRawDiscDelta(int ny, const int pos) const;
+
+	//! Get the material at a given coordinate, direction and type from CSX (internal use only)
+	virtual double GetMaterial(int ny, const double coords[3], int MatType, bool markAsUsed=true) 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;
+
 	//! Calc operator at certain \a pos
 	virtual void Calc_ECOperatorPos(int n, unsigned int* pos);
 

FDTD/operator_cylinder.cpp

@@ -46,6 +46,27 @@ void Operator_Cylinder::Init()
 	Operator_Multithread::Init();
 }
 
+double Operator_Cylinder::GetRawDiscDelta(int ny, const int pos) const
+{
+	if (CC_closedAlpha && ny==1 && pos==-1)
+	{
+//		cerr << (discLines[1][numLines[1]-2] - discLines[1][numLines[1]-3]) << " vs "  << Operator_Multithread::GetRawDiscDelta(ny,pos) << endl;
+		return (discLines[1][numLines[1]-2] - discLines[1][numLines[1]-3]);
+	}
+
+	return Operator_Multithread::GetRawDiscDelta(ny,pos);
+}
+
+double Operator_Cylinder::GetMaterial(int ny, const double* coords, int MatType, bool markAsUsed) const
+{
+	double l_coords[] = {coords[0],coords[1],coords[2]};
+	if (CC_closedAlpha && (coords[1]>GetDiscLine(1,0,false)+2*PI))
+		l_coords[1]-=2*PI;
+	if (CC_closedAlpha && (coords[1]<GetDiscLine(1,0,false)))
+		l_coords[1] += 2*PI;
+	return Operator_Multithread::GetMaterial(ny,l_coords,MatType,markAsUsed);
+}
+
 int Operator_Cylinder::CalcECOperator( DebugFlags debugFlags )
 {
 	// debugs only work with the native vector dumps
@@ -61,7 +82,7 @@ inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const
 {
 	//this is necessary for a correct field processing... cylindrical engine has to reset this by adding +1
 	if (CC_closedAlpha && ny==1)
-		return Operator_Multithread::GetNumberOfLines(ny)-1;
+		return Operator_Multithread::GetNumberOfLines(ny)-2;
 
 	return Operator_Multithread::GetNumberOfLines(ny);
 }
@@ -76,16 +97,12 @@ string Operator_Cylinder::GetDirName(int ny) const
 
 bool Operator_Cylinder::GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool dualMesh) const
 {
-	bool ret = Operator_Multithread::GetYeeCoords(ny,pos,coords,dualMesh);
+	if (CC_closedAlpha && (coords[1]>GetDiscLine(1,0,false)+2*PI))
-	if ((CC_closedAlpha==false) || (ny!=1))
-		return ret;
-
-	if (coords[1]>GetDiscLine(1,0,false)+2*PI)
-	{
 		coords[1]-=2*PI;
-	}
+	if (CC_closedAlpha && (coords[1]<GetDiscLine(1,0,false)))
+		coords[1]+=2*PI;
 
-	return ret;
+	return Operator_Multithread::GetYeeCoords(ny,pos,coords,dualMesh);
 }
 
 double Operator_Cylinder::GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh) const
@@ -98,6 +115,18 @@ double Operator_Cylinder::GetNodeWidth(int ny, const unsigned int pos[3], bool d
 	return width;
 }
 
+double Operator_Cylinder::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const
+{
+	if ( (pos[0]<0) || (pos[1]<0 && CC_closedAlpha==false) || (pos[2]<0) )
+		return 0.0;
+
+	unsigned int uiPos[]={pos[0],pos[1],pos[2]};
+	if (pos[1]<0 && CC_closedAlpha==true)
+		uiPos[1]+=numLines[1]-2;
+
+	return GetNodeWidth(ny, uiPos, dualMesh);
+}
+
 double Operator_Cylinder::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const
 {
 	if (pos[ny]>=numLines[ny]) return 0.0;
@@ -127,6 +156,18 @@ double Operator_Cylinder::GetNodeArea(int ny, const unsigned int pos[3], bool du
 	return Operator_Multithread::GetNodeArea(ny,pos,dualMesh);
 }
 
+double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
+{
+	if ( (pos[0]<0) || (pos[1]<0 && CC_closedAlpha==false) || (pos[2]<0) )
+		return 0.0;
+
+	unsigned int uiPos[]={pos[0],pos[1],pos[2]};
+	if (pos[1]<0 && CC_closedAlpha==true)
+		uiPos[1]+=numLines[1]-2;
+
+	return GetNodeArea(ny, uiPos, dualMesh);
+}
+
 double Operator_Cylinder::GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh) const
 {
 	double length = Operator_Multithread::GetEdgeLength(ny,pos,dualMesh);
@@ -143,31 +184,23 @@ double Operator_Cylinder::GetEdgeArea(int ny, const unsigned int pos[3], bool du
 	return GetEdgeLength(1,pos,!dualMesh) * GetEdgeLength(2,pos,!dualMesh);
 }
 
-bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
+bool Operator_Cylinder::SetupCSXGrid(CSRectGrid* grid)
 {
-	if (Operator_Multithread::SetGeometryCSX(geo)==false) return false;
+	unsigned int alphaNum;
+	double* alphaLines = NULL;
+	alphaLines = grid->GetLines(1,alphaLines,alphaNum,true);
 
-	double minmaxA = fabs(discLines[1][numLines[1]-1]-discLines[1][0]);
+	double minmaxA = fabs(alphaLines[alphaNum-1]-alphaLines[0]);
-	if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width
+	if (fabs(minmaxA-2*PI) < OPERATOR_CYLINDER_CLOSED_ALPHA_THRESHOLD)
 	{
-		cout << "Operator_Cylinder::SetGeometryCSX: Alpha is a full 2*PI => closed Cylinder..." << endl;
+		cout << "Operator_Cylinder::SetupCSXGrid: Alpha is a full 2*PI => closed Cylinder..." << endl;
 		CC_closedAlpha = true;
-		discLines[1][numLines[1]-1] = discLines[1][0] + 2*PI;
+		grid->SetLine(1,alphaNum-1,2*PI+alphaLines[0]);
-		cerr << "Operator_Cylinder::SetGeometryCSX: Warning, not handling the disc-line width and material averaging correctly yet for a closed cylinder..." << endl;
+		grid->AddDiscLine(1,2*PI+alphaLines[1]);
-		if (MainOp->GetIndexDelta(1,0)-MainOp->GetIndexDelta(1,numLines[1]-2) > (2*PI)/10/numLines[1])
-		{
-			cerr << "Operator_Cylinder::SetGeometryCSX: first and last angle delta must be the same... deviation to large..." << MainOp->GetIndexDelta(1,0) - MainOp->GetIndexDelta(1,numLines[1]-2) << endl;
-			exit(1);
-		}
-		if (MainOp->GetIndexDelta(1,0)-MainOp->GetIndexDelta(1,numLines[1]-2) > 0)
-		{
-			cerr << "Operator_Cylinder::SetGeometryCSX: first and last angle delta must be the same... auto correction of deviation: " << MainOp->GetIndexDelta(1,0) - MainOp->GetIndexDelta(1,numLines[1]-2) << endl;
-			discLines[1][numLines[1]-2] = discLines[1][numLines[1]-1]-MainOp->GetIndexDelta(1,0);
-		}
 	}
 	else if (minmaxA>2*PI)
 	{
-		cerr << "Operator_Cylinder::SetGeometryCSX: Alpha Max-Min must not be larger than 2*PI!!!" << endl;
+		cerr << "Operator_Cylinder::SetupCSXGrid: Alpha Max-Min must not be larger than 2*PI!!!" << endl;
 		Reset();
 		return false;
 	}
@@ -176,48 +209,27 @@ bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
 		CC_closedAlpha=false;
 	}
 
-	if (discLines[0][0]<0)
+	if (grid->GetLine(0,0)<0)
 	{
-		cerr << "Operator_Cylinder::SetGeometryCSX: r<0 not allowed in Cylinder Coordinates!!!" << endl;
+		cerr << "Operator_Cylinder::SetupCSXGrid: r<0 not allowed in Cylinder Coordinates!!!" << endl;
 		Reset();
 		return false;
 	}
-	else if (discLines[0][0]==0.0)
+	else if (grid->GetLine(0,0)==0.0)
 	{
-		cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl;
+		cout << "Operator_Cylinder::SetupCSXGrid: r=0 included..." << endl;
 		CC_R0_included= true;  //also needed for correct ec-calculation
 	}
 
+	if (Operator_Multithread::SetupCSXGrid(grid)==false)
+		return false;
+
 	if (CC_closedAlpha || CC_R0_included)
 		this->AddExtension(new Operator_Ext_Cylinder(this));
 
 	return true;
 }
 
-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)
-	{
-		// E in alpha direction ( aka volt[1][x][y][z] ) is not defined for r==0 --> always zero...
-		unsigned int pos[3] = {0,0,0};
-		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
-		{
-			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
-			{
-				SetVV(1,pos[0],pos[1],pos[2], 0 );
-				SetVI(1,pos[0],pos[1],pos[2], 0 );
-			}
-		}
-	}
-	Operator_Multithread::ApplyElectricBC(dirs);
-}
-
 void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
 {
 	if (dirs==NULL) return;

FDTD/operator_cylinder.h

@@ -18,6 +18,8 @@
 #ifndef OPERATOR_CYLINDER_H
 #define OPERATOR_CYLINDER_H
 
+#define OPERATOR_CYLINDER_CLOSED_ALPHA_THRESHOLD 1e-6
+
 #include "operator_multithread.h"
 
 //! This class creates an operator for a cylindrical FDTD.
@@ -32,11 +34,6 @@ public:
 	static Operator_Cylinder* New(unsigned int numThreads = 0);
 	virtual ~Operator_Cylinder();
 
-	virtual int CalcECOperator( DebugFlags debugFlags = None );
-
-	virtual bool SetGeometryCSX(ContinuousStructure* geo);
-
-	virtual void ApplyElectricBC(bool* dirs);
 	virtual void ApplyMagneticBC(bool* dirs);
 
 	virtual unsigned int GetNumberOfLines(int ny) const;
@@ -49,9 +46,13 @@ public:
 
 	//! Get the node width for a given direction \a n and a given mesh posisition \a pos
 	virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const;
+	//! Get the node width for a given direction \a n and a given mesh posisition \a pos
+	virtual double GetNodeWidth(int ny, const int pos[3], bool dualMesh = false) const;
 
 	//! Get the node area for a given direction \a n and a given mesh posisition \a pos
 	virtual double GetNodeArea(int n, const unsigned int* pos, bool dualMesh=false) const;
+	//! Get the node area for a given direction \a n and a given mesh posisition \a pos
+	virtual double GetNodeArea(int ny, const int pos[3], bool dualMesh = false) const;
 
 	//! Get the length of an FDTD edge, including radius corrected alpha-mesh width.
 	virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const;
@@ -75,6 +76,14 @@ protected:
 	//Calc timestep only internal use
 	virtual double CalcTimestep();
 
+	virtual bool SetupCSXGrid(CSRectGrid* grid);
+
+	virtual double GetRawDiscDelta(int ny, const int pos) const;
+
+	virtual double GetMaterial(int ny, const double coords[3], int MatType, bool markAsUsed=true) const;
+
+	virtual int CalcECOperator( DebugFlags debugFlags = None );
+
 	bool CC_closedAlpha;
 	bool CC_R0_included;
 };

FDTD/operator_cylindermultigrid.cpp

@@ -60,14 +60,14 @@ double Operator_CylinderMultiGrid::GetNumberCells() const
 	return 0;
 }
 
-bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
+bool Operator_CylinderMultiGrid::SetupCSXGrid(CSRectGrid* grid)
 {
-	if (Operator_Cylinder::SetGeometryCSX(geo)==false)
+	if (Operator_Cylinder::SetupCSXGrid(grid)==false)
 		return false;
 
-	if (numLines[1]%2 != 1)
+	if ((numLines[1]-CC_closedAlpha)%2 != 1)
 	{
-		cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, number of line in alpha direction must be odd... found: " << numLines[1] << endl;
+		cerr << "Operator_CylinderMultiGrid::SetupCSXGrid: Error, number of line in alpha direction must be odd... found: " << numLines[1] << endl;
 		exit(0);
 	}
 
@@ -77,7 +77,7 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
 	{
 		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;
+			cerr << "Operator_CylinderMultiGrid::SetupCSXGrid: Error, mesh has to be homogenous in alpha direction for multi grid engine, violation found at: " << n << endl;
 			exit(0);
 		}
 	}
@@ -88,16 +88,23 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
 		if (m_Split_Rad < discLines[0][n])
 		{
 			m_Split_Pos = n;
-			cout << "Operator_CylinderMultiGrid::SetGeometryCSX: Found mesh split position @" << m_Split_Pos << endl;
+			cout << "Operator_CylinderMultiGrid::SetupCSXGrid: 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;
+		cerr << "Operator_CylinderMultiGrid::SetupCSXGrid: Error, split invalid..." << endl;
 		return false;
 	}
+	return true;
+}
+
+bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
+{
+	if (Operator_Cylinder::SetGeometryCSX(geo)==false)
+		return false;
 
 	CSRectGrid* grid = geo->GetGrid();
 

FDTD/operator_cylindermultigrid.h

@@ -43,8 +43,6 @@ public:
 
 	virtual unsigned int GetSplitPos() const {return m_Split_Pos;}
 
-	virtual int CalcECOperator( DebugFlags debugFlags = None );
-
 	virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS);
 
 	virtual void SetBoundaryCondition(int* BCs);
@@ -67,6 +65,10 @@ protected:
 	void Delete();
 	virtual void Reset();
 
+	virtual bool SetupCSXGrid(CSRectGrid* grid);
+
+	virtual int CalcECOperator( DebugFlags debugFlags = None );
+
 	//! The material data storage in the sub-grid area's will not be filled by the base-operator. Check and do this here!
 	void FillMissingDataStorage();
 

FDTD/operator_multithread.h

@@ -38,8 +38,6 @@ public:
 	static Operator_Multithread* New(unsigned int numThreads = 0);
 	virtual ~Operator_Multithread();
 
-	virtual int CalcECOperator( DebugFlags debugFlags = None );
-
 	virtual void setNumThreads( unsigned int numThreads );
 
 	virtual Engine* CreateEngine() const;
@@ -55,6 +53,8 @@ protected:
 	unsigned int (*m_Nr_PEC_thread)[3]; //count PEC edges per thread
 	virtual bool CalcPEC(); //this method is using multi-threading
 
+	virtual int CalcECOperator( DebugFlags debugFlags = None );
+
 	//Calc_EC barrier
 	boost::barrier* m_CalcEC_Start;
 	boost::barrier* m_CalcEC_Stop;

FDTD/operator_sse_compressed.h

@@ -45,8 +45,6 @@ public:
 
 	virtual Engine* CreateEngine() const;
 
-	virtual int CalcECOperator( DebugFlags debugFlags = None );
-
 	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVV(n,x,y,z);}
 	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVI(n,x,y,z);}
 	inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetII(n,x,y,z);}
@@ -73,6 +71,8 @@ protected:
 
 	bool CompareOperators(unsigned int pos1[3], unsigned int pos2[3]);
 
+	virtual int CalcECOperator( DebugFlags debugFlags = None );
+
 	// engine needs access
 public:
 	unsigned int*** m_Op_index;