major operator revision

- now the equivalent circuits are calculated by first calculating the averaged material properties
- this approach should also be save for the cylindrical FDTD

This needs some further testing, especially for the cylindrical operator!!

FDTD/operator.cpp

@@ -145,7 +145,7 @@ double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
 {
 	int nyP = (ny+1)%3;
 	int nyPP = (ny+2)%3;
-	return GetMeshDelta(nyP,pos,!dualMesh) * GetMeshDelta(nyPP,pos,!dualMesh);
+	return GetNodeWidth(nyP,pos,dualMesh) * GetNodeWidth(nyPP,pos,dualMesh);
 }
 
 bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, bool* inside)
@@ -646,9 +646,51 @@ void Operator::ApplyMagneticBC(bool* dirs)
 	}
 }
 
-
+bool Operator::Calc_ECPos(int ny, const unsigned int* pos, double* EC) const
-bool Operator::Calc_ECPos(int n, const unsigned int* pos, double* inEC) const
 {
+	double EffMat[4];
+	Calc_EffMatPos(ny,pos,EffMat);
+
+
+	double delta = GetEdgeLength(ny,pos);
+	double area  = GetEdgeArea(ny,pos);
+
+//	if (isnan(EffMat[0]))
+//	{
+//		cerr << ny << " " << pos[0] << " " << pos[1] << " " << pos[2] << " : " << EffMat[0] << endl;
+//	}
+
+	if (delta)
+	{
+		EC[0] = EffMat[0] * area/delta;
+		EC[1] = EffMat[1] * area/delta;
+	}
+	else
+	{
+		EC[0] = 0;
+		EC[1] = 0;
+	}
+
+	delta = GetEdgeLength(ny,pos,true);
+	area  = GetEdgeArea(ny,pos,true);
+
+	if (delta)
+	{
+		EC[2] = EffMat[2] * area/delta;
+		EC[3] = EffMat[3] * area/delta;
+	}
+	else
+	{
+		EC[2] = 0;
+		EC[3] = 0;
+	}
+
+	return true;
+}
+
+bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const
+{
+	int n=ny;
 	double coord[3];
 	double shiftCoord[3];
 	int nP = (n+1)%3;
@@ -663,132 +705,158 @@ bool Operator::Calc_ECPos(int n, const unsigned int* pos, double* inEC) const
 	double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
 	double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
 
+	int loc_pos[3]={pos[0],pos[1],pos[2]};
+	double A_n;
+	double area = 0;
+
 	//******************************* 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;
+	A_n = GetNodeArea(ny,loc_pos,true);
+//	{
+//		cerr << ny << " " << pos[0] << " " <<  pos[1] << " " <<  pos[2] << ": " << A_n << endl;
+//		exit(0);
+//	}
 	CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] = mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
+		EffMat[0] = mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		inEC[1] = mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP);
+		EffMat[1] = mat->GetKappaWeighted(n,shiftCoord)*A_n;
 	}
 	else
 	{
-		inEC[0] = 1*fabs(deltaP*deltaPP);
+		EffMat[0] = 1*A_n;
-		inEC[1] = 0;
+		EffMat[1] = 0;
 	}
+	area+=A_n;
+
 	//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;
+
+	--loc_pos[nP];
+	A_n = GetNodeArea(ny,loc_pos,true);
+//	cerr << A_n << endl;
 	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP_M*deltaPP);
+		EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP_M*deltaPP);
+		EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n;
 	}
 	else
 	{
-		inEC[0] += 1*fabs(deltaP_M*deltaPP);
+		EffMat[0] += 1*A_n;
-		inEC[1] += 0;
+		EffMat[1] += 0;
 	}
+	area+=A_n;
 
 	//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;
+	++loc_pos[nP];
+	--loc_pos[nPP];
+	A_n = GetNodeArea(ny,loc_pos,true);
 	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP_M);
+		EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP_M);
+		EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n;
 	}
 	else
 	{
-		inEC[0] += 1*fabs(deltaP*deltaPP_M);
+		EffMat[0] += 1*A_n;
-		inEC[1] += 0;
+		EffMat[1] += 0;
 	}
+	area+=A_n;
 
 	//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;
+	--loc_pos[nP];
+	A_n = GetNodeArea(ny,loc_pos,true);
 	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP_M*deltaPP_M);
+		EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n;
-		inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP_M*deltaPP_M);
+		EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n;
 	}
 	else
 	{
-		inEC[0] += 1*fabs(deltaP_M*deltaPP_M);
+		EffMat[0] += 1*A_n;
-		inEC[1] += 0;
+		EffMat[1] += 0;
 	}
+	area+=A_n;
 
-	inEC[0]*=gridDelta/fabs(delta)/4.0*__EPS0__;
+	EffMat[0]*=__EPS0__/area;
-	inEC[1]*=gridDelta/fabs(delta)/4.0;
+	EffMat[1]/=area;
 
 	//******************************* mu,sigma averaging *****************************//
+	loc_pos[0]=pos[0];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,loc_pos,true);
 	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[2] = fabs(delta_M) / mat->GetMueWeighted(n,shiftCoord);
+		EffMat[2] = delta_ny / mat->GetMueWeighted(n,shiftCoord);
-		if (mat->GetSigma(n))
+		if (mat->GetSigmaWeighted(n,shiftCoord))
-			inEC[3] = fabs(delta_M) / mat->GetSigmaWeighted(n,shiftCoord);
+			EffMat[3] = delta_ny / mat->GetSigmaWeighted(n,shiftCoord);
 		else
-			inEC[3] = 0;
+			EffMat[3] = 0;
 	}
 	else
 	{
-		inEC[2] = fabs(delta_M);
+		EffMat[2] = delta_ny;
-		inEC[3] = 0;
+		EffMat[3] = 0;
 	}
+	length=delta_ny;
+
 	//shift up
 	shiftCoord[n] = coord[n]+delta*0.25;
 	shiftCoord[nP] = coord[nP]+deltaP*0.5;
 	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
+	++loc_pos[n];
+	delta_ny = GetNodeWidth(n,loc_pos,true);
 	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
 	if (prop)
 	{
 		CSPropMaterial* mat = prop->ToMaterial();
-		inEC[2] += fabs(delta)/mat->GetMueWeighted(n,shiftCoord);
+		EffMat[2] += delta_ny / mat->GetMueWeighted(n,shiftCoord);
 		if (mat->GetSigmaWeighted(n,shiftCoord))
-			inEC[3] += fabs(delta)/mat->GetSigmaWeighted(n,shiftCoord);
+			EffMat[3] += delta_ny/mat->GetSigmaWeighted(n,shiftCoord);
 		else
-			inEC[3] = 0;
+			EffMat[3] = 0;
 	}
 	else
 	{
-		inEC[2] += fabs(delta);
+		EffMat[2] += 1*delta_ny;
-		inEC[3] = 0;
+		EffMat[3] = 0;
 	}
+	length+=delta_ny;
 
-	inEC[2] = gridDelta * fabs(deltaP*deltaPP) * 2.0 * __MUE0__ / inEC[2];
+	EffMat[2] = length * __MUE0__ / EffMat[2];
-	if (inEC[3]) inEC[3]=gridDelta*fabs(deltaP*deltaPP) * 2.0 / inEC[3];
+	if (EffMat[3]) EffMat[3]=length / EffMat[3];
 
-	return true;
+	for (int n=0;n<4;++n)
-}
+		if (isnan(EffMat[n]) || isinf(EffMat[n]))
-
-bool Operator::Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const
 		{
-	this->Calc_ECPos(n,pos,inMat);
+			cerr << "Operator::Calc_EffMatPos: An effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
-
+			exit(0);
-	inMat[0] *= GetMeshDelta(n,pos)/GetNodeArea(n,pos);
+		}
-	inMat[1] *= GetMeshDelta(n,pos)/GetNodeArea(n,pos);
-
-	inMat[2] *= GetMeshDelta(n,pos,true)/GetNodeArea(n,pos,true);
-	inMat[3] *= GetMeshDelta(n,pos,true)/GetNodeArea(n,pos,true);
 
 	return true;
 }

FDTD/operator.h

@@ -44,6 +44,12 @@ public:
 
 	virtual int CalcECOperator();
 
+	//! 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;
+
+	//! 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;
+
 	inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; }
 	inline virtual FDTD_FLOAT& GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; }
 
@@ -83,11 +89,30 @@ 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 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 {return GetNodeWidth(ny,(const int*)pos,dualMesh);}
+	//! 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 {return GetMeshDelta(ny,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 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;
 
+	//! Get the length of an FDTD edge.
+	virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeLength(ny,(const int*)pos,dualMesh);}
+	//! Get the length of an FDTD edge.
+	virtual double GetEdgeLength(int ny, const int pos[3], bool dualMesh = false) const {return GetMeshDelta(ny,pos,dualMesh);}
+
+	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
+	/*!
+		This will return the area around an edge with a given direction, measured at the middle of the edge.
+		In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems.
+	*/
+	virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeArea(ny,(const int*)pos,dualMesh);}
+	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos \sa GetEdgeArea
+	virtual double GetEdgeArea(int ny, const int pos[3], bool dualMesh = false) const  {return GetNodeArea(ny,(const int*)pos,dualMesh);};
+
 	virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL);
 
 	virtual void AddExtension(Operator_Extension* op_ext);
@@ -136,8 +161,6 @@ protected:
 	//EC elements, internal only!
 	virtual void Init_EC();
 	virtual bool Calc_EC();
-	virtual bool Calc_ECPos(int n, const unsigned int* pos, double* inEC) const;
-	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

@@ -85,6 +85,20 @@ double Operator_Cylinder::GetMeshDelta(int n, const int* pos, bool dualMesh) con
 	return delta;
 }
 
+double Operator_Cylinder::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const
+{
+	if ((ny<0) || (ny>2)) return 0.0;
+	double width = 0;
+	if (dualMesh)
+		width = fabs(MainOp->GetIndexDelta(ny,pos[ny]))*gridDelta;
+	else
+		width = fabs(MainOp->GetIndexWidth(ny,pos[ny]))*gridDelta;
+	if (ny==1)
+		width *= GetDiscLine(0,pos[0],dualMesh);
+	return width;
+}
+
+
 double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
 {
 	if (ny==2)
@@ -105,9 +119,33 @@ double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) c
 			return da * pow(r2,2);
 		return da/2* (pow(r2,2) - pow(r1,2));
 	}
+
+	if (ny==0)
+	{
+		if (dualMesh)
+			return fabs(MainOp->GetIndexDelta(1,pos[1]) * MainOp->GetIndexDelta(2,pos[2]) * GetDiscLine(0,pos[0],true) * gridDelta * gridDelta);
+		else
+			return fabs(MainOp->GetIndexDelta(1,pos[1]) * MainOp->GetIndexDelta(2,pos[2]) * GetDiscLine(0,pos[0],false) * gridDelta * gridDelta);
+	}
 	return __OP_CYLINDER_BASE_CLASS__::GetNodeArea(ny,pos,dualMesh);
 }
 
+double Operator_Cylinder::GetEdgeLength(int ny, const int pos[3], bool dualMesh) const
+{
+	double length = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(ny,pos,dualMesh);
+	if (ny!=1)
+		return length;
+	return length * GetDiscLine(0,pos[0],dualMesh);
+}
+
+double Operator_Cylinder::GetEdgeArea(int ny, const int pos[3], bool dualMesh) const
+{
+	if (ny!=0)
+		return GetNodeArea(ny,pos,dualMesh);
+
+	return GetEdgeLength(1,pos,!dualMesh) * GetEdgeLength(2,pos,!dualMesh);
+}
+
 bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
 {
 	if (__OP_CYLINDER_BASE_CLASS__::SetGeometryCSX(geo)==false) return false;
@@ -188,235 +226,6 @@ void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
 	__OP_CYLINDER_BASE_CLASS__::ApplyMagneticBC(dirs);
 }
 
-bool Operator_Cylinder::Calc_ECPos(int n, const unsigned int* pos, double* inEC) const
-{
-	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]];
-	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=0,A_n=0;
-
-	//******************************* 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;
-	}
-
-	//******************************* 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 * fabs(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 * fabs(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->GetSigmaWeighted(n,shiftCoord))
-			inEC[3] += delta_n/mat->GetSigmaWeighted(n,shiftCoord);
-		else
-			inEC[3] = 0;
-	}
-	else
-	{
-		inEC[2] += delta_n;
-		inEC[3] = 0;
-	}
-
-	A_n = fabs(deltaP*deltaPP);
-	if (n==0) //x-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==1) && (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, const unsigned int* pos, double* inMat) const
-{
-	__OP_CYLINDER_BASE_CLASS__::Calc_EffMatPos(n, pos, inMat);
-
-	// 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)
 {
 	if (op_ext->IsCylinderCoordsSave())

FDTD/operator_cylinder.h

@@ -36,6 +36,13 @@ public:
 
 	virtual bool SetGeometryCSX(ContinuousStructure* geo);
 
+
+//	virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC) 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 void ApplyElectricBC(bool* dirs);
 	virtual void ApplyMagneticBC(bool* dirs);
 
@@ -47,19 +54,29 @@ 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 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 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;
 
+	//! Get the length of an FDTD edge.
+	virtual double GetEdgeLength(int ny, const int pos[3], bool dualMesh = false) const;
+
+	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
+	/*!
+		This will return the area around an edge with a given direction, measured at the middle of the edge.
+		In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems.
+	*/
+	virtual double GetEdgeArea(int ny, const int pos[3], 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, const unsigned int* pos, double* inEC) const;
-	virtual bool Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const;
-
 protected:
 	Operator_Cylinder();
 	virtual void Init();

FDTD/operator_ext_pml_sf.cpp

@@ -251,22 +251,22 @@ void Operator_Ext_PML_SF_Plane::SetPMLLength(int width)
 
 }
 
-double Operator_Ext_PML_SF_Plane::GetNodeArea(int ny, unsigned int pos[3], bool dualMesh) const
+double Operator_Ext_PML_SF_Plane::GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh) const
 {
 	unsigned int l_pos[] = {pos[0],pos[1],pos[2]};
 	l_pos[m_ny] = m_LineNr;
 
-	return m_Op->GetNodeArea(ny,l_pos,dualMesh);
+	return m_Op->GetEdgeArea(ny,l_pos,dualMesh);
 }
 
-double Operator_Ext_PML_SF_Plane::GetNodeLength(int ny, unsigned int pos[3], bool dualMesh) const
+double Operator_Ext_PML_SF_Plane::GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh) const
 {
 	if (ny==m_ny)
 		return m_pml_delta;
 
 	unsigned int l_pos[] = {pos[0],pos[1],pos[2]};
 	l_pos[m_ny] = m_LineNr;
-	return m_Op->GetMeshDelta(ny,l_pos,dualMesh);
+	return m_Op->GetEdgeLength(ny,l_pos,dualMesh);
 }
 
 double Operator_Ext_PML_SF_Plane::GetKappaGraded(double depth, double Zm) const
@@ -311,13 +311,13 @@ bool Operator_Ext_PML_SF_Plane::Calc_ECPos(int nP, int n, unsigned int* pos, dou
 		}
 	}
 
-	double geomFactor = GetNodeArea(n,pos) / GetNodeLength(n,pos);
+	double geomFactor = GetEdgeArea(n,pos) / GetEdgeLength(n,pos);
 	if (geomFactor<=0 || isnan(geomFactor) || isinf(geomFactor)) //check if geomFactor is positive, not zero and a valid number (necessary for cylindrical coords)
 		geomFactor = 0;
 	inEC[0] = inMat[0] * geomFactor;
 	inEC[1] = (inMat[1]+kappa) * geomFactor;
 
-	geomFactor = GetNodeArea(n,pos) / GetNodeLength(n,pos);
+	geomFactor = GetEdgeArea(n,pos,true) / GetEdgeLength(n,pos,true);
 	if (geomFactor<=0 || isnan(geomFactor) || isinf(geomFactor)) //check if geomFactor is positive, not zero and a valid number (necessary for cylindrical coords)
 		geomFactor = 0;
 	inEC[2] = inMat[2] * geomFactor;

FDTD/operator_ext_pml_sf.h

@@ -41,8 +41,8 @@ public:
 	inline virtual FDTD_FLOAT& GetII(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[nP][n][x][y][z]; }
 	inline virtual FDTD_FLOAT& GetIV(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[nP][n][x][y][z]; }
 
-	virtual double GetNodeArea(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;}
+	virtual double GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;}
-	virtual double GetNodeLength(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;}
+	virtual double GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;}
 
 	//! This will resturn the pml parameter grading
 	virtual double GetKappaGraded(double depth, double Zm) const {UNUSED(depth);UNUSED(Zm);return 0.0;}
@@ -105,8 +105,8 @@ public:
 	void SetDirection(int ny, bool top_ny);
 	void SetPMLLength(int width);
 
-	virtual double GetNodeArea(int ny, unsigned int pos[3], bool dualMesh = false) const;
+	virtual double GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh = false) const;
-	virtual double GetNodeLength(int ny, unsigned int pos[3], bool dualMesh = false) const;
+	virtual double GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh = false) const;
 
 	virtual double GetKappaGraded(double depth, double Zm) const;