From 3ca8609020ab315e9904cbf23a1a26ab437780ac Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <Thorsten.Liebig@gmx.de>
Date: Fri, 9 Apr 2010 15:58:15 +0200
Subject: [PATCH] fix: don't process current integration outside of simulation
 domain

---
 FDTD/operator.cpp       |  8 +++++---
 FDTD/operator.h         |  2 +-
 FDTD/processcurrent.cpp | 37 +++++++++++++++++++++++--------------
 FDTD/processing.cpp     |  4 ++--
 FDTD/processing.h       |  5 +++++
 5 files changed, 36 insertions(+), 20 deletions(-)

diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index cf1ebac..5c7eed1 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -102,15 +102,17 @@ unsigned int Operator::CalcNyquistNum(double fmax)
 	return floor(T0/2/dT);
 }
 
-bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower)
+bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, bool* inside)
 {
 	bool ok=true;
 	for (int n=0;n<3;++n)
 	{
+		if (inside) //set defaults
+			inside[n] = true;
 		uicoord[n]=0;
-		if (dcoord[n]<discLines[n][0]) {ok=false;uicoord[n]=0;}
+		if (dcoord[n]<discLines[n][0]) {ok=false;uicoord[n]=0; if (inside) inside[n] = false;}
 		else if (dcoord[n]==discLines[n][0]) {uicoord[n]=0;}
-		else if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
+		else if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2; if (inside) inside[n] = false; }
 		else if (dcoord[n]==discLines[n][numLines[n]-1]) {uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
 		else
 			for (unsigned int i=1;i<numLines[n];++i)
diff --git a/FDTD/operator.h b/FDTD/operator.h
index fd4a6eb..623cc2c 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -67,7 +67,7 @@ public:
 
 	virtual void Reset();
 
-	bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false);
+	bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL);
 
 protected:
 	//! use New() for creating a new Operator
diff --git a/FDTD/processcurrent.cpp b/FDTD/processcurrent.cpp
index a63a997..c951fd4 100644
--- a/FDTD/processcurrent.cpp
+++ b/FDTD/processcurrent.cpp
@@ -39,8 +39,8 @@ void ProcessCurrent::OpenFile(string outfile)
 
 void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
 {
-	if (Op->SnapToMesh(dstart,start,true)==false) cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl;
-	if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
+	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;
 }
 
 int ProcessCurrent::Process()
@@ -62,27 +62,36 @@ int ProcessCurrent::Process()
 			unsigned int help=start[n];
 			start[n]=stop[n];
 			stop[n]=help;
+			bool b_help=m_start_inside[n];
+			m_start_inside[n] = m_stop_inside[n];
+			m_stop_inside[n] = b_help;
 		}
 	}
 
 	//x-current
-	for (unsigned int i=start[0];i<stop[0];++i)
-		current+=curr[0][i][start[1]][start[2]];
+	if (m_start_inside[1] && m_start_inside[2])
+		for (unsigned int i=start[0];i<stop[0];++i)
+			current+=curr[0][i][start[1]][start[2]];
 	//y-current
-	for (unsigned int i=start[1];i<stop[1];++i)
-		current+=curr[1][stop[0]][i][start[2]];
+	if (m_stop_inside[0] && m_start_inside[2])
+		for (unsigned int i=start[1];i<stop[1];++i)
+			current+=curr[1][stop[0]][i][start[2]];
 	//z-current
-	for (unsigned int i=start[2];i<stop[2];++i)
-		current+=curr[2][stop[0]][stop[1]][i];
+	if (m_stop_inside[0] && m_stop_inside[1])
+		for (unsigned int i=start[2];i<stop[2];++i)
+			current+=curr[2][stop[0]][stop[1]][i];
 	//x-current
-	for (unsigned int i=start[0];i<stop[0];++i)
-		current-=curr[0][i][stop[1]][stop[2]];
+	if (m_stop_inside[1] && m_stop_inside[2])
+		for (unsigned int i=start[0];i<stop[0];++i)
+			current-=curr[0][i][stop[1]][stop[2]];
 	//y-current
-	for (unsigned int i=start[1];i<stop[1];++i)
-		current-=curr[1][start[0]][i][stop[2]];
+	if (m_start_inside[0] && m_stop_inside[2])
+		for (unsigned int i=start[1];i<stop[1];++i)
+			current-=curr[1][start[0]][i][stop[2]];
 	//z-current
-	for (unsigned int i=start[2];i<stop[2];++i)
-		current-=curr[2][start[0]][start[1]][i];
+	if (m_start_inside[0] && m_start_inside[1])
+		for (unsigned int i=start[2];i<stop[2];++i)
+			current-=curr[2][start[0]][start[1]][i];
 
 //	cerr << "ts: " << Eng->numTS << " i: " << current << endl;
 	v_current.push_back(current);
diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp
index 4f25a3f..699626d 100644
--- a/FDTD/processing.cpp
+++ b/FDTD/processing.cpp
@@ -85,8 +85,8 @@ void Processing::AddSteps(vector<unsigned int> steps)
 
 void Processing::DefineStartStopCoord(double* dstart, double* dstop)
 {
-	if (Op->SnapToMesh(dstart,start)==false) cerr << "Processing::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl;
-	if (Op->SnapToMesh(dstop,stop)==false) cerr << "Processing::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
+	if (Op->SnapToMesh(dstart,start)==false) cerr << "Processing::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
+	if (Op->SnapToMesh(dstop,stop)==false) cerr << "Processing::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
 }
 
 double Processing::CalcLineIntegral(unsigned int* start, unsigned int* stop, int field)
diff --git a/FDTD/processing.h b/FDTD/processing.h
index 9286307..cdf6259 100644
--- a/FDTD/processing.h
+++ b/FDTD/processing.h
@@ -58,9 +58,14 @@ protected:
 	size_t m_PS_pos; //! current position in list of processing steps
 	vector<unsigned int> m_ProcessSteps; //! list of processing steps
 
+	//! define/store snapped start/stop coords as mesh index
 	unsigned int start[3];
 	unsigned int stop[3];
 
+	//! define/store if snapped start/stop coords are inside the field domain
+	bool m_start_inside[3];
+	bool m_stop_inside[3];
+
 	double CalcLineIntegral(unsigned int* start, unsigned int* stop, int field);
 };