updated debug dumps

The debug dump routines do not work correctly in cylinder coords (--nativeFieldDumps may help)!

FDTD/operator.cpp

@@ -345,25 +345,47 @@ void Operator::DumpOperator2File(string filename)
 #endif
 
 	ofstream file(filename.c_str(),ios_base::out);
-	if (file.is_open()==false)
+	if (!file.is_open())
 	{
-		cerr << "Operator::DumpOperator2File: Can't open file: " << filename << endl;
+		cerr << "Operator::DumpOperator2File(): Can't open file: " << filename << endl;
 		return;
 	}
 
-	cout << "Operator::Dumping FDTD operator information to vtk file: " << filename << " ..." << flush ;
+	cout << "Operator: Dumping FDTD operator information to vtk file: " << filename << " ..." << flush;
 
 	FDTD_FLOAT**** exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
-	if (Exc) {
+	if (Exc)
+	{
 		for (unsigned int n=0;n<Exc->Volt_Count;++n)
 			exc[Exc->Volt_dir[n]][Exc->Volt_index[0][n]][Exc->Volt_index[1][n]][Exc->Volt_index[2][n]] = Exc->Volt_amp[n];
 	}
 
+	FDTD_FLOAT**** vv_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** vi_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** iv_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** ii_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+
+	unsigned int pos[3], n;
+	for (n=0; n<3; n++)
+		for (pos[0]=0; pos[0]<numLines[0]; pos[0]++)
+			for (pos[1]=0; pos[1]<numLines[1]; pos[1]++)
+				for (pos[2]=0; pos[2]<numLines[2]; pos[2]++)
+				{
+					vv_temp[n][pos[0]][pos[1]][pos[2]] = GetVV(n,pos);
+					vi_temp[n][pos[0]][pos[1]][pos[2]] = GetVI(n,pos);
+					iv_temp[n][pos[0]][pos[1]][pos[2]] = GetIV(n,pos);
+					ii_temp[n][pos[0]][pos[1]][pos[2]] = GetII(n,pos);
+				}
+
 	string names[] = {"vv", "vi", "iv" , "ii", "exc"};
-	FDTD_FLOAT**** array[] = {vv,vi,iv,ii,exc};
+	FDTD_FLOAT**** array[] = {vv_temp,vi_temp,iv_temp,ii_temp,exc};
 
 	ProcessFields::DumpMultiVectorArray2VTK(file, names , array , 5, discLines, numLines, 6, "Operator dump" , (ProcessFields::MeshType)m_MeshType, discLines_scaling);
 
+	Delete_N_3DArray(ii_temp,numLines);
+	Delete_N_3DArray(iv_temp,numLines);
+	Delete_N_3DArray(vi_temp,numLines);
+	Delete_N_3DArray(vv_temp,numLines);
 	Delete_N_3DArray(exc,numLines);
 
 	file.close();
@@ -377,12 +399,13 @@ void Operator::DumpOperator2File(string filename)
 void Operator::DumpPEC2File( string filename )
 {
 	ofstream file( filename.c_str() );
-	if (!file.is_open()) {
+	if (!file.is_open())
-		cerr << "Operator::DumpPEC2File: Can't open file: " << filename << endl;
+	{
+		cerr << "Operator::DumpPEC2File(): Can't open file: " << filename << endl;
 		return;
 	}
 
-	cout << "Dumping PEC information to vtk file: " << filename << " ..." << flush;
+	cout << "Operator: Dumping PEC information to vtk file: " << filename << " ..." << flush;
 
 	FDTD_FLOAT**** pec = Create_N_3DArray<FDTD_FLOAT>( numLines );
 	unsigned int pos[3];
@@ -399,19 +422,19 @@ void Operator::DumpPEC2File( string filename )
 				if ((pos[1] != 0) && (pos[2] != 0))
 				{
 					// PEC surrounds the computational area; do not output this
-					if ((GetVV(0,pos[0],pos[1],pos[2]) == 0) && (GetVI(0,pos[0],pos[1],pos[2]) == 0))
+					if ((GetVV(0,pos) == 0) && (GetVI(0,pos) == 0))
 						pec[0][pos[0]][pos[1]][pos[2]] = GetEdgeLength( 0, pos ) * scaling; // PEC-x found
 				}
 				if ((pos[0] != 0) && (pos[2] != 0))
 				{
 					// PEC surrounds the computational area; do not output this
-					if ((GetVV(1,pos[0],pos[1],pos[2]) == 0) && (GetVI(1,pos[0],pos[1],pos[2]) == 0))
+					if ((GetVV(1,pos) == 0) && (GetVI(1,pos) == 0))
 						pec[1][pos[0]][pos[1]][pos[2]] = GetEdgeLength( 1, pos ) * scaling; // PEC-y found
 				}
 				if ((pos[0] != 0) && (pos[1] != 0))
 				{
 					// PEC surrounds the computational area; do not output this
-					if ((GetVV(2,pos[0],pos[1],pos[2]) == 0) && (GetVI(2,pos[0],pos[1],pos[2]) == 0))
+					if ((GetVV(2,pos) == 0) && (GetVI(2,pos) == 0))
 						pec[2][pos[0]][pos[1]][pos[2]] = GetEdgeLength( 2, pos ) * scaling; // PEC-z found
 				}
 			}
@@ -463,25 +486,20 @@ void Operator::DumpMaterial2File(string filename)
 #endif
 
 	ofstream file(filename.c_str(),ios_base::out);
-	if (file.is_open()==false)
+	if (!file.is_open())
 	{
-		cerr << "Operator::DumpMaterial2File: Can't open file: " << filename << endl;
+		cerr << "Operator::DumpMaterial2File(): Can't open file: " << filename << endl;
 		return;
 	}
 
-	cout << "Dumping material information to vtk file: " << filename << " ..."  << flush;
+	cout << "Operator: Dumping material information to vtk file: " << filename << " ..."  << flush;
+
+	FDTD_FLOAT**** epsilon = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** mue     = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** kappa   = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** sigma   = Create_N_3DArray<FDTD_FLOAT>(numLines);
 
-	FDTD_FLOAT**** epsilon;
-	FDTD_FLOAT**** mue;
-	FDTD_FLOAT**** kappa;
-	FDTD_FLOAT**** sigma;
 	unsigned int pos[3];
-	double inMat[4];
-
-	epsilon = Create_N_3DArray<FDTD_FLOAT>( numLines);
-	mue = Create_N_3DArray<FDTD_FLOAT>( numLines);
-	kappa = Create_N_3DArray<FDTD_FLOAT>( numLines);
-	sigma = Create_N_3DArray<FDTD_FLOAT>( numLines);
 	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 	{
 		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
@@ -490,11 +508,12 @@ void Operator::DumpMaterial2File(string filename)
 			{
 				for (int n=0;n<3;++n)
 				{
+					double inMat[4];
 					Calc_EffMatPos(n, pos, inMat);
-					epsilon[n][pos[0]][pos[1]][pos[2]]	=inMat[0]/__EPS0__;
+					epsilon[n][pos[0]][pos[1]][pos[2]] = inMat[0]/__EPS0__;
-					mue[n][pos[0]][pos[1]][pos[2]]		=inMat[2]/__MUE0__;
+					mue[n][pos[0]][pos[1]][pos[2]]     = inMat[2]/__MUE0__;
-					kappa[n][pos[0]][pos[1]][pos[2]]	=inMat[1];
+					kappa[n][pos[0]][pos[1]][pos[2]]   = inMat[1];
-					sigma[n][pos[0]][pos[1]][pos[2]]	=inMat[3];
+					sigma[n][pos[0]][pos[1]][pos[2]]   = inMat[3];
 				}
 			}
 		}
@@ -502,11 +521,14 @@ void Operator::DumpMaterial2File(string filename)
 
 	string names[] = {"epsilon","mue","kappa","sigma"};
 	FDTD_FLOAT**** array[] = {epsilon,mue,kappa,sigma};
+
 	ProcessFields::DumpMultiVectorArray2VTK(file, names, array, 4, discLines, numLines,  6, "Material dump" , (ProcessFields::MeshType)m_MeshType, discLines_scaling);
+
 	Delete_N_3DArray(epsilon,numLines);
 	Delete_N_3DArray(mue,numLines);
 	Delete_N_3DArray(kappa,numLines);
 	Delete_N_3DArray(sigma,numLines);
+
 	file.close();
 
 	cout << " done!" << endl;
@@ -577,7 +599,7 @@ void Operator::Calc_ECOperatorPos(int n, unsigned int* pos)
 	}
 }
 
-int Operator::CalcECOperator()
+int Operator::CalcECOperator( DebugFlags debugFlags )
 {
 	Init_EC();
 
@@ -643,6 +665,13 @@ int Operator::CalcECOperator()
 	for (size_t n=0;n<m_Op_exts.size();++n)
 		m_Op_exts.at(n)->BuildExtension();
 
+	if (debugFlags & debugMaterial)
+		DumpMaterial2File( "material_dump.vtk" );
+	if (debugFlags & debugOperator)
+		DumpOperator2File( "operator_dump.vtk" );
+	if (debugFlags & debugPEC)
+		DumpPEC2File( "PEC_dump.vtk" );
+
 	//cleanup
 	for (int n=0;n<3;++n)
 	{

FDTD/operator.h

@@ -33,7 +33,9 @@ class Operator
 	friend class Engine;
 	friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice
 	friend class Operator_Ext_PML_SF_Plane;
-public:
+public:	
+	enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4};
+
 	//! Create a new operator
 	static Operator* New();
 	virtual ~Operator();
@@ -43,7 +45,7 @@ public:
 	virtual bool SetGeometryCSX(ContinuousStructure* geo);
 	virtual ContinuousStructure* GetGeometryCSX() {return CSX;}
 
-	virtual int CalcECOperator();
+	virtual int CalcECOperator( DebugFlags );
 
 	//! Calculate the FDTD equivalent circuit parameter for the given position and direction ny. \sa Calc_EffMat_Pos
 	virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC) const;
@@ -59,6 +61,12 @@ public:
 	inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[n][x][y][z]; }
 	inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[n][x][y][z]; }
 
+	// convenient access functions
+	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int pos[3] ) const { return GetVV(n,pos[0],pos[1],pos[2]); }
+	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int pos[3] ) const { return GetVI(n,pos[0],pos[1],pos[2]); }
+	inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int pos[3] ) const { return GetII(n,pos[0],pos[1],pos[2]); }
+	inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int pos[3] ) const { return GetIV(n,pos[0],pos[1],pos[2]); }
+
 	// the next four functions need to be reimplemented in a derived class
 	inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vv[n][x][y][z] = value; }
 	inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vi[n][x][y][z] = value; }
@@ -83,10 +91,6 @@ public:
 	virtual void ShowStat() const;
 	virtual void ShowExtStat() const;
 
-	virtual void DumpOperator2File(string filename);
-	virtual void DumpMaterial2File(string filename);
-	virtual void DumpPEC2File( string filename );
-
 	//! Get the name for the given direction: 0 -> x, 1 -> y, 2 -> z
 	virtual string GetDirName(int ny) const;
 
@@ -142,6 +146,11 @@ protected:
 	//! Calculate the field excitations.
 	virtual bool CalcFieldExcitation();
 
+	// debug
+	virtual void DumpOperator2File(string filename);
+	virtual void DumpMaterial2File(string filename);
+	virtual void DumpPEC2File( string filename );
+
 	unsigned int m_Nr_PEC[3]; //count PEC edges
 	virtual bool CalcPEC();
 	virtual void CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned int* counter);	//internal to CalcPEC
@@ -188,4 +197,7 @@ public:
 	Excitation* Exc;
 };
 
+inline Operator::DebugFlags operator|( Operator::DebugFlags a, Operator::DebugFlags b ) { return static_cast<Operator::DebugFlags>(static_cast<int>(a) | static_cast<int>(b)); }
+inline Operator::DebugFlags& operator|=( Operator::DebugFlags& a, const Operator::DebugFlags& b ) { return a = a | b; }
+
 #endif // OPERATOR_H

FDTD/operator_cylindermultigrid.cpp

@@ -146,25 +146,25 @@ void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, ve
 }
 
 
-int Operator_CylinderMultiGrid::CalcECOperator()
+int Operator_CylinderMultiGrid::CalcECOperator( DebugFlags debugFlags )
 {
 	int retCode=0;
 	if (dT)
 		m_InnerOp->SetTimestep(dT);
 
 	//calc inner child first
-	m_InnerOp->CalcECOperator();
+	m_InnerOp->CalcECOperator( debugFlags );
 
 	dT = m_InnerOp->GetTimestep();
 
-	retCode = Operator_Cylinder::CalcECOperator();
+	retCode = Operator_Cylinder::CalcECOperator( debugFlags );
 	if (GetTimestepValid()==false)
 	{
 		cerr << "Operator_CylinderMultiGrid::CalcECOperator(): Warning, timestep invalid... resetting..." << endl;
 		dT = opt_dT;
 		m_InnerOp->SetTimestep(dT);
-		m_InnerOp->CalcECOperator();
+		m_InnerOp->CalcECOperator( debugFlags );
-		return Operator_Cylinder::CalcECOperator();
+		return Operator_Cylinder::CalcECOperator( debugFlags );
 	}
 
 	return retCode;

FDTD/operator_cylindermultigrid.h

@@ -41,7 +41,7 @@ public:
 
 	virtual unsigned int GetSplitPos() const {return m_Split_Pos;}
 
-	virtual int CalcECOperator();
+	virtual int CalcECOperator( DebugFlags );
 
 	virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS);
 

FDTD/operator_multithread.cpp

@@ -94,7 +94,7 @@ void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector<u
 	}
 }
 
-int Operator_Multithread::CalcECOperator()
+int Operator_Multithread::CalcECOperator( DebugFlags debugFlags )
 {
 	if (m_numThreads == 0)
 		m_numThreads = boost::thread::hardware_concurrency();
@@ -118,7 +118,7 @@ int Operator_Multithread::CalcECOperator()
 		m_thread_group.add_thread( t );
 	}
 
-	return Operator_SSE_Compressed::CalcECOperator();
+	return Operator_SSE_Compressed::CalcECOperator( debugFlags );
 }
 
 bool Operator_Multithread::Calc_EC()

FDTD/operator_multithread.h

@@ -31,7 +31,7 @@ public:
 	static Operator_Multithread* New(unsigned int numThreads = 0);
 	virtual ~Operator_Multithread();
 
-	virtual int CalcECOperator();
+	virtual int CalcECOperator( DebugFlags );
 
 	virtual void setNumThreads( unsigned int numThreads );
 

FDTD/operator_sse.cpp

@@ -83,67 +83,3 @@ void Operator_sse::InitOperator()
 
 	numVectors =  ceil((double)numLines[2]/4.0);
 }
-
-void Operator_sse::DumpOperator2File(string filename)
-{
-#ifdef OUTPUT_IN_DRAWINGUNITS
-	double discLines_scaling = 1;
-#else
-	double discLines_scaling = GetGridDelta();
-#endif
-
-	ofstream file(filename.c_str(),ios_base::out);
-	if (file.is_open()==false)
-	{
-		cerr << "Operator_sse::DumpOperator2File: Can't open file: " << filename << endl;
-		return;
-	}
-
-	cout << "Operator_sse::Dumping FDTD operator information to vtk file: " << filename << " ..." << flush ;
-
-
-	FDTD_FLOAT**** exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
-	if (Exc) {
-		for (unsigned int n=0;n<Exc->Volt_Count;++n)
-			exc[Exc->Volt_dir[n]][Exc->Volt_index[0][n]][Exc->Volt_index[1][n]][Exc->Volt_index[2][n]] = Exc->Volt_amp[n];
-	}
-
-	vv = Create_N_3DArray<FDTD_FLOAT>(numLines);
-	vi = Create_N_3DArray<FDTD_FLOAT>(numLines);
-	iv = Create_N_3DArray<FDTD_FLOAT>(numLines);
-	ii = Create_N_3DArray<FDTD_FLOAT>(numLines);
-
-	unsigned int pos[3];
-	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
-	{
-		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
-		{
-			for (pos[2]=0;pos[2]<numLines[2];++pos[2])
-			{
-				for (int n=0;n<3;++n)
-				{
-					vv[n][pos[0]][pos[1]][pos[2]] = GetVV(n,pos[0],pos[1],pos[2]);
-					vi[n][pos[0]][pos[1]][pos[2]] = GetVI(n,pos[0],pos[1],pos[2]);
-					ii[n][pos[0]][pos[1]][pos[2]] = GetII(n,pos[0],pos[1],pos[2]);
-					iv[n][pos[0]][pos[1]][pos[2]] = GetIV(n,pos[0],pos[1],pos[2]);
-				}
-			}
-		}
-	}
-
-	string names[] = {"vv", "vi", "iv" , "ii", "exc"};
-	FDTD_FLOAT**** array[] = {vv,vi,iv,ii,exc};
-
-	ProcessFields::DumpMultiVectorArray2VTK(file, names , array , 5, discLines, numLines,  6, "Operator dump" , (ProcessFields::MeshType)m_MeshType, discLines_scaling);
-
-	Delete_N_3DArray(exc,numLines);
-	Delete_N_3DArray(vv,numLines);vv=NULL;
-	Delete_N_3DArray(vi,numLines);vi=NULL;
-	Delete_N_3DArray(iv,numLines);iv=NULL;
-	Delete_N_3DArray(ii,numLines);ii=NULL;
-
-	file.close();
-
-	cout << " done!" << endl;
-}
-

FDTD/operator_sse.h

@@ -30,8 +30,6 @@ public:
 
 	virtual Engine* CreateEngine() const;
 
-	virtual void DumpOperator2File(string filename);
-
 	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vv[n][x][y][z%numVectors].f[z/numVectors]; }
 	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vi[n][x][y][z%numVectors].f[z/numVectors]; }
 	inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_ii[n][x][y][z%numVectors].f[z/numVectors]; }

FDTD/operator_sse_compressed.cpp

@@ -54,21 +54,10 @@ Engine* Operator_SSE_Compressed::CreateEngine() const
 	return eng;
 }
 
-void Operator_SSE_Compressed::DumpOperator2File(string filename)
+int Operator_SSE_Compressed::CalcECOperator( DebugFlags debugFlags )
-{
-	if (m_Use_Compression)
-	{
-		cerr << "Operator_SSE_Compressed::DumpOperator2File: Warning: Operator dump not implemented for a compressed operator yet, try disabling operator compression." << endl;
-		return;
-	}
-	else
-		Operator_sse::DumpOperator2File(filename);
-}
-
-int Operator_SSE_Compressed::CalcECOperator()
 {
+	Operator_sse::CalcECOperator( debugFlags );
 	m_Use_Compression = false;
-	Operator_sse::CalcECOperator();
 	m_Use_Compression = CompressOperator();
 
 	return 0;

FDTD/operator_sse_compressed.h

@@ -45,9 +45,7 @@ public:
 
 	virtual Engine* CreateEngine() const;
 
-	virtual void DumpOperator2File(string filename);
+	virtual int CalcECOperator( DebugFlags );
-
-	virtual int CalcECOperator();
 
 	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVV(n,x,y,z);}
 	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVI(n,x,y,z);}

openems.cpp

@@ -385,7 +385,14 @@ int openEMS::SetupFDTD(const char* file)
 	if (timestep)
 		FDTD_Op->SetTimestep(timestep);
 
-	FDTD_Op->CalcECOperator();
+	Operator::DebugFlags debugFlags = Operator::None;
+	if (DebugMat)
+		debugFlags |= Operator::debugMaterial;
+	if (DebugOp)
+		debugFlags |= Operator::debugOperator;
+	if (m_debugPEC)
+		debugFlags |= Operator::debugPEC;
+	FDTD_Op->CalcECOperator( debugFlags );
 	
 	unsigned int maxTime_TS = (unsigned int)(maxTime/FDTD_Op->GetTimestep());
 	if ((maxTime_TS>0) && (maxTime_TS<NrTS))
@@ -394,14 +401,6 @@ int openEMS::SetupFDTD(const char* file)
 	if (!FDTD_Op->SetupExcitation( FDTD_Opts->FirstChildElement("Excitation"), NrTS ))
 		exit(2);
 
-	// create debug output, if requested
-	if (DebugMat)
-		FDTD_Op->DumpMaterial2File("material_dump.vtk");
-	if (DebugOp)
-		FDTD_Op->DumpOperator2File("operator_dump.vtk");
-	if (m_debugPEC)
-		FDTD_Op->DumpPEC2File("PEC_dump.vtk");
-
 	timeval OpDoneTime;
 	gettimeofday(&OpDoneTime,NULL);