CylinderCoords: include r=0 and many fixes & necessary changes

- r=0 case included... needs lots of testing...
 - field processing can't access mesh directly --> use operator methods

FDTD/engine_cylinder.cpp

@@ -46,7 +46,28 @@ inline void Engine_Cylinder::CloseAlphaVoltages()
 			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]];
+		}
 	}
 }
 
@@ -74,6 +95,10 @@ bool Engine_Cylinder::IterateTS(unsigned int iterTS)
 	for (unsigned int iter=0;iter<iterTS;++iter)
 	{
 		UpdateVoltages();
+
+		if (cyl_Op->GetR0Included())
+			R0IncludeVoltages();
+
 		ApplyVoltageExcite();
 
 		CloseAlphaVoltages();

FDTD/engine_cylinder.h

@@ -35,6 +35,8 @@ protected:
 	virtual inline void CloseAlphaVoltages();
 	virtual inline void CloseAlphaCurrents();
 
+	virtual inline void R0IncludeVoltages();
+
 	const Operator_Cylinder* cyl_Op;
 };
 

FDTD/operator.cpp

@@ -102,6 +102,42 @@ unsigned int Operator::CalcNyquistNum(double fmax)
 	return floor(T0/2/dT);
 }
 
+double Operator::GetMeshDelta(int n, unsigned int* pos, bool dualMesh) const
+{
+	if ((n<0) || (n>2)) return 0.0;
+	int i_pos[] = {pos[0],pos[1],pos[2]};
+	return GetMeshDelta(n,i_pos,dualMesh);
+}
+
+double Operator::GetMeshDelta(int n, int* pos, bool dualMesh) const
+{
+	if ((n<0) || (n>2)) return 0.0;
+	if (dualMesh==false)
+		return fabs(MainOp->GetIndexDelta(n,pos[n]))*gridDelta;
+	else
+		return fabs(MainOp->GetIndexWidth(n,pos[n]))*gridDelta;
+}
+
+double Operator::GetDiscLine(int n, unsigned int pos, bool dualMesh) const
+{
+	return GetDiscLine(n,(int)pos,dualMesh);
+}
+
+double Operator::GetDiscLine(int n, int pos, bool dualMesh) const
+{
+	if ((n<0) || (n>2)) return 0.0;
+	if ((pos<0) || (pos>=(int)numLines[n])) return 0.0;
+	if (dualMesh==false)
+		return discLines[n][pos];
+	else
+	{
+		if (pos<(int)numLines[n]-1)
+			return 0.5*(discLines[n][pos+1]+discLines[n][pos]);
+		else
+			return 0.5*(discLines[n][pos]+discLines[n][pos-1]);
+	}
+}
+
 bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, bool* inside)
 {
 	bool ok=true;
@@ -451,6 +487,7 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo)
 	for (int n=0;n<3;++n)
 	{
 		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]);
@@ -821,11 +858,16 @@ double Operator::CalcTimestep()
 					ipos_PPM= MainOp->Shift(nPP,-1);
 					MainOp->ResetShift();
 					newT = 2/sqrt( ( 4/EC_L[nP][ipos] + 4/EC_L[nP][ipos_PPM] + 4/EC_L[nPP][ipos] + 4/EC_L[nPP][ipos_PM]) / EC_C[n][ipos] );
-					if (newT<dT) dT=newT;
+					if ((newT<dT) && (newT>0.0)) dT=newT;
 				}
 			}
 		}
 	}
+	if (dT==0)
+	{
+		cerr << "Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit!" << endl;
+		exit(3);
+	}
 //	cerr << "Operator Timestep: " << dT << endl;
 	return 0;
 }
@@ -852,13 +894,11 @@ bool Operator::CalcEFieldExcitation()
 			for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 			{
 				delta[0]=fabs(MainOp->GetIndexDelta(0,pos[0]));
+				for (int n=0;n<3;++n)
+				{
 					coord[0] = discLines[0][pos[0]];
 					coord[1] = discLines[1][pos[1]];
 					coord[2] = discLines[2][pos[2]];
-//				CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::ELECTRODE | CSProperties::METAL));
-				CSProperties* prop = NULL;
-				for (int n=0;n<3;++n)
-				{
 					coord[n]+=delta[n]*0.5;
 					CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::ELECTRODE));
 					if (prop)
@@ -866,9 +906,9 @@ bool Operator::CalcEFieldExcitation()
 						CSPropElectrode* elec = prop->ToElectrode();
 						if (elec!=NULL)
 						{
-							if ((elec->GetActiveDir(n)) && (pos[n]<numLines[n]-1))
+							if ((elec->GetActiveDir(n)) )//&& (pos[n]<numLines[n]-1))
 							{
-								amp = elec->GetWeightedExcitation(n,coord)*delta[n]*gridDelta;
+								amp = elec->GetWeightedExcitation(n,coord)*GetMeshDelta(n,pos);// delta[n]*gridDelta;
 								if (amp!=0)
 								{
 									vExcit.push_back(amp);
@@ -886,7 +926,6 @@ bool Operator::CalcEFieldExcitation()
 							}
 						}
 					}
-					coord[n]-=delta[n]*0.5;
 				}
 			}
 		}

FDTD/operator.h

@@ -59,7 +59,7 @@ public:
 	virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];}
 
 	void SetNyquistNum(unsigned int nyquist) {m_nyquistTS=nyquist;}
-	unsigned int GetNyquistNum() const {return m_nyquistTS;};
+	unsigned int GetNyquistNum() const {return m_nyquistTS;}
 	unsigned int CalcNyquistNum(double fmax);
 
 	void ShowStat() const;
@@ -67,15 +67,22 @@ public:
 	void DumpOperator2File(string filename);
 	void DumpMaterial2File(string filename);
 
-	virtual void Reset();
+	virtual double GetGridDelta() const {return gridDelta;}
+	//! Get the mesh delta times the grid delta for a 3D position
+	virtual double GetMeshDelta(int n, int* pos, bool dualMesh=false) const;
+	virtual double GetMeshDelta(int n, unsigned int* pos, bool dualMesh=false) const;
 
-	bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL);
+	//! Get the disc line in n direction
+	virtual double GetDiscLine(int n, int pos, bool dualMesh=false) const;
+	virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const;
+	virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL);
 
 protected:
 	//! use New() for creating a new Operator
 	Operator();
 
 	virtual void Init();
+	virtual void Reset();
 	virtual void InitOperator();
 
 	struct Grid_Path
@@ -111,14 +118,13 @@ protected:
 	double* EC_R[3];
 
 	unsigned int numLines[3];
-
-	// engine/post-proc needs access
-public:
 	double* discLines[3];
 	double gridDelta;
 	AdrOp* MainOp;
 	AdrOp* DualOp;
 
+	// engine/post-proc needs access
+public:
 	//EC operator
 	FDTD_FLOAT**** vv; //calc new voltage from old voltage
 	FDTD_FLOAT**** vi; //calc new voltage from old current

FDTD/operator_cylinder.cpp

@@ -37,12 +37,26 @@ 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();
 }
 
 inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const
@@ -54,13 +68,23 @@ inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const
 	return numLines[ny];
 }
 
+double Operator_Cylinder::GetMeshDelta(int n, int* pos, bool dualMesh) const
+{
+	double delta = Operator::GetMeshDelta(n,pos,dualMesh);
+	if (delta==0) return delta;
+	if (n==1)
+	{
+		return delta * GetDiscLine(n,pos[0],dualMesh);
+	}
+	return delta;
+}
+
 
 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
 	{
 		cout << "Operator_Cylinder::SetGeometryCSX: Alpha is a full 2*PI => closed Cylinder..." << endl;
@@ -87,13 +111,41 @@ bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
 	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;
 }
 
+int Operator_Cylinder::CalcECOperator()
+{
+	int val = Operator::CalcECOperator();
+	if (val)
+		return val;
+
+	if (CC_R0_included==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]-CC_closedAlpha;++pos[1])
+		{
+			Calc_ECPos(2,pos,inEC);
+			C+=inEC[0];
+			G+=inEC[1];
+		}
+		cerr << C << " and " << G << endl;
+		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;
+}
+
 inline void Operator_Cylinder::Calc_ECOperatorPos(int n, unsigned int* pos)
 {
 	unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]);
@@ -118,6 +170,12 @@ inline void Operator_Cylinder::Calc_ECOperatorPos(int n, unsigned int* pos)
 		ii[n][pos[0]][pos[1]][pos[2]] = 0;
 		iv[n][pos[0]][pos[1]][pos[2]] = 0;
 	}
+
+	if (CC_R0_included && (n==2) && (pos[0]==0))
+	{
+		vv[n][pos[0]][pos[1]][pos[2]] = 1;
+		vi[n][pos[0]][pos[1]][pos[2]] = 0;
+	}
 }
 
 void Operator_Cylinder::ApplyElectricBC(bool* dirs)
@@ -165,7 +223,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	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;
+	double geom_factor=0,A_n=0;
 
 	//******************************* epsilon,kappa averaging *****************************//
 	//shift up-right
@@ -176,13 +234,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	switch (n)
 	{
 	case 0:
-		geom_factor = fabs(deltaPP*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		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;
+		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;
@@ -206,13 +264,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	switch (n)
 	{
 	case 0:
-		geom_factor = fabs(deltaPP*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		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;
+		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;
@@ -236,13 +294,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	switch (n)
 	{
 	case 0:
-		geom_factor = fabs(deltaPP_M*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		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;
+		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;
@@ -266,13 +324,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	switch (n)
 	{
 	case 0:
-		geom_factor = fabs(deltaPP_M*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25;
+		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;
+		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;
@@ -288,15 +346,12 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 		inEC[1] += 0;
 	}
 
-	if (CC_R0_included && (n>0) && (coord[0]==0))
+	if (CC_R0_included && (n==1) && (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;
@@ -306,7 +361,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	double delta_n = fabs(delta_M);
 	if (n==1)
 	{
-		delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius
+		delta_n = delta_n * fabs(coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius
 	}
 	if (prop)
 	{
@@ -330,7 +385,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	delta_n = fabs(delta);
 	if (n==1)
 	{
-		delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius
+		delta_n = delta_n * fabs(coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius
 	}
 	if (prop)
 	{
@@ -360,7 +415,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	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))
+//	if ((n==1) && (pos[1]==0) && (pos[2]==0))
 //		cerr <<  inEC[2]/(coord[0]) <<  endl;
 //	cerr << n << " -> " <<  pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[2] << endl;
 
@@ -369,28 +424,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 
 bool Operator_Cylinder::Calc_EffMatPos(int n, unsigned int* pos, double* inMat)
 {
-	return false; //fixme
+	cerr << "Operator_Cylinder::Calc_EffMatPos: Warning! method not implemented yet..." << endl;
-
+	return false;
-//	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

@@ -28,19 +28,24 @@ public:
 	
 	virtual bool SetGeometryCSX(ContinuousStructure* geo);
 
+	virtual int CalcECOperator();
+
 	virtual void ApplyElectricBC(bool* dirs);
 	virtual void ApplyMagneticBC(bool* dirs);
 
-	virtual void Reset();
-
 	virtual unsigned int GetNumberOfLines(int ny) const;
 
+	//! Get the mesh delta times the grid delta for a 3D position, including radius corrected alpha-mesh width
+	virtual double GetMeshDelta(int n, int* pos, bool dualMesh=false) const;
+
 	bool GetClosedAlpha() const {return CC_closedAlpha;}
 	bool GetR0Included() const {return CC_R0_included;}
 
 protected:
 	Operator_Cylinder();
 	virtual void Init();
+	virtual void InitOperator();
+	virtual void Reset();
 
 	virtual inline void Calc_ECOperatorPos(int n, unsigned int* pos);
 
@@ -49,6 +54,12 @@ protected:
 
 	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/processfields.cpp

@@ -128,7 +128,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
 			lines.clear();
 			for (unsigned int i=start[n];i<=stop[n];i+=subSample[n])
 			{
-				lines.push_back(Op->discLines[n][i]);
+				lines.push_back(Op->GetDiscLine(n,i));//Op->discLines[n][i]);
 			}
 			numLines[n] = lines.size();
 			delete[] discLines[n];
@@ -156,7 +156,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
 			lines.clear();
 			for (unsigned int i=start[n];i<stop[n];i+=subSample[n])
 			{
-				lines.push_back(0.5*(Op->discLines[n][i+1] +  Op->discLines[n][i]));
+				lines.push_back(Op->GetDiscLine(n,i,true));//0.5*(Op->discLines[n][i+1] +  Op->discLines[n][i]));
 			}		
 			numDLines[n] = lines.size();
 			delete[] discDLines[n];
@@ -280,7 +280,7 @@ void ProcessFields::WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT*
 		{
 			for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 			{
-				file << array[pos[0]][pos[1]][pos[2]] << " ";
+				file << setprecision(precision) << array[pos[0]][pos[1]][pos[2]] << " ";
 				++count;
 				if (count%10==0)
 					file << endl;

FDTD/processfields_td.cpp

@@ -52,17 +52,26 @@ void ProcessFieldsTD::DumpCellInterpol(string filename)
 				{
 					OpPos[2]=start[2]+pos[2]*subSample[0];
 					//in x
-					delta  = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]];
+					delta  = Op->GetMeshDelta(0,OpPos); //Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]];
+					if (delta)
+					{
 						E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]+1];
-					E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta);
+						E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta);
+					}
 					//in y
-					delta  = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]];
+					delta  = Op->GetMeshDelta(1,OpPos); //Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]];
+					if (delta)
+					{
 						E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]+1];
-					E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta);
+						E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta);
+					}
 					//in z
-					delta  = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]];
+					delta  = Op->GetMeshDelta(2,OpPos); //Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]];
+					if (delta)
+					{
 						E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]+1][OpPos[2]];
-					E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta);
+						E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta);
+					}
 				}
 			}
 		}
@@ -104,18 +113,30 @@ void ProcessFieldsTD::DumpCellInterpol(string filename)
 				{
 					OpPos[2]=start[2]+pos[2]*subSample[2];
 					//in x
-					if (OpPos[0]==0) delta  = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]];
+//					if (OpPos[0]==0) delta  = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]];
-					else delta = 0.5* (Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]-1]);
+//					else delta = 0.5* (Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]-1]);
+					delta = Op->GetDiscLine(0,OpPos[0],true);
+					if (delta)
+					{
 						H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + curr[0][OpPos[0]+1][OpPos[1]][OpPos[2]];
-					H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta);
+						H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta);
+					}
 					//in y
-					delta  = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]];
+//					delta  = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]];
+					delta = Op->GetDiscLine(1,OpPos[1],true);
+					if (delta)
+					{
 						H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]] + curr[1][OpPos[0]][OpPos[1]+1][OpPos[2]];
-					H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta);
+						H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta);
+					}
 					//in z
-					delta  = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]];
+//					delta  = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]];
+					delta = Op->GetDiscLine(2,OpPos[2],true);
+					if (delta)
+					{
 						H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]] + curr[2][OpPos[0]][OpPos[1]][OpPos[2]+1];
-					H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta);
+						H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta);
+					}
 				}
 			}
 		}
@@ -154,18 +175,21 @@ void ProcessFieldsTD::DumpNoInterpol(string filename)
 		for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 		{
 			OpPos[0]=start[0]+pos[0]*subSample[0];
-			delta[0]=fabs(Op->MainOp->GetIndexDelta(0,OpPos[0]));
+			delta[0]=Op->GetMeshDelta(0,OpPos);//fabs(Op->MainOp->GetIndexDelta(0,OpPos[0]));
 			for (pos[1]=0;pos[1]<numLines[1];++pos[1])
 			{
 				OpPos[1]=start[1]+pos[1]*subSample[1];
-				delta[1]=fabs(Op->MainOp->GetIndexDelta(1,OpPos[1]));
+				delta[1]=Op->GetMeshDelta(1,OpPos);//fabs(Op->MainOp->GetIndexDelta(1,OpPos[1]));
 				for (pos[2]=0;pos[2]<numLines[2];++pos[2])
 				{
 					OpPos[2]=start[2]+pos[2]*subSample[2];
-					delta[2]=fabs(Op->MainOp->GetIndexDelta(2,OpPos[2]));
+					delta[2]=Op->GetMeshDelta(2,OpPos);//fabs(Op->MainOp->GetIndexDelta(2,OpPos[2]));
-					E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0]/Op->gridDelta;
+					if (delta[0])
-					E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1]/Op->gridDelta;
+						E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0];// /Op->gridDelta;
-					E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2]/Op->gridDelta;
+					if (delta[1])
+						E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1];// /Op->gridDelta;
+					if (delta[2])
+						E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2];// /Op->gridDelta;
 				}
 			}
 		}
@@ -196,19 +220,22 @@ void ProcessFieldsTD::DumpNoInterpol(string filename)
 		for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 		{
 			OpPos[0]=start[0]+pos[0]*subSample[0];
-			delta[0]=fabs(Op->MainOp->GetIndexWidth(0,OpPos[0]));
+			delta[0]=Op->GetMeshDelta(0,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(0,OpPos[0]));
 			for (pos[1]=0;pos[1]<numLines[1];++pos[1])
 			{
 				OpPos[1]=start[1]+pos[1]*subSample[1];
-				delta[1]=fabs(Op->MainOp->GetIndexWidth(1,OpPos[1]));
+				delta[1]=Op->GetMeshDelta(1,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(1,OpPos[1]));
 				for (pos[2]=0;pos[2]<numLines[2];++pos[2])
 				{
 					OpPos[2]=start[2]+pos[2]*subSample[2];
-					delta[2]=fabs(Op->MainOp->GetIndexWidth(2,OpPos[2]));
+					delta[2]=Op->GetMeshDelta(2,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(2,OpPos[2]));
 					//in x
-					H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0]/Op->gridDelta;
+					if (delta[0])
-					H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1]/Op->gridDelta;
+						H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0];// /Op->gridDelta;
-					H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2]/Op->gridDelta;
+					if (delta[1])
+						H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1];// /Op->gridDelta;
+					if (delta[2])
+						H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2];// /Op->gridDelta;
 				}
 			}
 		}

openems.cpp

@@ -237,6 +237,7 @@ int openEMS::SetupFDTD(const char* file)
 	if (CylinderCoords)
 	{
 		FDTD_Op = Operator_Cylinder::New();
+		CSX.SetCoordInputType(1); //tell CSX to use cylinder-coords
 	}
 	else
 	{