added processing steps and fix in end criteria check, now checking max excite TS

FDTD/operator.cpp

@@ -230,6 +230,21 @@ void Operator::ShowStat() const
 	cout << "-----------------------------------" << endl;
 }
 
+unsigned int Operator::GetMaxExcitationTimestep() const
+{
+	FDTD_FLOAT maxAmp=0;
+	unsigned int maxStep=0;
+	for (unsigned int n=1;n<ExciteLength+1;++n)
+	{
+		if (fabs(ExciteSignal[n])>maxAmp)
+		{
+			maxAmp = fabs(ExciteSignal[n]);
+			maxStep = n;
+		}
+	}
+	return maxStep;
+}
+
 unsigned int Operator::CalcGaussianPulsExcitation(double f0, double fc)
 {
 	if (dT==0) return 0;

FDTD/operator.h

@@ -44,6 +44,9 @@ public:
 	//! Calculate a step excitation \return number of Nyquist timesteps
 	virtual unsigned int CalcStepExcitation();
 
+	//! Get the excitation timestep with the (first) max amplitude
+	virtual unsigned int GetMaxExcitationTimestep() const;
+
 	virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
 	virtual void ApplyMagneticBC(bool* dirs);
 

FDTD/processing.cpp

@@ -22,22 +22,65 @@ Processing::Processing(Operator* op, Engine* eng)
 	Op=op;
 	Eng=eng;
 	Enabled = true;
+	m_PS_pos = 0;
+	ProcessInterval=0;
 }
 
 Processing::~Processing()
 {
 }
 
+void Processing::Reset()
+{
+	m_PS_pos=0;
+}
+
 bool Processing::CheckTimestep()
+{
+	if (m_ProcessSteps.size()>m_PS_pos)
+	{
+		if (m_ProcessSteps.at(m_PS_pos)==Eng->GetNumberOfTimesteps())
+		{
+			++m_PS_pos;
+			return true;
+		}
+	}
+	if (ProcessInterval)
 	{
 		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0) return true;
+	}
 	return false;
 }
 
-int Processing::GetNextInterval()
+int Processing::GetNextInterval() const
 {
 	if (Enabled==false) return -1;
-	return ProcessInterval - Eng->GetNumberOfTimesteps()%ProcessInterval;
+	unsigned int next=-1;
+	if (m_ProcessSteps.size()>m_PS_pos)
+	{
+		next = m_ProcessSteps.at(m_PS_pos)-Eng->GetNumberOfTimesteps();
+	}
+	if (ProcessInterval==0) return next;
+	unsigned int next_Interval = ProcessInterval - Eng->GetNumberOfTimesteps()%ProcessInterval;
+	if (next_Interval<next)
+		next = next_Interval;
+	return next;
+}
+
+void Processing::AddStep(unsigned int step)
+{
+	if (m_ProcessSteps.size()==0)
+		m_ProcessSteps.push_back(step);
+	else if (find(m_ProcessSteps.begin(), m_ProcessSteps.end(),step)==m_ProcessSteps.end())
+		m_ProcessSteps.push_back(step);
+}
+
+void Processing::AddSteps(vector<unsigned int> steps)
+{
+	for (size_t n=0;n<steps.size();++n)
+	{
+		AddStep(steps.at(n));
+	}
 }
 
 void Processing::DefineStartStopCoord(double* dstart, double* dstop)
@@ -84,6 +127,14 @@ void ProcessingArray::AddProcessing(Processing* proc)
 	ProcessArray.push_back(proc);
 }
 
+void ProcessingArray::Reset()
+{
+	for (size_t i=0;i<ProcessArray.size();++i)
+	{
+		ProcessArray.at(i)->Reset();
+	}
+}
+
 void ProcessingArray::DeleteAll()
 {
 	for (size_t i=0;i<ProcessArray.size();++i)

FDTD/processing.h

@@ -29,15 +29,22 @@ public:
 	Processing(Operator* op, Engine* eng);
 	virtual ~Processing();
 
+	virtual void Reset();
+
 	virtual void DefineStartStopCoord(double* dstart, double* dstop);
 
 	void SetProcessInterval(unsigned int interval) {ProcessInterval=max((unsigned int)1,interval);}
+
+	void AddStep(unsigned int step);
+	void AddSteps(vector<unsigned int> steps);
+
+	bool CheckTimestep();
 	virtual int Process() {return GetNextInterval();}
 
 	//! If Disabled Process() will do nothing...
 	virtual void SetEnable(bool val) {Enabled=val;}
 	//! If Disabled Process() will do nothing...
-	virtual bool GetEnable() {return Enabled;}
+	virtual bool GetEnable() const {return Enabled;}
 
 protected:
 	Operator* Op;
@@ -45,10 +52,12 @@ protected:
 
 	bool Enabled;
 
-	bool CheckTimestep();
-	int GetNextInterval();
+	int GetNextInterval() const;
 	unsigned int ProcessInterval;
 
+	size_t m_PS_pos; //! current position in list of processing steps
+	vector<unsigned int> m_ProcessSteps; //! list of processing steps
+
 	unsigned int start[3];
 	unsigned int stop[3];
 
@@ -63,6 +72,8 @@ public:
 
 	void AddProcessing(Processing* proc);
 
+	void Reset();
+
 	//! Deletes all given processing's, can be helpful, but use carefull!!!
 	void DeleteAll();
 

openems.cpp

@@ -358,8 +358,16 @@ void openEMS::RunFDTD()
 {
 	cout << "Running FDTD engine... this may take a while... grab a cup of coffee?!?" << endl;
 
-	ProcessFields ProcField(FDTD_Op,FDTD_Eng);
+	//special handling of a field processing, needed to realize the end criteria...
+	ProcessFields* ProcField = new ProcessFields(FDTD_Op,FDTD_Eng);
+	PA->AddProcessing(ProcField);
 	double maxE=0,currE=0;
+
+	//add all timesteps to end-crit field processing with max excite amplitude
+	unsigned int maxExcite = FDTD_Op->GetMaxExcitationTimestep();
+	for (unsigned int n=0;n<FDTD_Op->E_Exc_Count;++n)
+		ProcField->AddStep(FDTD_Op->E_Exc_delay[n]+maxExcite);
+
 	double change=1;
 	int prevTS=0,currTS=0;
 	double speed = FDTD_Op->GetNumberCells()/1e6;
@@ -378,6 +386,14 @@ void openEMS::RunFDTD()
 	{
 		FDTD_Eng->IterateTS(step);
 		step=PA->Process();
+
+		if (ProcField->CheckTimestep())
+		{
+			currE = ProcField->CalcTotalEnergy();
+			if (currE>maxE)
+				maxE=currE;
+		}
+
 //		cout << " do " << step << " steps; current: " << eng.GetNumberOfTimesteps() << endl;
 		currTS = FDTD_Eng->GetNumberOfTimesteps();
 		if ((step<0) || (step>(int)(NrTS - currTS))) step=NrTS - currTS;
@@ -387,7 +403,7 @@ void openEMS::RunFDTD()
 		t_diff = CalcDiffTime(currTime,prevTime);
 		if (t_diff>4)
 		{
-			currE = ProcField.CalcTotalEnergy();
+			currE = ProcField->CalcTotalEnergy();
 			if (currE>maxE)
 				maxE=currE;
 			cout << "[@" << setw(8) << (int)CalcDiffTime(currTime,startTime)  <<  "s] Timestep: " << setw(12)  << currTS << " (" << setw(6) << setprecision(2) << std::fixed << (double)currTS/(double)NrTS*100.0  << "%)" ;