Operator: GetNodeArea & Update/fix in Calc_EffMatPos method

- new method GetNodeArea will return the area of a node for a given direction
- methods Calc_ECPos & Calc_EffMatPos now const
- Calc_EffMatPos in Operator class updated to use new functions GetMeshDelta & GetNodeArea
- Calc_EffMatPos introduced (fixed) in Operator_Cylinder (need some testing)
- treatment of E_alpha=0 at r==0 moved from Calc_ECPos to Electric-BC

FDTD/operator.cpp

@@ -138,12 +138,14 @@ double Operator::GetDiscLine(int n, int pos, bool dualMesh) const
 	if (dualMesh==false)
 		return discLines[n][pos];
 	else
-	{
+		return (discLines[n][pos] + 0.5*fabs(MainOp->GetIndexDelta(n,pos)));
-		if (pos<(int)numLines[n]-1)
+}
-			return 0.5*(discLines[n][pos+1]+discLines[n][pos]);
+
-		else
+double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
-			return 0.5*(discLines[n][pos]+discLines[n][pos-1]);
+{
-	}
+	int nyP = (ny+1)%3;
+	int nyPP = (ny+2)%3;
+	return GetMeshDelta(nyP,pos,!dualMesh) * GetMeshDelta(nyPP,pos,!dualMesh);
 }
 
 bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, bool* inside)
@@ -645,7 +647,7 @@ void Operator::ApplyMagneticBC(bool* dirs)
 }
 
 
-bool Operator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
+bool Operator::Calc_ECPos(int n, const unsigned int* pos, double* inEC) const
 {
 	double coord[3];
 	double shiftCoord[3];
@@ -778,26 +780,15 @@ bool Operator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	return true;
 }
 
-bool Operator::Calc_EffMatPos(int n, unsigned int* pos, double* inMat)
+bool Operator::Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const
 {
-	int nP = (n+1)%3;
-	int nPP = (n+2)%3;
-
-	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] *= fabs(delta)/(0.25*(fabs(deltaP_M) + fabs(deltaP))*(fabs(deltaPP_M) + fabs(deltaPP)))/gridDelta;
+	inMat[0] *= GetMeshDelta(n,pos)/GetNodeArea(n,pos);
-	inMat[1] *= fabs(delta)/(0.25*(fabs(deltaP_M) + fabs(deltaP))*(fabs(deltaPP_M) + fabs(deltaPP)))/gridDelta;
+	inMat[1] *= GetMeshDelta(n,pos)/GetNodeArea(n,pos);
 
-	inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
+	inMat[2] *= GetMeshDelta(n,pos,true)/GetNodeArea(n,pos,true);
-	inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
+	inMat[3] *= GetMeshDelta(n,pos,true)/GetNodeArea(n,pos,true);
 
 	return true;
 }

FDTD/operator.h

@@ -83,6 +83,11 @@ public:
 	//! Get the disc line in \a n direction (in drawing units)
 	virtual double GetDiscLine(int n, 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 unsigned int pos[3], bool dualMesh = false) const {return GetNodeArea(ny,(const int*)pos,dualMesh);}
+	//! 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;
+
 	virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL);
 
 	virtual void AddExtension(Operator_Extension* op_ext);
@@ -131,8 +136,8 @@ protected:
 	//EC elements, internal only!
 	virtual void Init_EC();
 	virtual bool Calc_EC();
-	virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
+	virtual bool Calc_ECPos(int n, const unsigned int* pos, double* inEC) const;
-	virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
+	virtual bool Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const;
 	double* EC_C[3];
 	double* EC_G[3];
 	double* EC_L[3];

FDTD/operator_cylinder.cpp

@@ -88,6 +88,28 @@ double Operator_Cylinder::GetMeshDelta(int n, const int* pos, bool dualMesh) con
 	return delta;
 }
 
+double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
+{
+	if (ny==2)
+	{
+		double da = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(1,pos,dualMesh)/gridDelta;
+		double r1,r2;
+		if (!dualMesh)
+		{
+			r1 = (discLines[0][pos[0]] - fabs(MainOp->GetIndexDelta(0,pos[0]-1))/2.0)*gridDelta;
+			r2 = (discLines[0][pos[0]] + fabs(MainOp->GetIndexDelta(0,pos[0]  ))/2.0)*gridDelta;
+		}
+		else
+		{
+			r1 =  discLines[0][pos[0]]*gridDelta;
+			r2 = (discLines[0][pos[0]] + fabs(MainOp->GetIndexDelta(0,pos[0])))*gridDelta;
+		}
+		if (r1<0)
+			return da * pow(r2,2);
+		return da/2* (pow(r2,2) - pow(r1,2));
+	}
+	return __OP_CYLINDER_BASE_CLASS__::GetNodeArea(ny,pos,dualMesh);
+}
 
 bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
 {
@@ -138,8 +160,16 @@ void Operator_Cylinder::ApplyElectricBC(bool* dirs)
 	}
 	if (CC_R0_included)
 	{
-		// no special treatment necessary
+		// E in alpha direction ( aka volt[1][x][y][z] ) is not defined for r==0 --> always zero...
-		// operator for z-direction at r=0 will be calculated and set separately
+		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])
+			{
+				GetVV(1,pos[0],pos[1],pos[2]) = 0;
+				GetVI(1,pos[0],pos[1],pos[2]) = 0;
+			}
+		}
 	}
 	__OP_CYLINDER_BASE_CLASS__::ApplyElectricBC(dirs);
 }
@@ -158,7 +188,7 @@ void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
 	__OP_CYLINDER_BASE_CLASS__::ApplyMagneticBC(dirs);
 }
 
-bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
+bool Operator_Cylinder::Calc_ECPos(int n, const unsigned int* pos, double* inEC) const
 {
 	double coord[3];
 	double shiftCoord[3];
@@ -296,12 +326,6 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 		inEC[1] += 0;
 	}
 
-	if (CC_R0_included && (n==1) && (pos[0]==0))
-	{
-		inEC[0]=0;
-		inEC[1]=0;
-	}
-
 	//******************************* mu,sigma averaging *****************************//
 	//shift down
 	shiftCoord[n] = coord[n]-delta_M*0.25;
@@ -372,10 +396,25 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC)
 	return true;
 }
 
-bool Operator_Cylinder::Calc_EffMatPos(int /*n*/, unsigned int* /*pos*/, double* /*inMat*/)
+bool Operator_Cylinder::Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const
 {
-	cerr << "Operator_Cylinder::Calc_EffMatPos: Warning! method not implemented yet..." << endl;
+	__OP_CYLINDER_BASE_CLASS__::Calc_EffMatPos(n, pos, inMat);
-	return false;
+
+	// H_rho is not defined at position r==0
+	if (CC_R0_included && (n==0) && (pos[0]==0))
+	{
+		inMat[2] = 0;
+		inMat[3] = 0;
+	}
+
+	// E_alpha is not defined at position r==0
+	if (CC_R0_included && (n==1) && (pos[0]==0))
+	{
+		inMat[0]=0;
+		inMat[1]=0;
+	}
+
+	return true;
 }
 
 void Operator_Cylinder::AddExtension(Operator_Extension* op_ext)

FDTD/operator_cylinder.h

@@ -47,13 +47,18 @@ public:
 	//! Get the mesh delta times the grid delta for a 3D position, including radius corrected alpha-mesh width
 	virtual double GetMeshDelta(int n, const 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 unsigned int pos[3], bool dualMesh = false) const {return GetNodeArea(ny,(const int*)pos,dualMesh);}
+	//! Get the node area for a given direction \a n and a given mesh posisition \a pos
+	virtual double GetNodeArea(int n, const int* pos, bool dualMesh=false) const;
+
 	bool GetClosedAlpha() const {return CC_closedAlpha;}
 	bool GetR0Included() const {return CC_R0_included;}
 
 	virtual void AddExtension(Operator_Extension* op_ext);
 
-	virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
+	virtual bool Calc_ECPos(int n, const unsigned int* pos, double* inEC) const;
-	virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
+	virtual bool Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const;
 
 protected:
 	Operator_Cylinder();