Revision: integral parameter processing (voltage & current calc)

revision necessary to keep the object oriented concept consistent for upcoming new integral parameter

FDTD/processcurrent.cpp

@@ -22,33 +22,16 @@
 
 ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : ProcessIntegral(op, eng)
 {
+	m_TimeShift = op->GetTimestep()/2.0;
+	m_dualMesh = true;
 }
 
 ProcessCurrent::~ProcessCurrent()
 {
 }
 
-void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
+double ProcessCurrent::CalcIntegral()
 {
-	if (Op->SnapToMesh(dstart,start,true,m_start_inside)==false)
-		cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
-	if (Op->SnapToMesh(dstop,stop,true,m_stop_inside)==false)
-		cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
-
-	if (g_settings.showProbeDiscretization()) {
-		cerr << m_Name << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], true ) << ","
-				<< Op->GetDiscLine( 1, start[1], true ) << "," << Op->GetDiscLine( 2, start[2], true ) << ") -> ("
-				<< Op->GetDiscLine( 0, stop[0], true ) << ","<< Op->GetDiscLine( 1, stop[1], true ) << ","
-				<< Op->GetDiscLine( 2, stop[2], true ) << ")";
-		cerr << "   [" << start[0] << "," << start[1] << "," << start[2] << "] -> ["
-				<< stop[0] << "," << stop[1] << "," << stop[2] << "]" << endl;
-	}
-}
-
-int ProcessCurrent::Process()
-{
-	if (Enabled==false) return -1;
-	if (CheckTimestep()==false) return GetNextInterval();
 	FDTD_FLOAT current=0;
 
 	int Dump_Dim = 0;
@@ -78,7 +61,7 @@ int ProcessCurrent::Process()
 	{
 		cerr << "ProcessCurrent::Process(): Warning Current Integration Box \"" << m_filename << "\" is not a surface (found dimension: " << Dump_Dim << ") --> i = 0" << endl;
 		current = 0;
-		return -1;
+		return 0.0;
 	}
 
 	switch (NormDir)
@@ -139,45 +122,9 @@ int ProcessCurrent::Process()
 		break;
 	default:
 		//this cannot happen...
-		return -2;
+		return 0.0;
 		break;
 	}
 
-	//	cerr << "ts: " << Eng->numTS << " i: " << current << endl;
+	return current;
-	current*=m_weight;
-
-
-	if (ProcessInterval)
-	{
-		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
-		{
-			TD_Values.push_back(current);
-			//current is sampled half a timestep later then the voltages
-			file  << setprecision(m_precision) << (0.5 + (double)Eng->GetNumberOfTimesteps())*Op->GetTimestep() << "\t" << current << endl;
-		}
-	}
-
-	if (m_FD_Interval)
-	{
-		if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0)
-		{
-			double T = ((double)Eng->GetNumberOfTimesteps() + 0.5) * Op->GetTimestep();
-			for (size_t n=0;n<m_FD_Samples.size();++n)
-			{
-				FD_Values.at(n) += current * cexp( -2.0 * 1.0i * M_PI * m_FD_Samples.at(n) * T );
-			}
-			++m_FD_SampleCount;
-			if (m_Flush)
-					FlushData();
-			m_Flush = false;
-		}
-	}
-
-	return GetNextInterval();
 }
-
-void ProcessCurrent::DumpBox2File( string vtkfilenameprefix, bool /*dualMesh*/ ) const
-{
-	Processing::DumpBox2File( vtkfilenameprefix, true );
-}
-

FDTD/processcurrent.h

@@ -26,11 +26,8 @@ public:
 	ProcessCurrent(Operator* op, Engine* eng);
 	virtual ~ProcessCurrent();
 
-	virtual void DefineStartStopCoord(double* dstart, double* dstop);
+	//! Integrate currents flowing through an area
-
+	virtual double CalcIntegral();
-	virtual int Process();
-
-	virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh = false ) const; //!< dump geometry to file
 
 protected:
 };

FDTD/processing.cpp

@@ -34,6 +34,7 @@ Processing::Processing(Operator* op, Engine* eng)
 	m_FD_Interval=0;
 	m_weight=1;
 	m_Flush = false;
+	m_dualMesh = false;
 }
 
 Processing::~Processing()
@@ -141,16 +142,16 @@ void Processing::AddFrequency(vector<double>  freqs)
 
 void Processing::DefineStartStopCoord(double* dstart, double* dstop)
 {
-	if (Op->SnapToMesh(dstart,start)==false)
+	if (Op->SnapToMesh(dstart,start,m_dualMesh,m_start_inside)==false)
 		cerr << "Processing::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
-	if (Op->SnapToMesh(dstop,stop)==false)
+	if (Op->SnapToMesh(dstop,stop,m_dualMesh,m_stop_inside)==false)
 		cerr << "Processing::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
 
 	if (g_settings.showProbeDiscretization()) {
-		cerr << m_Name << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], false ) << ","
+		cerr << m_Name << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], m_dualMesh ) << ","
-				<< Op->GetDiscLine( 1, start[1], false ) << "," << Op->GetDiscLine( 2, start[2], false ) << ") -> ("
+				<< Op->GetDiscLine( 1, start[1], m_dualMesh ) << "," << Op->GetDiscLine( 2, start[2], m_dualMesh ) << ") -> ("
-				<< Op->GetDiscLine( 0, stop[0], false ) << ","<< Op->GetDiscLine( 1, stop[1], false ) << ","
+				<< Op->GetDiscLine( 0, stop[0], m_dualMesh ) << ","<< Op->GetDiscLine( 1, stop[1], m_dualMesh ) << ","
-				<< Op->GetDiscLine( 2, stop[2], false ) << ")";
+				<< Op->GetDiscLine( 2, stop[2], m_dualMesh ) << ")";
 		cerr << "   [" << start[0] << "," << start[1] << "," << start[2] << "] -> ["
 				<< stop[0] << "," << stop[1] << "," << stop[2] << "]" << endl;
 	}

FDTD/processing.h

@@ -63,7 +63,11 @@ public:
 	//! Set the dump precision
 	void SetPrecision(unsigned int val) {m_precision = val;}
 
-	virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh = false ) const; //!< dump geometry to file
+	//! Dump probe geometry to file (will obay main or dual mesh property)
+	virtual void DumpBox2File( string vtkfilenameprefix) const {DumpBox2File(vtkfilenameprefix,m_dualMesh);}
+
+	//! Dump probe geometry to file
+	virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh) const;
 
 protected:
 	Operator* Op;
@@ -93,6 +97,9 @@ protected:
 
 	void Dump_FD_Data(vector<_Complex double> value, double factor, string filename);
 
+	//! define if given coords are on main or dualMesh (default is false)
+	bool m_dualMesh;
+
 	//! define/store snapped start/stop coords as mesh index
 	unsigned int start[3];
 	unsigned int stop[3];

FDTD/processintegral.cpp

@@ -16,9 +16,12 @@
 */
 
 #include "processintegral.h"
+#include <complex.h>
+#include <iomanip>
 
 ProcessIntegral::ProcessIntegral(Operator* op, Engine* eng)  : Processing(op, eng)
 {
+	m_TimeShift = 0.0;
 }
 
 ProcessIntegral::~ProcessIntegral()
@@ -42,3 +45,41 @@ void ProcessIntegral::FlushData()
 		Dump_FD_Data(FD_Values,1.0/(double)m_FD_SampleCount,m_filename + "_FD");
 }
 
+int ProcessIntegral::Process()
+{
+	if (Enabled==false) return -1;
+	if (CheckTimestep()==false) return GetNextInterval();
+
+	FDTD_FLOAT integral=CalcIntegral();
+	integral*=m_weight;
+
+	double time = (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() + m_TimeShift;
+
+	if (ProcessInterval)
+	{
+		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
+		{
+			TD_Values.push_back(integral);
+			file << setprecision(m_precision) << time << "\t" << integral << endl;
+		}
+	}
+
+	if (m_FD_Interval)
+	{
+		if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0)
+		{
+			double T = time;
+			for (size_t n=0;n<m_FD_Samples.size();++n)
+			{
+				FD_Values.at(n) += integral * cexp( -2.0 * 1.0i * M_PI * m_FD_Samples.at(n) * T );
+			}
+			++m_FD_SampleCount;
+			if (m_Flush)
+					FlushData();
+			m_Flush = false;
+		}
+	}
+
+	return GetNextInterval();
+}
+

FDTD/processintegral.h

@@ -27,13 +27,23 @@ public:
 	virtual ~ProcessIntegral();
 
 	virtual void InitProcess();
+
+	//! Flush FD data to disk
 	virtual void FlushData();
 
+	//! This method should calculate the integral parameter and must be overloaded for each derived class
+	virtual double CalcIntegral() {return 0;}
+
+	//! This method will write the TD and FD dump files using CalcIntegral() to calculate the integral parameter
+	virtual int Process();
+
 protected:
 	ProcessIntegral(Operator* op, Engine* eng);
 
-	vector<FDTD_FLOAT> TD_Values;
+	//! timeshift to be used in TD and FD data, e.g. 0.5*dT in case of current based parameter
+	double m_TimeShift;
 
+	vector<FDTD_FLOAT> TD_Values;
 	vector<_Complex double> FD_Values;
 };
 

FDTD/processvoltage.cpp

@@ -27,38 +27,8 @@ ProcessVoltage::~ProcessVoltage()
 {
 }
 
-int ProcessVoltage::Process()
+double ProcessVoltage::CalcIntegral()
 {
-	if (Enabled==false) return -1;
+	//integrate voltages from start to stop on a line
-	if (CheckTimestep()==false) return GetNextInterval();
+	return CalcLineIntegral(start,stop,0);
-
-	FDTD_FLOAT voltage=CalcLineIntegral(start,stop,0);
-	voltage*=m_weight;
-
-	if (ProcessInterval)
-	{
-		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
-		{
-			TD_Values.push_back(voltage);
-			file << setprecision(m_precision) << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl;
-		}
-	}
-
-	if (m_FD_Interval)
-	{
-		if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0)
-		{
-			double T = (double)Eng->GetNumberOfTimesteps() * Op->GetTimestep();
-			for (size_t n=0;n<m_FD_Samples.size();++n)
-			{
-				FD_Values.at(n) += voltage * cexp( -2.0 * 1.0i * M_PI * m_FD_Samples.at(n) * T );
-			}
-			++m_FD_SampleCount;
-			if (m_Flush)
-					FlushData();
-			m_Flush = false;
-		}
-	}
-
-	return GetNextInterval();
 }

FDTD/processvoltage.h

@@ -27,7 +27,7 @@ public:
 	ProcessVoltage(Operator* op, Engine* eng);
 	virtual ~ProcessVoltage();
 
-	virtual int Process();
+	virtual double CalcIntegral();
 
 protected:
 };