first hdf5-file dump implementation

FDTD/processfields.cpp

@@ -17,10 +17,14 @@
 
 #include "processfields.h"
 
+#include "H5Cpp.h"
+
 ProcessFields::ProcessFields(Operator* op, Engine* eng) : Processing(op, eng)
 {
-	DumpMode=0;
+	m_DumpMode = NO_INTERPOLATION;
-	DumpType = 0;
+	m_DumpType = E_FIELD_DUMP;
+	// vtk-file is default
+	m_fileType = VTK_FILETYPE;
 //	SetSubSampling(1);
 
 	for (int n=0;n<3;++n)
@@ -42,9 +46,68 @@ ProcessFields::~ProcessFields()
 	}
 }
 
+void ProcessFields::InitProcess()
+{
+	string names[] = {"x","y","z"};
+	if (m_fileType==HDF5_FILETYPE)
+	{
+		unsigned int* NrLines;
+		double** Lines;
+
+		if (m_DumpMode==CELL_INTERPOLATE)
+		{
+			NrLines = numDLines;
+			Lines = discDLines;
+		}
+		else if (m_DumpMode==NO_INTERPOLATION)
+		{
+			NrLines = numLines;
+			Lines = discLines;
+		}
+		else
+			return;
+
+		m_fileName+= ".h5";
+		H5::H5File* file = new H5::H5File( m_fileName , H5F_ACC_TRUNC );
+
+		H5::Group* group = new H5::Group( file->createGroup( "/Mesh" ));
+		for (int n=0;n<3;++n)
+		{
+			hsize_t dimsf[1];              // dataset dimensions
+			dimsf[0] = NrLines[n];
+			H5::DataSpace dataspace( 1, dimsf );
+			H5::FloatType datatype( H5::PredType::NATIVE_FLOAT );
+			H5::DataSet dataset = group->createDataSet( names[n].c_str(), datatype, dataspace );
+			//convert to float...
+			float* array = new float[NrLines[n]];
+			for (int i=0;i<NrLines[n];++i)
+				array[i] = Lines[n][i];
+			//write to dataset
+			dataset.write( array, H5::PredType::NATIVE_FLOAT );
+		}
+		delete group;
+
+		group = new H5::Group( file->createGroup( "/FieldData" ));
+		delete group;
+		delete file;
+	}
+}
+
+string ProcessFields::GetFieldNameByType(DumpType type)
+{
+	switch (type)
+	{
+		case E_FIELD_DUMP:
+			return "E-Field";
+		case H_FIELD_DUMP:
+			return "H-Field";
+	}
+	return "unknown field";
+}
+
 void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
 {
-	if (DumpMode==0)
+	if (m_DumpMode==NO_INTERPOLATION)
 	{
 		if (Op->SnapToMesh(dstart,start)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl;
 		if (Op->SnapToMesh(dstop,stop)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
@@ -69,7 +132,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
 			}
 		}
 	}
-	else if (DumpMode==2)
+	else if (m_DumpMode==CELL_INTERPOLATE)
 	{
 		if (Op->SnapToMesh(dstart,start,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl;
 		if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
@@ -86,7 +149,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
 			}
 			++stop[n];
 			numDLines[n]=stop[n]-start[n];
-	//		cerr << " number of lines " << numDLines[n] << endl;
+//			cerr << " number of lines " << numDLines[n] << endl;
 			delete[] discDLines[n];
 			discDLines[n] = new double[numDLines[n]];
 			for (unsigned int i=0;i<numDLines[n];++i)
@@ -238,4 +301,47 @@ bool ProcessFields::DumpMultiScalarArray2VTK(ofstream &file, string names[], FDT
 	return true;
 }
 
+bool ProcessFields::DumpVectorArray2HDF5(string filename, string name, FDTD_FLOAT**** array, unsigned int* numLines)
+{
+	const H5std_string FILE_NAME(filename);
+	const H5std_string DATASET_NAME( name );
+
+	H5::H5File file( FILE_NAME, H5F_ACC_RDWR );
+
+	H5::Group group( file.openGroup( "/FieldData" ));
+
+	hsize_t dimsf[4];              // dataset dimensions
+
+	dimsf[0] = 3;
+	dimsf[1] = numLines[2];
+	dimsf[2] = numLines[1];
+	dimsf[3] = numLines[0];
+
+	H5::DataSpace dataspace( 4, dimsf );
+
+	H5::FloatType datatype( H5::PredType::NATIVE_FLOAT );
+//	datatype.setOrder( H5T_ORDER_LE );
+	H5::DataSet dataset = group.createDataSet( DATASET_NAME, datatype, dataspace );
+
+	// I have not the slightest idea why this array-copy action is necessary...  but it's the only way hdf5 does what it is supposed to do anyway!!
+	// at least it is save in case FDTD_FLOAT was defined as double...
+	// why does hdf5 write the dimensions backwards??? or matlab???
+	float hdf5array[3][numLines[2]][numLines[1]][numLines[0]];
+	for (int n=0;n<3;++n)
+	{
+		for (unsigned int i=0;i<numLines[0];++i)
+		{
+			for (unsigned int j=0;j<numLines[1];++j)
+			{
+				for (unsigned int k=0;k<numLines[2];++k)
+				{
+					hdf5array[n][k][j][i] = array[n][i][j][k];
+				}
+			}
+		}
+	}
+	dataset.write( hdf5array, H5::PredType::NATIVE_FLOAT );
+	return true;
+}
+
 

FDTD/processfields.h

@@ -27,37 +27,53 @@ public:
 	ProcessFields(Operator* op, Engine* eng);
 	virtual ~ProcessFields();
 
+	enum FileType { VTK_FILETYPE, HDF5_FILETYPE};
+	enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP};
+	enum DumpMode { NO_INTERPOLATION, NODE_INTERPOLATE, CELL_INTERPOLATE};
+
+	virtual void InitProcess();
+
 	virtual void DefineStartStopCoord(double* dstart, double* dstop);
 
 //	virtual void SetSubSampling(unsigned int subSampleRate, int dir=-1);
 
-	//! Used file pattern e.g. pattern="tmp/efield_" --> "tmp/efield_000045.vtk" for timestep 45 or "tmp/efield_2.40000e9.vtk" for 2.4GHz E-field dump.
+	//! Used file pattern e.g. pattern="tmp/efield_" --> "tmp/efield_000045.vtk" for timestep 45 or "tmp/efield_2.40000e9.vtk" for 2.4GHz E-field dump. (VTK FileType only) \sa SetFileType
 	void SetFilePattern(string fp) {filePattern=fp;}
 
+	//! Set the filename for a hdf5 data group file (HDF5 FileType only) \sa SetFileType
+	void SetFileName(string fn) {m_fileName=fn;}
+
 	//! Define the Dump-Mode
-	void SetDumpMode(int mode) {DumpMode=mode;}
+	void SetDumpMode(DumpMode mode) {m_DumpMode=mode;}
 	//! This methode will dump all fields in the center of a main cell (dual-node) using 4 E-field and 2 H-fields per direction. (default)
-	void SetDumpMode2Cell() {DumpMode=2;}
+	void SetDumpMode2Cell() {m_DumpMode=CELL_INTERPOLATE;}
 
 	//! 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(int type) {DumpType=type;}
+	void SetDumpType(DumpType type) {m_DumpType=type;}
 
 	static bool DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT**** array, double** discLines, unsigned int* numLines);
 	static bool DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT**** array[], unsigned int numFields, double** discLines, unsigned int* numLines);
 	static bool DumpScalarArray2VTK(ofstream &file, string name, FDTD_FLOAT*** array, double** discLines, unsigned int* numLines);
 	static bool DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT*** array[], unsigned int numFields, double** discLines, unsigned int* numLines);
 
+	static bool DumpVectorArray2HDF5(string filename, string name, FDTD_FLOAT**** array, unsigned int* numLines);
+
 	double CalcTotalEnergy();
 
+	void SetFileType(FileType fileType) {m_fileType=fileType;}
+
 //	virtual void Process();
 protected:
 	static void WriteVTKHeader(ofstream &file, double** discLines, unsigned int* numLines);
 	static void WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT**** array, unsigned int* numLines);
 	static void WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT*** array, unsigned int* numLines);
 
-	int DumpMode;
+	static string GetFieldNameByType(DumpType type);
-	int DumpType;
+	DumpMode m_DumpMode;
+	DumpType m_DumpType;
 	string filePattern;
+	string m_fileName;
+	FileType m_fileType;
 
 //	unsigned int subSample[3];
 

FDTD/processfields_td.cpp

@@ -29,12 +29,12 @@ ProcessFieldsTD::~ProcessFieldsTD()
 {
 }
 
-void ProcessFieldsTD::DumpCellInterpol(ofstream &file)
+void ProcessFieldsTD::DumpCellInterpol(string filename)
 {
 	FDTD_FLOAT**** volt = Eng->GetVoltages();
 	FDTD_FLOAT**** curr = Eng->GetCurrents();
 
-	if (DumpType==0)
+	if (m_DumpType==E_FIELD_DUMP)
 	{
 		//create array
 		FDTD_FLOAT**** E_T = Create_N_3DArray(numDLines);
@@ -66,12 +66,27 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file)
 				}
 			}
 		}
-		DumpVectorArray2VTK(file,string("E-Field"),E_T,discDLines,numDLines);
+
+		if (m_fileType==VTK_FILETYPE)
+		{
+			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);
+			file.close();
+		}
+		else if (m_fileType==HDF5_FILETYPE)
+		{
+			stringstream ss;
+			ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps();
+			DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),E_T,numDLines);
+		}
+		else
+			cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl;
 		Delete_N_3DArray(E_T,numDLines);
 		E_T = NULL;
 	}
 
-	if (DumpType==1)
+	if (m_DumpType==1)
 	{
 		//create array
 		FDTD_FLOAT**** H_T = Create_N_3DArray(numDLines);
@@ -104,20 +119,34 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file)
 				}
 			}
 		}
-		DumpVectorArray2VTK(file,string("H-Field"),H_T,discDLines,numDLines);
+		if (m_fileType==VTK_FILETYPE)
+		{
+			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);
+			file.close();
+		}
+		else if (m_fileType==HDF5_FILETYPE)
+		{
+			stringstream ss;
+			ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps();
+			DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),H_T,numDLines);
+		}
+		else
+			cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl;
 		Delete_N_3DArray(H_T,numDLines);
 		H_T = NULL;
 	}
 }
 
-void ProcessFieldsTD::DumpNoInterpol(ofstream &file)
+void ProcessFieldsTD::DumpNoInterpol(string filename)
 {
 	FDTD_FLOAT**** volt = Eng->GetVoltages();
 	FDTD_FLOAT**** curr = Eng->GetCurrents();
 
 	unsigned int pos[3];
 	double delta[3];
-	if (DumpType==0)
+	if (m_DumpType==E_FIELD_DUMP)
 	{
 		//create array
 		FDTD_FLOAT**** E_T = Create_N_3DArray(numLines);
@@ -136,12 +165,27 @@ void ProcessFieldsTD::DumpNoInterpol(ofstream &file)
 				}
 			}
 		}
-		DumpVectorArray2VTK(file,string("E-Field"),E_T,discLines,numLines);
+		if (m_fileType==VTK_FILETYPE)
+		{
+			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);
+			file.close();
+		}
+		else if (m_fileType==HDF5_FILETYPE)
+		{
+			stringstream ss;
+			ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps();
+			DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),E_T,numLines);
+		}
+		else
+			cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl;
+
 		Delete_N_3DArray(E_T,numLines);
 		E_T = NULL;
 	}
 
-	if (DumpType==1)
+	if (m_DumpType==H_FIELD_DUMP)
 	{
 		//create array
 		FDTD_FLOAT**** H_T = Create_N_3DArray(numLines);
@@ -161,7 +205,22 @@ void ProcessFieldsTD::DumpNoInterpol(ofstream &file)
 				}
 			}
 		}
-		DumpVectorArray2VTK(file,string("H-Field"),H_T,discLines,numLines);
+		if (m_fileType==VTK_FILETYPE)
+		{
+			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);
+			file.close();
+		}
+		else if (m_fileType==HDF5_FILETYPE)
+		{
+			stringstream ss;
+			ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps();
+			DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),H_T,numLines);
+		}
+		else
+			cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl;
+
 		Delete_N_3DArray(H_T,numLines);
 		H_T = NULL;
 	}
@@ -173,16 +232,22 @@ int ProcessFieldsTD::Process()
 	if (filePattern.empty()) return -1;
 	if (CheckTimestep()==false) return GetNextInterval();
 	stringstream ss;
-	ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps();
+	ss << filePattern << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps();
 
-	string filename = filePattern + ss.str() + ".vtk";
+	if (m_fileType==VTK_FILETYPE)
-	ofstream file(filename.c_str());
+	{
-	if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl; return GetNextInterval();};
+		ss << ".vtk";
-
+		if (m_DumpMode==NO_INTERPOLATION)
-	if (DumpMode==0)
+			DumpNoInterpol(ss.str());
-		DumpNoInterpol(file);
+		if (m_DumpMode==CELL_INTERPOLATE)
-	if (DumpMode==2)
+			DumpCellInterpol(ss.str());
-		DumpCellInterpol(file);
+	}
-	file.close();
+	else if (m_fileType==HDF5_FILETYPE)
+	{
+		if (m_DumpMode==NO_INTERPOLATION)
+			DumpNoInterpol(m_fileName);
+		if (m_DumpMode==CELL_INTERPOLATE)
+			DumpCellInterpol(m_fileName);
+	}
 	return GetNextInterval();
 }

FDTD/processfields_td.h

@@ -34,8 +34,8 @@ public:
 protected:
 	int pad_length;
 
-	void DumpNoInterpol(ofstream &file);
+	void DumpNoInterpol(string filename);
-	void DumpCellInterpol(ofstream &file);
+	void DumpCellInterpol(string filename);
 };
 
 #endif // PROCESSFIELDS_TD_H

matlab/examples/PlaneWave.m

@@ -53,13 +53,13 @@ CSX = AddExcitation(CSX,'excite',0,[0 1 0]);
 CSX = AddBox(CSX,'excite',0 ,start,stop);
  
 %dump
-CSX = AddDump(CSX,'Et_',0,0);
+CSX = AddDump(CSX,'Et',0,0,1);
 start = [mesh.x(1) , 0 , mesh.z(1)];
 stop = [mesh.x(end) , 0 , mesh.z(end)];
-CSX = AddBox(CSX,'Et_',0 , start,stop);
+CSX = AddBox(CSX,'Et',0 , start,stop);
 
-CSX = AddDump(CSX,'Ht_',1,0);
+CSX = AddDump(CSX,'Ht',1,0,1);
-CSX = AddBox(CSX,'Ht_',0,start,stop);
+CSX = AddBox(CSX,'Ht',0,start,stop);
 
 %Write openEMS compatoble xml-file
 WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);

openEMS.pro

@@ -15,10 +15,12 @@ LIBS += -L../CSXCAD \
     -lfparser \
     -L../tinyxml \
 	-ltinyxml \
-	-lboost_thread
+	-lboost_thread \
+	-lhdf5_cpp
 QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../CSXCAD\'
 QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../fparser\'
 QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../tinyxml\'
+
 SOURCES += main.cpp \
     tools/ErrorMsg.cpp \
     tools/AdrOp.cpp \

openems.cpp

@@ -346,10 +346,13 @@ int openEMS::SetupFDTD(const char* file)
 			CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox();
 			if (db)
 			{
-				ProcTD->SetDumpType(db->GetDumpType());
+				ProcTD->SetDumpType((ProcessFields::DumpType)db->GetDumpType());
-				ProcTD->SetDumpMode(db->GetDumpMode());
+				ProcTD->SetDumpMode((ProcessFields::DumpMode)db->GetDumpMode());
+				ProcTD->SetFileType((ProcessFields::FileType)db->GetFileType());
 				ProcTD->SetFilePattern(db->GetName());
+				ProcTD->SetFileName(db->GetName());
 				ProcTD->DefineStartStopCoord(start,stop);
+				ProcTD->InitProcess();
 				PA->AddProcessing(ProcTD);
 			}
 			else