Added frequency domain probe support

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

FDTD/excitation.cpp

@@ -16,6 +16,7 @@
 */
 
 #include "tools/array_ops.h"
+#include "tools/useful.h"
 #include "fparser.hh"
 #include "tinyxml.h"
 #include "excitation.h"
@@ -100,13 +101,6 @@ bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS )
 	return true;
 }
 
-unsigned int Excitation::CalcNyquistNum(double fmax) const
-{
-	if (dT==0) return 1;
-	double T0 = 1/fmax;
-	return floor(T0/2/dT);
-}
-
 unsigned int Excitation::GetMaxExcitationTimestep() const
 {
 	FDTD_FLOAT maxAmp=0;
@@ -147,7 +141,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
 		Signal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
 	}
 
-	SetNyquistNum( CalcNyquistNum(f0+fc) );
+	SetNyquistNum( CalcNyquistNum(f0+fc,dT) );
 }
 
 void Excitation::CalcDiracPulsExcitation()
@@ -216,7 +210,7 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
 		Signal_curr[n] = fParse.Eval(vars);
 	}
 
-	SetNyquistNum( CalcNyquistNum(f0) );
+	SetNyquistNum( CalcNyquistNum(f0,dT) );
 }
 
 void Excitation::CalcSinusExcitation(double f0, int nTS)
@@ -240,7 +234,7 @@ void Excitation::CalcSinusExcitation(double f0, int nTS)
 		Signal_curr[n] = sin(2.0*PI*f0*t);
 	}
 
-	SetNyquistNum( CalcNyquistNum(f0) );
+	SetNyquistNum( CalcNyquistNum(f0,dT) );
 }
 
 void Excitation::setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit,

FDTD/excitation.h

@@ -34,7 +34,6 @@ public:
 
 	//! Get the excitation timestep with the (first) max amplitude
 	virtual unsigned int GetMaxExcitationTimestep() const;
-	unsigned int CalcNyquistNum(double fmax) const;
 
 	void setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit,
 								 vector<unsigned int> const& volt_vDelay, vector<unsigned int> const& volt_vDir );

FDTD/processcurrent.cpp

@@ -18,6 +18,7 @@
 #include "tools/global.h"
 #include "processcurrent.h"
 #include <iomanip>
+#include <complex.h>
 
 ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng)
 {
@@ -25,6 +26,7 @@ ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng)
 
 ProcessCurrent::~ProcessCurrent()
 {
+	ProcessCurrent::FlushData();
 }
 
 void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
@@ -44,6 +46,14 @@ void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
 	}
 }
 
+void ProcessCurrent::Init()
+{
+	FD_currents.clear();
+	for (size_t n=0;n<m_FD_Samples.size();++n)
+		FD_currents.push_back(0);
+}
+
+
 int ProcessCurrent::Process()
 {
 	if (Enabled==false) return -1;
@@ -144,9 +154,33 @@ int ProcessCurrent::Process()
 
 	//	cerr << "ts: " << Eng->numTS << " i: " << current << endl;
 	current*=m_weight;
+
+	if (ProcessInterval)
+	{
+		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
+		{
 			v_current.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_currents.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();
 }
 
@@ -154,3 +188,8 @@ void ProcessCurrent::DumpBox2File( string vtkfilenameprefix, bool /*dualMesh*/ )
 {
 	Processing::DumpBox2File( vtkfilenameprefix, true );
 }
+
+void ProcessCurrent::FlushData()
+{
+	Dump_FD_Data(FD_currents,1.0/(double)m_FD_SampleCount,m_filename + "_FD");
+}

FDTD/processcurrent.h

@@ -28,11 +28,16 @@ public:
 
 	virtual void DefineStartStopCoord(double* dstart, double* dstop);
 
+	virtual void Init();
+	virtual void FlushData();
 	virtual int Process();
 
 	virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh = false ) const; //!< dump geometry to file
 protected:
 	vector<FDTD_FLOAT> v_current;
+
+	vector<_Complex double> FD_currents;
+	void WriteFDCurrents();
 };
 
 #endif // PROCESSCURRENT_H

FDTD/processing.cpp

@@ -16,7 +16,10 @@
 */
 
 #include "tools/global.h"
+#include "tools/useful.h"
 #include "processing.h"
+#include <complex.h>
+#include <climits>
 
 Processing::Processing(Operator* op, Engine* eng)
 {
@@ -26,7 +29,10 @@ Processing::Processing(Operator* op, Engine* eng)
 	m_PS_pos = 0;
 	SetPrecision(12);
 	ProcessInterval=0;
+	m_FD_SampleCount=0;
+	m_FD_Interval=0;
 	m_weight=1;
+	m_Flush = false;
 }
 
 Processing::~Processing()
@@ -53,21 +59,36 @@ bool Processing::CheckTimestep()
 	{
 		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0) return true;
 	}
+
+	if (m_FD_Interval)
+	{
+		if (Eng->GetNumberOfTimesteps()%m_FD_Interval==0) return true;
+	}
 	return false;
 }
 
 int Processing::GetNextInterval() const
 {
 	if (Enabled==false) return -1;
-	unsigned int next=-1;
+	int next=INT_MAX;
 	if (m_ProcessSteps.size()>m_PS_pos)
 	{
-		next = m_ProcessSteps.at(m_PS_pos)-Eng->GetNumberOfTimesteps();
+		next = (int)m_ProcessSteps.at(m_PS_pos)-(int)Eng->GetNumberOfTimesteps();
 	}
-	if (ProcessInterval==0) return next;
+	if (ProcessInterval!=0)
-	unsigned int next_Interval = ProcessInterval - Eng->GetNumberOfTimesteps()%ProcessInterval;
+	{
+		int next_Interval = (int)ProcessInterval - (int)Eng->GetNumberOfTimesteps()%ProcessInterval;
 		if (next_Interval<next)
 			next = next_Interval;
+	}
+
+	//check for FD sample interval
+	if (m_FD_Interval!=0)
+	{
+		int next_Interval = (int)m_FD_Interval - (int)Eng->GetNumberOfTimesteps()%m_FD_Interval;
+		if (next_Interval<next)
+			next = next_Interval;
+	}
 	return next;
 }
 
@@ -87,6 +108,36 @@ void Processing::AddSteps(vector<unsigned int> steps)
 	}
 }
 
+void Processing::AddFrequency(double freq)
+{
+	unsigned int nyquistTS = CalcNyquistNum(freq,Op->GetTimestep());
+
+	if (nyquistTS == 0)
+	{
+		cerr << "Processing::AddFrequency: Requested frequency " << freq << " is too high for the current timestep used... skipping..." << endl;
+		return;
+	}
+	else if (nyquistTS<Op->Exc->GetNyquistNum())
+	{
+		cerr << "Processing::AddFrequency: Warning: Requested frequency " << freq << " is higher than maximum excited frequency..." << endl;
+	}
+
+	if (m_FD_Interval==0)
+		m_FD_Interval = Op->Exc->GetNyquistNum();
+	if (m_FD_Interval>nyquistTS)
+		m_FD_Interval = nyquistTS;
+
+	m_FD_Samples.push_back(freq);
+}
+
+void Processing::AddFrequency(vector<double>  freqs)
+{
+	for (size_t n=0;n<freqs.size();++n)
+	{
+		AddFrequency(freqs.at(n));
+	}
+}
+
 void Processing::DefineStartStopCoord(double* dstart, double* dstop)
 {
 	if (Op->SnapToMesh(dstart,start)==false)
@@ -222,11 +273,33 @@ void Processing::DumpBox2File( string vtkfilenameprefix, bool dualMesh ) const
 	file.close();
 }
 
+void Processing::Dump_FD_Data(vector<_Complex double> value, double factor, string filename)
+{
+	ofstream file;
+	file.open( filename.c_str() );
+	if (!file.is_open())
+		cerr << "Can't open file: " << filename << endl;
+
+	for (size_t n=0;n<value.size();++n)
+	{
+		file << m_FD_Samples.at(n) << "\t" << 2.0 * creal(value.at(n))*factor << "\t" << 2.0 * cimag(value.at(n))*factor << "\n";
+	}
+	file.close();
+}
+
 void ProcessingArray::AddProcessing(Processing* proc)
 {
 	ProcessArray.push_back(proc);
 }
 
+void ProcessingArray::FlushNext()
+{
+	for (size_t i=0;i<ProcessArray.size();++i)
+	{
+		ProcessArray.at(i)->FlushNext();
+	}
+}
+
 void ProcessingArray::Reset()
 {
 	for (size_t i=0;i<ProcessArray.size();++i)

FDTD/processing.h

@@ -20,6 +20,8 @@
 
 #include <iostream>
 #include <fstream>
+#include <complex>
+#include <cmath>
 #include "operator.h"
 #include "engine.h"
 
@@ -29,6 +31,7 @@ public:
 	Processing(Operator* op, Engine* eng);
 	virtual ~Processing();
 
+	virtual void Init() {};
 	virtual void Reset();
 
 	virtual void DefineStartStopCoord(double* dstart, double* dstop);
@@ -38,6 +41,9 @@ public:
 	void AddStep(unsigned int step);
 	void AddSteps(vector<unsigned int> steps);
 
+	void AddFrequency(double freq);
+	void AddFrequency(vector<double>  freqs);
+
 	bool CheckTimestep();
 	virtual int Process() {return GetNextInterval();}
 
@@ -49,17 +55,24 @@ public:
 	virtual void SetWeight(double weight) {m_weight=weight;}
 	virtual double GetWeight() {return m_weight;}
 
+	//! Invoke this flag to flush all stored data to disk
+	virtual void FlushNext() {m_Flush = true;}
+	virtual void FlushData() {};
+
 	//! Set the dump precision
 	void SetPrecision(unsigned int val) {m_precision = val;}
 
 	virtual void OpenFile( string outfile );
 
 	virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh = false ) const; //!< dump geometry to file
+
 protected:
 	Operator* Op;
 	Engine* Eng;
 	unsigned int m_precision;
 
+	bool m_Flush;
+
 	double m_weight;
 
 	bool Enabled;
@@ -70,6 +83,15 @@ protected:
 	size_t m_PS_pos; //! current position in list of processing steps
 	vector<unsigned int> m_ProcessSteps; //! list of processing steps
 
+	//! Vector of frequency samples
+	vector<double> m_FD_Samples;
+	//! Number of samples already processed
+	unsigned int m_FD_SampleCount;
+	//! Sampling interval needed for the FD_Samples
+	unsigned int m_FD_Interval;
+
+	void Dump_FD_Data(vector<_Complex double> value, double factor, string filename);
+
 	//! define/store snapped start/stop coords as mesh index
 	unsigned int start[3];
 	unsigned int stop[3];
@@ -94,6 +116,9 @@ public:
 
 	void AddProcessing(Processing* proc);
 
+	//! Invoke this flag to flush all stored data to disk for all processings on next Process()
+	void FlushNext();
+
 	void Reset();
 
 	//! Deletes all given processing's, can be helpful, but use carefull!!!

FDTD/processvoltage.cpp

@@ -16,6 +16,7 @@
 */
 
 #include "processvoltage.h"
+#include <complex.h>
 #include <iomanip>
 
 ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng)
@@ -24,16 +25,53 @@ ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng)
 
 ProcessVoltage::~ProcessVoltage()
 {
+	ProcessVoltage::FlushData();
+}
+
+void ProcessVoltage::Init()
+{
+	FD_voltages.clear();
+	for (size_t n=0;n<m_FD_Samples.size();++n)
+		FD_voltages.push_back(0);
 }
 
 int ProcessVoltage::Process()
 {
 	if (Enabled==false) return -1;
 	if (CheckTimestep()==false) return GetNextInterval();
+
 	FDTD_FLOAT voltage=CalcLineIntegral(start,stop,0);
-//	cerr << voltage << endl;
 	voltage*=m_weight;
+
+	if (ProcessInterval)
+	{
+		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
+		{
 			voltages.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_voltages.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();
 }
+
+void ProcessVoltage::FlushData()
+{
+	Dump_FD_Data(FD_voltages,1.0/(double)m_FD_SampleCount,m_filename + "_FD");
+}

FDTD/processvoltage.h

@@ -27,10 +27,16 @@ public:
 	ProcessVoltage(Operator* op, Engine* eng);
 	virtual ~ProcessVoltage();
 
+	virtual void Init();
+	virtual void FlushData();
+
 	virtual int Process();
 
 protected:
 	vector<FDTD_FLOAT> voltages;
+
+	vector<_Complex double> FD_voltages;
+	void WriteFDVoltages();
 };
 
 #endif // PROCESSVOLTAGE_H

openEMS.pro

@@ -56,7 +56,8 @@ SOURCES += main.cpp \
     FDTD/operator_sse_compressed.cpp \
     FDTD/engine_sse_compressed.cpp \
     FDTD/operator_multithread.cpp \
-    tools/global.cpp
+    tools/global.cpp \
+    tools/useful.cpp
 HEADERS += tools/ErrorMsg.h \
     tools/AdrOp.h \
     tools/constants.h \
@@ -83,7 +84,8 @@ HEADERS += tools/ErrorMsg.h \
     FDTD/operator_sse_compressed.h \
     FDTD/engine_sse_compressed.h \
     FDTD/operator_multithread.h \
-    tools/global.h
+    tools/global.h \
+    tools/useful.h
 QMAKE_CXXFLAGS_RELEASE = -O3 \
     -g \
 	-march=native

openems.cpp

@@ -330,8 +330,10 @@ int openEMS::SetupFDTD(const char* file)
 					continue;
 				}
 				proc->SetProcessInterval(Nyquist/m_OverSampling);
+				proc->AddFrequency(pb->GetFDSamples());
 				proc->DefineStartStopCoord(start,stop);
 				proc->SetWeight(pb->GetWeighting());
+				proc->Init();
 				PA->AddProcessing(proc);
 				prim->SetPrimitiveUsed(true);
 			}
@@ -467,6 +469,8 @@ void openEMS::RunFDTD()
 			cout << " --- Energy: ~" << setw(6) << setprecision(2) << std::scientific << currE << " (decrement: " << setw(6)  << setprecision(2) << std::fixed << fabs(10.0*log10(change)) << "dB)" << endl;
 			prevTime=currTime;
 			prevTS=currTS;
+
+			PA->FlushNext();
 		}
 	}
 

tools/useful.cpp

@@ -0,0 +1,31 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "useful.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <climits>
+
+unsigned int CalcNyquistNum(double fmax, double dT)
+{
+	if (fmax==0) return UINT_MAX;
+	if (dT==0) return 1;
+	double T0 = 1/fmax;
+	return floor(T0/2/dT);
+}
+

tools/useful.h

@@ -0,0 +1,24 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef USEFUL_H
+#define USEFUL_H
+
+//! Calc the nyquist number of timesteps for a given frequency and timestep
+unsigned int CalcNyquistNum(double fmax, double dT);
+
+#endif // USEFUL_H