replaced GetMeshDelta by GetEdgeLength & partially removed AdrOp

Common/operator_base.h

@@ -44,9 +44,6 @@ public:
 	//! Get the grid drawing unit in m
 	virtual double GetGridDelta() const =0;
 
-	//! Get the mesh delta times the grid delta for a 3D position (unit is meter)
-	virtual double GetMeshDelta(int n, const unsigned int* pos, bool dualMesh=false) const =0;
-
 	//! Get the disc line in \a n direction (in drawing units)
 	virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const =0;
 

Common/process_efield.cpp

@@ -91,7 +91,7 @@ int ProcessEField::Process()
 			file << setprecision(m_precision) << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep();
 			for (int n=0; n<3; n++)
 			{
-				FDTD_FLOAT field = Eng->GetVolt(n,start) / Op->GetMeshDelta(n,start);
+				FDTD_FLOAT field = Eng->GetVolt(n,start) / Op->GetEdgeLength(n,start);
 				field *= m_weight;
 //				TD_Values.push_back(voltage);
 				file << "\t" << field;
@@ -108,7 +108,7 @@ int ProcessEField::Process()
 			double T = (double)Eng->GetNumberOfTimesteps() * Op->GetTimestep();
 			for (int pol=0; pol<3; pol++)
 			{
-				FDTD_FLOAT field = Eng->GetVolt(pol,start) / Op->GetMeshDelta(pol,start);
+				FDTD_FLOAT field = Eng->GetVolt(pol,start) / Op->GetEdgeLength(pol,start);
 				field *= m_weight;
 				for (size_t n=0; n<m_FD_Samples.size(); ++n)
 				{

Common/process_hfield.cpp

@@ -75,7 +75,7 @@ int ProcessHField::Process()
 			file << setprecision(m_precision) << ((double)Eng->GetNumberOfTimesteps()+0.5)*Op->GetTimestep();
 			for (int n=0; n<3; n++)
 			{
-				FDTD_FLOAT field = Eng->GetCurr(n,start) / Op->GetMeshDelta(n,start,true);
+				FDTD_FLOAT field = Eng->GetCurr(n,start) / Op->GetEdgeLength(n,start,true);
 				field *= m_weight;
 //				TD_Values.push_back(voltage);
 				file << "\t" << field;
@@ -92,7 +92,7 @@ int ProcessHField::Process()
 			double T = (double)Eng->GetNumberOfTimesteps() * Op->GetTimestep();
 			for (int pol=0; pol<3; pol++)
 			{
-				FDTD_FLOAT field = Eng->GetCurr(pol,start) / Op->GetMeshDelta(pol,start,true);
+				FDTD_FLOAT field = Eng->GetCurr(pol,start) / Op->GetEdgeLength(pol,start,true);
 				field *= m_weight;
 				for (size_t n=0; n<m_FD_Samples.size(); ++n)
 				{

Common/processing.cpp

@@ -209,7 +209,7 @@ void Processing::DumpBox2File( string vtkfilenameprefix, bool dualMesh ) const
 		// line are not displayed correctly -> enlarge
 		for (int i=0; i<3; i++)
 		{
-			double delta = min( Op->GetMeshDelta( i, start,dualMesh ), Op->GetMeshDelta( i, stop,dualMesh ) ) / Op->GetGridDelta() / 4.0;
+			double delta = min( Op->GetEdgeLength( i, start,dualMesh ), Op->GetEdgeLength( i, stop,dualMesh ) ) / Op->GetGridDelta() / 4.0;
 			s1[i] -= delta;
 			s2[i] += delta;
 		}

FDTD/engine_interface_fdtd.cpp

@@ -48,7 +48,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
 	case NODE_INTERPOLATE:
 		for (int n=0; n<3; ++n)
 		{
-			delta = m_Op->GetMeshDelta(n,iPos);
+			delta = m_Op->GetEdgeLength(n,iPos);
 			out[n] = m_Eng->GetVolt(n,iPos);
 			if (delta==0)
 			{
@@ -61,7 +61,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
 				continue;
 			}
 			--iPos[n];
-			double deltaDown = m_Op->GetMeshDelta(n,iPos);
+			double deltaDown = m_Op->GetEdgeLength(n,iPos);
 			double deltaRel = delta / (delta+deltaDown);
 			out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetVolt(n,iPos)/deltaDown*deltaRel;
 			++iPos[n];
@@ -77,19 +77,19 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
 				out[n] = 0; //electric field outside the field domain is always zero
 				continue;
 			}
-			delta = m_Op->GetMeshDelta(n,iPos);
+			delta = m_Op->GetEdgeLength(n,iPos);
 			if (delta)
 				out[n]=m_Eng->GetVolt(n,iPos)/delta;
 			++iPos[nP];
-			delta = m_Op->GetMeshDelta(n,iPos);
+			delta = m_Op->GetEdgeLength(n,iPos);
 			if (delta)
 				out[n]+=m_Eng->GetVolt(n,iPos)/delta;
 			++iPos[nPP];
-			delta = m_Op->GetMeshDelta(n,iPos);
+			delta = m_Op->GetEdgeLength(n,iPos);
 			if (delta)
 				out[n]+=m_Eng->GetVolt(n,iPos)/delta;
 			--iPos[nP];
-			delta = m_Op->GetMeshDelta(n,iPos);
+			delta = m_Op->GetEdgeLength(n,iPos);
 			if (delta)
 				out[n]+=m_Eng->GetVolt(n,iPos)/delta;
 			--iPos[nPP];
@@ -123,13 +123,13 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 				out[n] = 0;
 				continue;
 			}
-			out[n]=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true);
+			out[n]=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
 			--iPos[nP];
-			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true);
+			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
 			--iPos[nPP];
-			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true);
+			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
 			++iPos[nP];
-			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true);
+			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
 			++iPos[nPP];
 			out[n]/=4;
 		}
@@ -137,7 +137,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 	case CELL_INTERPOLATE:
 		for (int n=0; n<3; ++n)
 		{
-			delta = m_Op->GetMeshDelta(n,iPos,true);
+			delta = m_Op->GetEdgeLength(n,iPos,true);
 			out[n] = m_Eng->GetCurr(n,iPos);
 			if ((pos[n]>=m_Op->GetNumberOfLines(n)-1))
 			{
@@ -145,7 +145,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 				continue;
 			}
 			++iPos[n];
-			double deltaUp = m_Op->GetMeshDelta(n,iPos,true);
+			double deltaUp = m_Op->GetEdgeLength(n,iPos,true);
 			double deltaRel = delta / (delta+deltaUp);
 			out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetCurr(n,iPos)/deltaUp*deltaRel;
 			--iPos[n];

FDTD/extensions/operator_ext_cylinder.cpp

@@ -64,8 +64,8 @@ bool Operator_Ext_Cylinder::BuildExtension()
 		for (pos[1]=0; pos[1]<m_Op->GetOriginalNumLines(1)-1; ++pos[1])
 		{
 			m_Op_Cyl->Calc_ECPos(2,pos,inEC);
-			C+=inEC[0]*0.5;
+			C+=inEC[0];
-			G+=inEC[1]*0.5;
+			G+=inEC[1];
 		}
 		m_Op->SetVV(2,0,0,pos[2], 1);
 		vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C);

FDTD/extensions/operator_ext_mur_abc.cpp

@@ -124,7 +124,7 @@ bool Operator_Ext_Mur_ABC::BuildExtension()
 	double dT = m_Op->GetTimestep();
 	unsigned int pos[] = {0,0,0};
 	pos[m_ny] = m_LineNr;
-	double delta = fabs(m_Op->GetMeshDelta(m_ny,pos));
+	double delta = fabs(m_Op->GetEdgeLength(m_ny,pos));
 	double coord[] = {0,0,0};
 	coord[0] = m_Op->GetDiscLine(0,pos[0]);
 	coord[1] = m_Op->GetDiscLine(1,pos[1]);

FDTD/extensions/operator_ext_pml_sf.cpp

@@ -227,7 +227,7 @@ void Operator_Ext_PML_SF_Plane::SetDirection(int ny, bool top_ny)
 	m_LineNr = (unsigned int)((int)m_top * (int)(m_Op->GetNumberOfLines(m_ny)-1));
 	pos[m_ny] = m_LineNr;
 
-	m_pml_delta = m_Op->GetMeshDelta(m_ny,pos);
+	m_pml_delta = m_Op->GetEdgeLength(m_ny,pos);
 }
 
 void Operator_Ext_PML_SF_Plane::SetPMLLength(int width)

FDTD/extensions/operator_ext_upml.cpp

@@ -277,17 +277,17 @@ void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm,
 			depth = width - (m_Op->GetDiscLine(n,pos[n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta();
 
 			if (n==ny)
-				depth-=m_Op->GetMeshDelta(n,pos)/2;
+				depth-=m_Op->GetEdgeLength(n,pos)/2;
 			double vars[5] = {depth, width/m_Size[2*n], width, Zm, m_Size[2*n]};
 			if (depth>0)
 				kappa_v[n] = m_GradingFunction->Eval(vars);
 			else
 				kappa_v[n]=0;
 			if (n==ny)
-				depth+=m_Op->GetMeshDelta(n,pos)/2;
+				depth+=m_Op->GetEdgeLength(n,pos)/2;
 
 			if (n!=ny)
-				depth-=m_Op->GetMeshDelta(n,pos)/2;
+				depth-=m_Op->GetEdgeLength(n,pos)/2;
 			if (depth<0)
 				depth=0;
 			vars[0]=depth;
@@ -302,17 +302,17 @@ void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm,
 			depth = width - (m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta();
 
 			if (n==ny)
-				depth+=m_Op->GetMeshDelta(n,pos)/2;
+				depth+=m_Op->GetEdgeLength(n,pos)/2;
 			double vars[5] = {depth, width/(m_Size[2*n]), width, Zm, m_Size[2*n]};
 			if (depth>0)
 				kappa_v[n] = m_GradingFunction->Eval(vars);
 			else
 				kappa_v[n]=0;
 			if (n==ny)
-				depth-=m_Op->GetMeshDelta(n,pos)/2;
+				depth-=m_Op->GetEdgeLength(n,pos)/2;
 
 			if (n!=ny)
-				depth+=m_Op->GetMeshDelta(n,pos)/2;
+				depth+=m_Op->GetEdgeLength(n,pos)/2;
 			if (depth>width)
 				depth=0;
 			vars[0]=depth;

FDTD/operator.cpp

@@ -97,24 +97,43 @@ void Operator::Reset()
 	Operator_Base::Reset();
 }
 
-double Operator::GetMeshDelta(int n, const unsigned int* pos, bool dualMesh) const
-{
-	if ((n<0) || (n>2)) return 0.0;
-	if (pos[n]>=numLines[n]) 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
 {
 	if ((n<0) || (n>2)) return 0.0;
 	if (pos>=numLines[n]) return 0.0;
 	if (dualMesh==false)
 		return discLines[n][pos];
+
+	// return dual mesh node
+	if (pos<numLines[n]-1)
+		return 0.5*(discLines[n][pos] + discLines[n][pos+1]);
+
+	// dual node for the last line (outside the field domain)
+	return discLines[n][pos] + 0.5*(discLines[n][pos] - discLines[n][pos-1]);
+
+}
+
+double Operator::GetEdgeLength(int n, const unsigned int* pos, bool dualMesh) const
+{
+	if ((n<0) || (n>2)) return 0.0;
+	if (pos[n]>=numLines[n]) return 0.0;
+	double delta=0;
+	if (dualMesh==false)
+	{
+		if (pos[n]<numLines[n]-1)
+			delta = GetDiscLine(n,pos[n]+1,false) - GetDiscLine(n,pos[n],false);
+		else
+			delta = GetDiscLine(n,pos[n],false) - GetDiscLine(n,pos[n]-1,false);
+		return delta*gridDelta;
+	}
 	else
-		return (discLines[n][pos] + 0.5*fabs(MainOp->GetIndexDelta(n,pos)));
+	{
+		if (pos[n]>0)
+			delta = GetDiscLine(n,pos[n],true) - GetDiscLine(n,pos[n]-1,true);
+		else
+			delta = GetDiscLine(n,1,false) - GetDiscLine(n,0,false);
+		return delta*gridDelta;
+	}
 }
 
 double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const
@@ -1242,7 +1261,7 @@ bool Operator::CalcFieldExcitation()
 						{
 							if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))//&& (pos[n]<numLines[n]-1))
 							{
-								amp = elec->GetWeightedExcitation(n,volt_coord)*GetMeshDelta(n,pos);// delta[n]*gridDelta;
+								amp = elec->GetWeightedExcitation(n,volt_coord)*GetEdgeLength(n,pos);// delta[n]*gridDelta;
 								if (amp!=0)
 								{
 									volt_vExcit.push_back(amp);
@@ -1284,7 +1303,7 @@ bool Operator::CalcFieldExcitation()
 						{
 							if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==2) || (elec->GetExcitType()==3) ))
 							{
-								amp = elec->GetWeightedExcitation(n,curr_coord)*GetMeshDelta(n,pos,true);// delta[n]*gridDelta;
+								amp = elec->GetWeightedExcitation(n,curr_coord)*GetEdgeLength(n,pos,true);// delta[n]*gridDelta;
 								if (amp!=0)
 								{
 									curr_vExcit.push_back(amp);

FDTD/operator.h

@@ -93,20 +93,18 @@ public:
 	virtual void ShowExtStat() const;
 
 	virtual double GetGridDelta() const {return gridDelta;}
-	//! Get the mesh delta times the grid delta for a 3D position (unit is meter)
-	virtual double GetMeshDelta(int n, const unsigned int* pos, bool dualMesh=false) const;
 
 	//! 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 position \a pos
-	virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetMeshDelta(ny,pos,!dualMesh);}
+	virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeLength(ny,pos,!dualMesh);}
 
 	//! 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 length of an FDTD edge (unit is meter).
-	virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetMeshDelta(ny,pos,dualMesh);}
+	virtual double GetEdgeLength(int ny, const unsigned 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
 	/*!

FDTD/operator_cylinder.cpp

@@ -85,68 +85,48 @@ string Operator_Cylinder::GetDirName(int ny) const
 	return "";
 }
 
-double Operator_Cylinder::GetMeshDelta(int n, const unsigned int* pos, bool dualMesh) const
-{
-	double delta = Operator_Multithread::GetMeshDelta(n,pos,dualMesh);
-	if (delta==0) return delta;
-	if (n==1)
-	{
-		return delta * GetDiscLine(0,pos[0],dualMesh);
-	}
-	return delta;
-}
-
 double Operator_Cylinder::GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh) const
 {
 	if ((ny<0) || (ny>2)) return 0.0;
 	if (pos[ny]>=numLines[ny]) return 0.0;
-	double width = 0;
+	double width = Operator_Multithread::GetEdgeLength(ny,pos,!dualMesh);
-	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 unsigned int pos[3], bool dualMesh) const
 {
 	if (pos[ny]>=numLines[ny]) return 0.0;
 	if (pos[0]>=numLines[0]) return 0.0;
 	if (ny==2)
 	{
-		double da = Operator_Multithread::GetMeshDelta(1,pos,dualMesh)/gridDelta;
+		double da = Operator_Multithread::GetEdgeLength(1,pos,dualMesh)/gridDelta;
 		double r1,r2;
-		if (!dualMesh)
+
+		if (dualMesh)
 		{
-			r1 = (discLines[0][pos[0]] - fabs(MainOp->GetIndexDelta(0,pos[0]-1))/2.0)*gridDelta;
+			r1 = GetDiscLine(0,pos[0],false)*gridDelta;
-			r2 = (discLines[0][pos[0]] + fabs(MainOp->GetIndexDelta(0,pos[0]  ))/2.0)*gridDelta;
+			r2 = r1 + GetEdgeLength(0,pos,false);
 		}
 		else
 		{
-			r1 =  discLines[0][pos[0]]*gridDelta;
+			r2 = GetDiscLine(0,pos[0],!dualMesh)*gridDelta;
-			r2 = (discLines[0][pos[0]] + fabs(MainOp->GetIndexDelta(0,pos[0])))*gridDelta;
+			r1 = r2 - GetEdgeLength(0,pos,true);
 		}
-		if (r1<0)
+
-			return da * pow(r2,2);
+		if (r1<=0)
-		return da/2* (pow(r2,2) - pow(r1,2));
+			return da/2 * pow(r2,2);
+		else
+			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 Operator_Multithread::GetNodeArea(ny,pos,dualMesh);
 }
 
 double Operator_Cylinder::GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh) const
 {
-	double length = Operator_Multithread::GetMeshDelta(ny,pos,dualMesh);
+	double length = Operator_Multithread::GetEdgeLength(ny,pos,dualMesh);
 	if (ny!=1)
 		return length;
 	return length * GetDiscLine(0,pos[0],dualMesh);

FDTD/operator_cylinder.h

@@ -44,16 +44,13 @@ public:
 	//! Get the name for the given direction: 0 -> rho, 1 -> alpha, 2 -> z
 	virtual string GetDirName(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, const 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;
 
 	//! 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 length of an FDTD edge.
+	//! 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;
 
 	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos