From 610f4e62da076d7002c97921a7b7cda1046ce05d Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <Thorsten.Liebig@gmx.de>
Date: Wed, 2 Jun 2010 14:18:25 +0200
Subject: [PATCH] processfields: allow different mesh types for vtk dump, e.g.
 cylindrical mesh dump

---
 FDTD/processfields.cpp    | 52 ++++++++++++++++++++++++++++++++-------
 FDTD/processfields.h      | 19 +++++++++-----
 FDTD/processfields_td.cpp | 12 ++++-----
 openems.cpp               |  2 ++
 4 files changed, 64 insertions(+), 21 deletions(-)

diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp
index cd01c6a..3115fe5 100644
--- a/FDTD/processfields.cpp
+++ b/FDTD/processfields.cpp
@@ -26,6 +26,7 @@ ProcessFields::ProcessFields(Operator* op, Engine* eng) : Processing(op, eng)
 	m_DumpType = E_FIELD_DUMP;
 	// vtk-file is default
 	m_fileType = VTK_FILETYPE;
+	m_Mesh_Type = CARTESIAN_MESH;
 	SetSubSampling(1);
 	SetPrecision(6);
 
@@ -254,7 +255,17 @@ void  ProcessFields::SetSubSampling(unsigned int subSampleRate, int dir)
 	else subSample[dir]=subSampleRate;
 }
 
-void ProcessFields::WriteVTKHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info)
+void ProcessFields::WriteVTKHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT)
+{
+	if (meshT==CARTESIAN_MESH)
+		WriteVTKCartesianGridHeader(file, discLines, numLines, precision, header_info);
+	else if (meshT==CYLINDRICAL_MESH)
+		WriteVTKCylindricalGridHeader(file, discLines, numLines, precision, header_info);
+	else
+		cerr << "ProcessFields::WriteVTKHeader: Warning: unknown mesh type, skipping header -> file will be invalid..." << endl;
+}
+
+void ProcessFields::WriteVTKCartesianGridHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info)
 {
 	file << "# vtk DataFile Version 2.0" << endl;
 	file << "Rectilinear Grid openEMS_ProcessFields";
@@ -279,6 +290,29 @@ void ProcessFields::WriteVTKHeader(ofstream &file, double const* const* discLine
 	file << "POINT_DATA " << numLines[0]*numLines[1]*numLines[2] << endl;
 }
 
+void ProcessFields::WriteVTKCylindricalGridHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info)
+{
+	file << "# vtk DataFile Version 3.0" << endl;
+	file << "Structured Grid from openEMS_ProcessFields";
+	if (!header_info.empty())
+		file << " " << header_info;
+	file << endl;
+	file << "ASCII" << endl;
+	file << "DATASET STRUCTURED_GRID " << endl;
+	file << "DIMENSIONS " << numLines[0] << " " << numLines[1] << " " << numLines[2] << endl;
+	file << "POINTS " << numLines[0]*numLines[1]*numLines[2] << " float" << endl;
+	for (unsigned int k=0;k<numLines[2];++k)
+		for (unsigned int j=0;j<numLines[1];++j)
+			for (unsigned int i=0;i<numLines[0];++i)
+			{
+				file << setprecision(precision) << discLines[0][i] * cos(discLines[1][j]) << " " << discLines[0][i] * sin(discLines[1][j])  << " " << discLines[2][k] << endl;
+			}
+	file << endl;
+	file << endl << endl;
+	file << "POINT_DATA " << numLines[0]*numLines[1]*numLines[2] << endl;
+}
+
+
 void ProcessFields::WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, unsigned int precision)
 {
 	file << "VECTORS " << name << " float " << endl;
@@ -302,16 +336,16 @@ void ProcessFields::WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT
 }
 
 
-bool ProcessFields::DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const* const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info)
+bool ProcessFields::DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const* const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT)
 {
-	WriteVTKHeader(file, discLines, numLines, precision, header_info);
+	WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT);
 	WriteVTKVectorArray(file, name, array, numLines, precision);
 	return true;
 }
 
-bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info)
+bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT)
 {
-	WriteVTKHeader(file, discLines, numLines, precision, header_info);
+	WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT);
 	for (unsigned int n=0;n<numFields;++n)
 	{
 		WriteVTKVectorArray(file, names[n], array[n], numLines, precision);
@@ -341,16 +375,16 @@ void ProcessFields::WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT
 	}
 }
 
-bool ProcessFields::DumpScalarArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info)
+bool ProcessFields::DumpScalarArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT)
 {
-	WriteVTKHeader(file, discLines, numLines, precision, header_info);
+	WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT);
 	WriteVTKScalarArray(file, name, array, numLines, precision);
 	return true;
 }
 
-bool ProcessFields::DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info)
+bool ProcessFields::DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT)
 {
-	WriteVTKHeader(file, discLines, numLines, precision, header_info);
+	WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT);
 	for (unsigned int n=0;n<numFields;++n)
 	{
 		WriteVTKScalarArray(file, names[n], array[n], numLines, precision);
diff --git a/FDTD/processfields.h b/FDTD/processfields.h
index 7c9c1a4..a282d9b 100644
--- a/FDTD/processfields.h
+++ b/FDTD/processfields.h
@@ -30,6 +30,7 @@ public:
 	enum FileType { VTK_FILETYPE, HDF5_FILETYPE};
 	enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP};
 	enum DumpMode { NO_INTERPOLATION, NODE_INTERPOLATE, CELL_INTERPOLATE};
+	enum MeshType { CARTESIAN_MESH, CYLINDRICAL_MESH};
 
 	virtual void InitProcess();
 
@@ -54,23 +55,28 @@ public:
 	//! Set dump type: 0 for E-fields, 1 for H-fields, 2 for D-fields, 3 for B-fields, 4 for J-fields, etc...
 	void SetDumpType(DumpType type) {m_DumpType=type;}
 
-	//! Write a vtk header to an already open file
-	static void WriteVTKHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string());
+	//! Write a vtk header to an already open file with given mesh-type
+	static void WriteVTKHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string(), MeshType meshT = CARTESIAN_MESH);
+	//! Write a vtk header to an already open file (cartesian grid)
+	static void WriteVTKCartesianGridHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string());
+	//! Write a vtk header to an already open file (cylindrical grid)
+	static void WriteVTKCylindricalGridHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string());
 	//! Append a vtk vector array to an already open vtk file, write a header first! \sa WriteVTKHeader
 	static void WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, unsigned int precision=12);
 	//! Append a vtk scalar array to an already open vtk file, write a header first! \sa WriteVTKHeader
 	static void WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT const* const* const* array, unsigned int const* numLines, unsigned int precision=12);
 
-	static bool DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const * const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string());
-	static bool DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string());
-	static bool DumpScalarArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string());
-	static bool DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string());
+	static bool DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const * const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string(), MeshType meshT = CARTESIAN_MESH);
+	static bool DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string(), MeshType meshT = CARTESIAN_MESH);
+	static bool DumpScalarArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string(), MeshType meshT = CARTESIAN_MESH);
+	static bool DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision=12, string header_info = string(), MeshType meshT = CARTESIAN_MESH);
 
 	static bool DumpVectorArray2HDF5(string filename, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, float time=0);
 
 	double CalcTotalEnergy() const;
 
 	void SetFileType(FileType fileType) {m_fileType=fileType;}
+	void SetMeshType(MeshType meshType) {m_Mesh_Type=meshType;}
 
 	static string GetFieldNameByType(DumpType type);
 	static string GetInterpolationNameByType(DumpMode mode);
@@ -79,6 +85,7 @@ public:
 protected:
 	DumpMode m_DumpMode;
 	DumpType m_DumpType;
+	MeshType m_Mesh_Type;
 	string filePattern;
 	string m_fileName;
 	FileType m_fileType;
diff --git a/FDTD/processfields_td.cpp b/FDTD/processfields_td.cpp
index 22f2099..c0cb0cb 100644
--- a/FDTD/processfields_td.cpp
+++ b/FDTD/processfields_td.cpp
@@ -77,7 +77,7 @@ void ProcessFieldsTD::DumpNodeInterpol(string filename)
 		{
 			ofstream file(filename.c_str());
 			if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;};
-			DumpVectorArray2VTK(file,string("H-Field"),H_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode));
+			DumpVectorArray2VTK(file,string("H-Field"),H_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type);
 			file.close();
 		}
 		else if (m_fileType==HDF5_FILETYPE)
@@ -133,7 +133,7 @@ void ProcessFieldsTD::DumpNodeInterpol(string filename)
 		{
 			ofstream file(filename.c_str());
 			if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;};
-			DumpVectorArray2VTK(file,string("E-Field"),E_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode));
+			DumpVectorArray2VTK(file,string("E-Field"),E_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type);
 			file.close();
 		}
 		else if (m_fileType==HDF5_FILETYPE)
@@ -197,7 +197,7 @@ void ProcessFieldsTD::DumpCellInterpol(string filename)
 		{
 			ofstream file(filename.c_str());
 			if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;};
-			DumpVectorArray2VTK(file,string("E-Field"),E_T,discDLines,numDLines,m_precision,GetInterpolationNameByType(m_DumpMode));
+			DumpVectorArray2VTK(file,string("E-Field"),E_T,discDLines,numDLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type);
 			file.close();
 		}
 		else if (m_fileType==HDF5_FILETYPE)
@@ -253,7 +253,7 @@ void ProcessFieldsTD::DumpCellInterpol(string filename)
 		{
 			ofstream file(filename.c_str());
 			if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;};
-			DumpVectorArray2VTK(file,string("H-Field"),H_T,discDLines,numDLines,m_precision,GetInterpolationNameByType(m_DumpMode));
+			DumpVectorArray2VTK(file,string("H-Field"),H_T,discDLines,numDLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type);
 			file.close();
 		}
 		else if (m_fileType==HDF5_FILETYPE)
@@ -303,7 +303,7 @@ void ProcessFieldsTD::DumpNoInterpol(string filename)
 		{
 			ofstream file(filename.c_str());
 			if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;};
-			DumpVectorArray2VTK(file,string("E-Field"),E_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode));
+			DumpVectorArray2VTK(file,string("E-Field"),E_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type);
 			file.close();
 		}
 		else if (m_fileType==HDF5_FILETYPE)
@@ -349,7 +349,7 @@ void ProcessFieldsTD::DumpNoInterpol(string filename)
 		{
 			ofstream file(filename.c_str());
 			if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;};
-			DumpVectorArray2VTK(file,string("H-Field"),H_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode));
+			DumpVectorArray2VTK(file,string("H-Field"),H_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type);
 			file.close();
 		}
 		else if (m_fileType==HDF5_FILETYPE)
diff --git a/openems.cpp b/openems.cpp
index d117545..2f4cc63 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -353,6 +353,8 @@ int openEMS::SetupFDTD(const char* file)
 				ProcTD->SetDumpType((ProcessFields::DumpType)db->GetDumpType());
 				ProcTD->SetDumpMode((ProcessFields::DumpMode)db->GetDumpMode());
 				ProcTD->SetFileType((ProcessFields::FileType)db->GetFileType());
+				if (CylinderCoords)
+					ProcTD->SetMeshType(ProcessFields::CYLINDRICAL_MESH);
 				for (int n=0;n<3;++n)
 					ProcTD->SetSubSampling(db->GetSubSampling(n),n);
 				ProcTD->SetFilePattern(db->GetName());