zcy
/
openEMS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

reorganized vtk writer and new hdf5 file writer & reader

pull/1/head
Thorsten Liebig 2012-02-02 10:20:49 +01:00
parent b788b0bf23
commit 807786b2d2
15 changed files with 1078 additions and 496 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

244
Common/processfields.cpp
View File

@ -18,7 +18,8 @@
#include <iomanip> #include <iomanip>
#include <H5Cpp.h> #include <H5Cpp.h>
#include "tools/global.h" #include "tools/global.h"
#include "tools/vtk_file_io.h" #include "tools/vtk_file_writer.h"
#include "tools/hdf5_file_writer.h"
#include "processfields.h" #include "processfields.h"
#include "FDTD/engine_interface_fdtd.h" #include "FDTD/engine_interface_fdtd.h"
@ -28,7 +29,8 @@ ProcessFields::ProcessFields(Engine_Interface_Base* eng_if) : Processing(eng_if)
// vtk-file is default // vtk-file is default
m_fileType = VTK_FILETYPE; m_fileType = VTK_FILETYPE;
m_SampleType = NONE; m_SampleType = NONE;
m_Dump_File = NULL; m_Vtk_Dump_File = NULL;
m_HDF5_Dump_File = NULL;
SetPrecision(6); SetPrecision(6);
m_dualTime = false; m_dualTime = false;
@ -47,8 +49,8 @@ ProcessFields::ProcessFields(Engine_Interface_Base* eng_if) : Processing(eng_if)
ProcessFields::~ProcessFields() ProcessFields::~ProcessFields()
{ {
delete m_Dump_File; delete m_Vtk_Dump_File;
m_Dump_File = NULL; m_Vtk_Dump_File = NULL;
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
{ {
delete[] posLines[n]; delete[] posLines[n];
@ -84,16 +86,28 @@ void ProcessFields::InitProcess()
if (m_fileType==VTK_FILETYPE) if (m_fileType==VTK_FILETYPE)
{ {
delete m_Dump_File; delete m_Vtk_Dump_File;
m_Dump_File = new VTK_File_IO(m_filename,(int)m_Mesh_Type); m_Vtk_Dump_File = new VTK_File_Writer(m_filename,(int)m_Mesh_Type);
#ifdef OUTPUT_IN_DRAWINGUNITS #ifdef OUTPUT_IN_DRAWINGUNITS
double discScaling = 1; double discScaling = 1;
#else #else
double discScaling = Op->GetGridDelta(); double discScaling = Op->GetGridDelta();
#endif #endif
m_Dump_File->SetMeshLines(discLines,numLines,discScaling); m_Vtk_Dump_File->SetMeshLines(discLines,numLines,discScaling);
m_Dump_File->SetNativeDump(g_settings.NativeFieldDumps()); m_Vtk_Dump_File->SetNativeDump(g_settings.NativeFieldDumps());
}
if (m_fileType==HDF5_FILETYPE)
{
delete m_HDF5_Dump_File;
m_HDF5_Dump_File = new HDF5_File_Writer(m_filename+".h5");
#ifdef OUTPUT_IN_DRAWINGUNITS
double discScaling = 1;
#else
double discScaling = Op->GetGridDelta();
#endif
m_HDF5_Dump_File->WriteRectMesh(numLines,discLines,(int)m_Mesh_Type,discScaling);
} }
} }
@ -216,220 +230,6 @@ void ProcessFields::CalcMeshPos()
} }
} }
bool ProcessFields::WriteMesh2HDF5(string filename, string groupName, unsigned int const* numLines, double const* const* discLines, MeshType meshT, double discLines_scaling)
{
H5::H5File file( filename, H5F_ACC_RDWR );
H5::Group hdf_group( file.openGroup( groupName ));
string names[] = {"x","y","z"};
if (meshT==CYLINDRICAL_MESH)
{
names[0]="rho";
names[1]="alpha";
}
H5::Group* group = new H5::Group( hdf_group.createGroup( "/Mesh" ));
for (int n=0; n<3; ++n)
{
hsize_t dimsf[1]; // dataset dimensions
dimsf[0] = numLines[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[numLines[n]];
for (unsigned int i=0; i<numLines[n]; ++i)
{
#ifdef OUTPUT_IN_DRAWINGUNITS
array[i] = Lines[n][i];
#else
if ((meshT==CYLINDRICAL_MESH) && (n==1)) //check for alpha-direction
array[i] = discLines[n][i];
else
array[i] = discLines[n][i] * discLines_scaling;
#endif
}
//write to dataset
dataset.write( array, H5::PredType::NATIVE_FLOAT );
delete[] array;
}
delete group;
return true;
}
bool ProcessFields::DumpVectorArray2HDF5(string filename, string groupName, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, float time)
{
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( groupName ));
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 );
hsize_t t_dimsf[] = {1};
H5::DataSpace t_dataspace( 1, t_dimsf );
H5::Attribute attr = dataset.createAttribute("time",H5::PredType::NATIVE_FLOAT,t_dataspace);
attr.write( H5::PredType::NATIVE_FLOAT , &time);
// 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???
unsigned long pos = 0;
float *hdf5array = new float[3*numLines[0]*numLines[1]*numLines[2]];
for (int n=0; n<3; ++n)
{
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)
{
hdf5array[pos++] = array[n][i][j][k];
}
}
}
}
dataset.write( hdf5array, H5::PredType::NATIVE_FLOAT );
delete[] hdf5array;
return true;
}
bool ProcessFields:: DumpScalarArray2HDF5(string filename, string groupName, string name, FDTD_FLOAT const* const* const* array, unsigned int const* numLines, string attr_name, float attr_value)
{
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( groupName ));
hsize_t dimsf[3]; // dataset dimensions
dimsf[0] = numLines[2];
dimsf[1] = numLines[1];
dimsf[2] = numLines[0];
H5::DataSpace dataspace( 3, dimsf );
H5::FloatType datatype( H5::PredType::NATIVE_FLOAT );
// datatype.setOrder( H5T_ORDER_LE );
H5::DataSet dataset = group.createDataSet( DATASET_NAME, datatype, dataspace );
if (!attr_name.empty())
{
hsize_t t_dimsf[] = {1};
H5::DataSpace t_dataspace( 1, t_dimsf );
H5::Attribute attr = dataset.createAttribute(attr_name,H5::PredType::NATIVE_FLOAT,t_dataspace);
attr.write( H5::PredType::NATIVE_FLOAT , &attr_value);
}
// 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???
unsigned long pos = 0;
float *hdf5array = new float[numLines[0]*numLines[1]*numLines[2]];
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)
{
hdf5array[pos++] = array[i][j][k];
}
}
}
dataset.write( hdf5array, H5::PredType::NATIVE_FLOAT );
delete[] hdf5array;
return true;
}
bool ProcessFields::DumpVectorArray2HDF5(string filename, string groupName, string name, std::complex<float> const* const* const* const* array, unsigned int const* numLines, float weight, float frequency)
{
const H5std_string FILE_NAME(filename);
const H5std_string DATASET_NAME_RE( name + "_real");
const H5std_string DATASET_NAME_IM( name + "_imag");
H5::H5File file( FILE_NAME, H5F_ACC_RDWR );
H5::Group group( file.openGroup( groupName ));
hsize_t t_dimsf[] = {1};
H5::DataSpace t_dataspace( 1, t_dimsf );
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 );
//create and write real part
H5::DataSet dataset = group.createDataSet( DATASET_NAME_RE, datatype, dataspace );
H5::Attribute attr = dataset.createAttribute("frequency",H5::PredType::NATIVE_FLOAT,t_dataspace);
attr.write( H5::PredType::NATIVE_FLOAT , &frequency);
// 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???
unsigned long pos = 0;
float *hdf5array = new float[3*numLines[0]*numLines[1]*numLines[2]];
for (int n=0; n<3; ++n)
{
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)
{
hdf5array[pos++] = array[n][i][j][k].real() * weight;
}
}
}
}
dataset.write( hdf5array, H5::PredType::NATIVE_FLOAT );
//create and write imaginary part
dataset = group.createDataSet( DATASET_NAME_IM, datatype, dataspace );
attr = dataset.createAttribute("frequency",H5::PredType::NATIVE_FLOAT,t_dataspace);
attr.write( H5::PredType::NATIVE_FLOAT , &frequency);
// 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???
pos=0;
for (int n=0; n<3; ++n)
{
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)
{
hdf5array[pos++] = array[n][i][j][k].imag() * weight;
}
}
}
}
dataset.write( hdf5array, H5::PredType::NATIVE_FLOAT );
delete[] hdf5array;
return true;
}
FDTD_FLOAT**** ProcessFields::CalcField() FDTD_FLOAT**** ProcessFields::CalcField()
{ {
unsigned int pos[3]; unsigned int pos[3];

18
Common/processfields.h
View File

@ -23,7 +23,8 @@
#define __VTK_DATA_TYPE__ "double" #define __VTK_DATA_TYPE__ "double"
class Base_File_IO; class VTK_File_Writer;
class HDF5_File_Writer;
class ProcessFields : public Processing class ProcessFields : public Processing
{ {
@ -67,18 +68,6 @@ 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... //! 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;} void SetDumpType(DumpType type) {m_DumpType=type;}
//! Write a mesh information to the given hdf5-group
static bool WriteMesh2HDF5(string filename, string groupName, unsigned int const* numLines, double const* const* discLines, MeshType meshT = CARTESIAN_MESH, double discLines_scaling = 1);
//! Dump a time-domain vector dump to an HDF5 file
static bool DumpVectorArray2HDF5(string filename, string groupName, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, float time=0);
//! Dump a scalar field to an HDF5 file
static bool DumpScalarArray2HDF5(string filename, string groupName, string name, FDTD_FLOAT const* const* const* array, unsigned int const* numLines, string attr_name=string(), float attr_value=0);
//! Dump a frequency-domain complex-vector dump to an HDF5 file
static bool DumpVectorArray2HDF5(string filename, string groupName, string name, std::complex<float> const* const* const* const* array, unsigned int const* numLines, float weight, float frequency);
double CalcTotalEnergyEstimate() const; double CalcTotalEnergyEstimate() const;
void SetFileType(FileType fileType) {m_fileType=fileType;} void SetFileType(FileType fileType) {m_fileType=fileType;}
@ -89,7 +78,8 @@ protected:
DumpType m_DumpType; DumpType m_DumpType;
FileType m_fileType; FileType m_fileType;
Base_File_IO* m_Dump_File; VTK_File_Writer* m_Vtk_Dump_File;
HDF5_File_Writer* m_HDF5_Dump_File;
enum SampleType {NONE, SUBSAMPLE, OPT_RESOLUTION} m_SampleType; enum SampleType {NONE, SUBSAMPLE, OPT_RESOLUTION} m_SampleType;
virtual void CalcMeshPos(); virtual void CalcMeshPos();

69
Common/processfields_fd.cpp
View File

@ -17,7 +17,8 @@
#include "processfields_fd.h" #include "processfields_fd.h"
#include "Common/operator_base.h" #include "Common/operator_base.h"
#include "tools/vtk_file_io.h" #include "tools/vtk_file_writer.h"
#include "tools/hdf5_file_writer.h"
#include <H5Cpp.h> #include <H5Cpp.h>
#include <iomanip> #include <iomanip>
#include <sstream> #include <sstream>
@ -50,35 +51,14 @@ void ProcessFieldsFD::InitProcess()
//setup the hdf5 file //setup the hdf5 file
ProcessFields::InitProcess(); ProcessFields::InitProcess();
if (m_Dump_File) if (m_Vtk_Dump_File)
m_Dump_File->SetHeader(string("openEMS FD Field Dump -- Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString()); m_Vtk_Dump_File->SetHeader(string("openEMS FD Field Dump -- Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString());
if (m_fileType==HDF5_FILETYPE) if (m_HDF5_Dump_File)
{ {
//create hdf5 file & necessary groups m_HDF5_Dump_File->SetCurrentGroup("/FieldData/FD");
m_filename+= ".h5"; int numFreq = m_FD_Samples.size();
H5::H5File* file = new H5::H5File( m_filename , H5F_ACC_TRUNC ); m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD","Number_of_Frequencies",&numFreq,1,H5T_NATIVE_INT);
H5::Group* group = new H5::Group( file->createGroup( "/FieldData" ));
delete group;
group = new H5::Group( file->createGroup( "/FieldData/FD" ));
//set number of frequencies
hsize_t t_dimsf[] = {1};
H5::DataSpace t_dataspace( 1, t_dimsf );
H5::Attribute attr = group->createAttribute("Number_of_Frequencies",H5::PredType::NATIVE_INT,t_dataspace);
int count = m_FD_Samples.size();
attr.write( H5::PredType::NATIVE_INT , &count);
delete group;
delete file;
//write mesh information in main root-group
#ifdef OUTPUT_IN_DRAWINGUNITS
double discScaling = 1;
#else
double discScaling = Op->GetGridDelta();
#endif
ProcessFields::WriteMesh2HDF5(m_filename,"/",numLines,discLines,m_Mesh_Type, discScaling);
} }
//create data structures... //create data structures...
@ -164,10 +144,10 @@ void ProcessFieldsFD::DumpFDData()
stringstream ss; stringstream ss;
ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_p=" << std::setw( 3 ) << std::setfill( '0' ) <<(int)(angle * 180 / M_PI); ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_p=" << std::setw( 3 ) << std::setfill( '0' ) <<(int)(angle * 180 / M_PI);
m_Dump_File->SetFilename(ss.str()); m_Vtk_Dump_File->SetFilename(ss.str());
m_Dump_File->ClearAllFields(); m_Vtk_Dump_File->ClearAllFields();
m_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field,numLines); m_Vtk_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field);
if (m_Dump_File->Write()==false) if (m_Vtk_Dump_File->Write()==false)
cerr << "ProcessFieldsFD::Process: can't dump to file... abort! " << endl; cerr << "ProcessFieldsFD::Process: can't dump to file... abort! " << endl;
} }
@ -187,10 +167,10 @@ void ProcessFieldsFD::DumpFDData()
} }
stringstream ss; stringstream ss;
ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_abs"; ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_abs";
m_Dump_File->SetFilename(ss.str()); m_Vtk_Dump_File->SetFilename(ss.str());
m_Dump_File->ClearAllFields(); m_Vtk_Dump_File->ClearAllFields();
m_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field,numLines); m_Vtk_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field);
if (m_Dump_File->Write()==false) if (m_Vtk_Dump_File->Write()==false)
cerr << "ProcessFieldsFD::Process: can't dump to file... abort! " << endl; cerr << "ProcessFieldsFD::Process: can't dump to file... abort! " << endl;
} }
@ -210,10 +190,10 @@ void ProcessFieldsFD::DumpFDData()
} }
stringstream ss; stringstream ss;
ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_arg"; ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_arg";
m_Dump_File->SetFilename(ss.str()); m_Vtk_Dump_File->SetFilename(ss.str());
m_Dump_File->ClearAllFields(); m_Vtk_Dump_File->ClearAllFields();
m_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field,numLines); m_Vtk_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field);
if (m_Dump_File->Write()==false) if (m_Vtk_Dump_File->Write()==false)
cerr << "ProcessFieldsFD::Process: can't dump to file... abort! " << endl; cerr << "ProcessFieldsFD::Process: can't dump to file... abort! " << endl;
} }
} }
@ -227,10 +207,15 @@ void ProcessFieldsFD::DumpFDData()
{ {
stringstream ss; stringstream ss;
ss << "f" << n; ss << "f" << n;
DumpVectorArray2HDF5(m_filename.c_str(), "/FieldData/FD", ss.str(), m_FD_Fields.at(n),numLines,1.0,m_FD_Samples.at(n)); size_t datasize[]={numLines[0],numLines[1],numLines[2]};
if (m_HDF5_Dump_File->WriteVectorField(ss.str(), m_FD_Fields.at(n), datasize)==false)
cerr << "ProcessFieldsFD::Process: can't dump to file...! " << endl;
float freq[1]={m_FD_Samples.at(n)};
if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false)
cerr << "ProcessFieldsFD::Process: can't dump to file...! " << endl;
} }
return; return;
} }
cerr << "ProcessFieldsTD::Process: unknown File-Type" << endl; cerr << "ProcessFieldsFD::Process: unknown File-Type" << endl;
} }

20
Common/processfields_sar.cpp
View File

@ -17,7 +17,8 @@
#include "processfields_sar.h" #include "processfields_sar.h"
#include "operator_base.h" #include "operator_base.h"
#include "tools/vtk_file_io.h" #include "tools/vtk_file_writer.h"
#include "tools/hdf5_file_writer.h"
ProcessFieldsSAR::ProcessFieldsSAR(Engine_Interface_Base* eng_if) : ProcessFieldsFD(eng_if) ProcessFieldsSAR::ProcessFieldsSAR(Engine_Interface_Base* eng_if) : ProcessFieldsFD(eng_if)
{ {
@ -109,17 +110,22 @@ void ProcessFieldsSAR::DumpFDData()
stringstream ss; stringstream ss;
ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n); ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n);
m_Dump_File->SetFilename(ss.str()); m_Vtk_Dump_File->SetFilename(ss.str());
m_Dump_File->ClearAllFields(); m_Vtk_Dump_File->ClearAllFields();
m_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType),SAR,numLines); m_Vtk_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType),SAR);
if (m_Dump_File->Write()==false) if (m_Vtk_Dump_File->Write()==false)
cerr << "ProcessFieldsSAR::Process: can't dump to file... abort! " << endl; cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
} }
else if (m_fileType==HDF5_FILETYPE) else if (m_fileType==HDF5_FILETYPE)
{ {
stringstream ss; stringstream ss;
ss << "f" << n; ss << "f" << n;
DumpScalarArray2HDF5(m_filename.c_str(), "/FieldData/FD", ss.str(), SAR,numLines , "frequency", m_FD_Samples.at(n)); size_t datasize[]={numLines[0],numLines[1],numLines[2]};
if (m_HDF5_Dump_File->WriteScalarField(ss.str(), SAR, datasize)==false)
cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
float freq[1]={m_FD_Samples.at(n)};
if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false)
cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
} }
else else
cerr << "ProcessFieldsSAR::Process: unknown File-Type" << endl; cerr << "ProcessFieldsSAR::Process: unknown File-Type" << endl;

53
Common/processfields_td.cpp
View File

@ -17,7 +17,8 @@
#include "processfields_td.h" #include "processfields_td.h"
#include "Common/operator_base.h" #include "Common/operator_base.h"
#include "tools/vtk_file_io.h" #include "tools/vtk_file_writer.h"
#include "tools/hdf5_file_writer.h"
#include <H5Cpp.h> #include <H5Cpp.h>
#include <iomanip> #include <iomanip>
#include <sstream> #include <sstream>
@ -38,29 +39,11 @@ void ProcessFieldsTD::InitProcess()
ProcessFields::InitProcess(); ProcessFields::InitProcess();
if (m_Dump_File) if (m_Vtk_Dump_File)
m_Dump_File->SetHeader(string("openEMS TD Field Dump -- Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString()); m_Vtk_Dump_File->SetHeader(string("openEMS TD Field Dump -- Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString());
#ifdef OUTPUT_IN_DRAWINGUNITS if (m_HDF5_Dump_File)
double discScaling = 1; m_HDF5_Dump_File->SetCurrentGroup("/FieldData/TD");
#else
double discScaling = Op->GetGridDelta();
#endif
if (m_fileType==HDF5_FILETYPE)
{
//create hdf5 file & necessary groups
m_filename+= ".h5";
H5::H5File* file = new H5::H5File( m_filename, H5F_ACC_TRUNC );
H5::Group* group = new H5::Group( file->createGroup( "/FieldData" ));
delete group;
group = new H5::Group( file->createGroup( "/FieldData/TD" ));
delete group;
delete file;
//write mesh information in main root-group
ProcessFields::WriteMesh2HDF5(m_filename,"/",numLines,discLines,m_Mesh_Type, discScaling);
}
} }
int ProcessFieldsTD::Process() int ProcessFieldsTD::Process()
@ -71,25 +54,37 @@ int ProcessFieldsTD::Process()
string filename = m_filename; string filename = m_filename;
float**** field = CalcField(); float**** field = CalcField();
bool success = true;
if (m_fileType==VTK_FILETYPE) if (m_fileType==VTK_FILETYPE)
{ {
m_Dump_File->SetTimestep(m_Eng_Interface->GetNumberOfTimesteps()); m_Vtk_Dump_File->SetTimestep(m_Eng_Interface->GetNumberOfTimesteps());
m_Dump_File->ClearAllFields(); m_Vtk_Dump_File->ClearAllFields();
m_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field,numLines); m_Vtk_Dump_File->AddVectorField(GetFieldNameByType(m_DumpType),field);
if (m_Dump_File->Write()==false) success &= m_Vtk_Dump_File->Write();
cerr << "ProcessFieldsTD::Process: can't dump to file... abort! " << endl;
} }
else if (m_fileType==HDF5_FILETYPE) else if (m_fileType==HDF5_FILETYPE)
{ {
stringstream ss; stringstream ss;
ss << std::setw( pad_length ) << std::setfill( '0' ) << m_Eng_Interface->GetNumberOfTimesteps(); ss << std::setw( pad_length ) << std::setfill( '0' ) << m_Eng_Interface->GetNumberOfTimesteps();
DumpVectorArray2HDF5(filename.c_str(), "/FieldData/TD", string( ss.str() ), field, numLines, m_Eng_Interface->GetTime(m_dualTime)); size_t datasize[]={numLines[0],numLines[1],numLines[2]};
success &= m_HDF5_Dump_File->WriteVectorField(ss.str(), field, datasize);
float time[1]={m_Eng_Interface->GetTime(m_dualTime)};
success &= m_HDF5_Dump_File->WriteAtrribute("/FieldData/TD/"+ss.str(),"time",time,1);
} }
else else
{
success = false;
cerr << "ProcessFieldsTD::Process: unknown File-Type" << endl; cerr << "ProcessFieldsTD::Process: unknown File-Type" << endl;
}
Delete_N_3DArray<FDTD_FLOAT>(field,numLines); Delete_N_3DArray<FDTD_FLOAT>(field,numLines);
if (success==false)
{
SetEnable(false);
cerr << "ProcessFieldsTD::Process: can't dump to file... disabled! " << endl;
}
return GetNextInterval(); return GetNextInterval();
} }

52
FDTD/operator.cpp
View File

@ -23,7 +23,7 @@
#include "extensions/operator_ext_excitation.h" #include "extensions/operator_ext_excitation.h"
#include "Common/processfields.h" #include "Common/processfields.h"
#include "tools/array_ops.h" #include "tools/array_ops.h"
#include "tools/vtk_file_io.h" #include "tools/vtk_file_writer.h"
#include "fparser.hh" #include "fparser.hh"
Operator* Operator::New() Operator* Operator::New()
@ -493,24 +493,24 @@ void Operator::DumpOperator2File(string filename)
ii_temp[n][pos[0]][pos[1]][pos[2]] = GetII(n,pos); ii_temp[n][pos[0]][pos[1]][pos[2]] = GetII(n,pos);
} }
VTK_File_IO* vtk_io = new VTK_File_IO(filename.c_str(), m_MeshType); VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType);
vtk_io->SetMeshLines(discLines,numLines,discLines_scaling); vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling);
vtk_io->SetHeader("openEMS - Operator dump"); vtk_Writer->SetHeader("openEMS - Operator dump");
vtk_io->SetNativeDump(true); vtk_Writer->SetNativeDump(true);
vtk_io->AddVectorField("vv",vv_temp,numLines); vtk_Writer->AddVectorField("vv",vv_temp);
Delete_N_3DArray(vv_temp,numLines); Delete_N_3DArray(vv_temp,numLines);
vtk_io->AddVectorField("vi",vi_temp,numLines); vtk_Writer->AddVectorField("vi",vi_temp);
Delete_N_3DArray(vi_temp,numLines); Delete_N_3DArray(vi_temp,numLines);
vtk_io->AddVectorField("iv",iv_temp,numLines); vtk_Writer->AddVectorField("iv",iv_temp);
Delete_N_3DArray(iv_temp,numLines); Delete_N_3DArray(iv_temp,numLines);
vtk_io->AddVectorField("ii",ii_temp,numLines); vtk_Writer->AddVectorField("ii",ii_temp);
Delete_N_3DArray(ii_temp,numLines); Delete_N_3DArray(ii_temp,numLines);
vtk_io->AddVectorField("exc",exc,numLines); vtk_Writer->AddVectorField("exc",exc);
Delete_N_3DArray(exc,numLines); Delete_N_3DArray(exc,numLines);
if (vtk_io->Write()==false) if (vtk_Writer->Write()==false)
cerr << "Operator::DumpOperator2File: Error: Can't write file... skipping!" << endl; cerr << "Operator::DumpOperator2File: Error: Can't write file... skipping!" << endl;
cout << " done!" << endl; cout << " done!" << endl;
@ -590,16 +590,16 @@ void Operator::DumpPEC2File( string filename )
scaling = GetGridDelta(); scaling = GetGridDelta();
#endif #endif
VTK_File_IO* vtk_io = new VTK_File_IO(filename.c_str(), m_MeshType); VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType);
vtk_io->SetMeshLines(discLines,numLines,scaling); vtk_Writer->SetMeshLines(discLines,numLines,scaling);
vtk_io->SetHeader("openEMS - PEC dump"); vtk_Writer->SetHeader("openEMS - PEC dump");
vtk_io->SetNativeDump(true); vtk_Writer->SetNativeDump(true);
vtk_io->AddVectorField("PEC",pec,numLines); vtk_Writer->AddVectorField("PEC",pec);
Delete_N_3DArray(pec,numLines); Delete_N_3DArray(pec,numLines);
if (vtk_io->Write()==false) if (vtk_Writer->Write()==false)
cerr << "Operator::DumpPEC2File: Error: Can't write file... skipping!" << endl; cerr << "Operator::DumpPEC2File: Error: Can't write file... skipping!" << endl;
cout << " done!" << endl; cout << " done!" << endl;
@ -640,22 +640,22 @@ void Operator::DumpMaterial2File(string filename)
} }
} }
VTK_File_IO* vtk_io = new VTK_File_IO(filename.c_str(), m_MeshType); VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType);
vtk_io->SetMeshLines(discLines,numLines,discLines_scaling); vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling);
vtk_io->SetHeader("openEMS - material dump"); vtk_Writer->SetHeader("openEMS - material dump");
vtk_io->SetNativeDump(true); vtk_Writer->SetNativeDump(true);
vtk_io->AddVectorField("epsilon",epsilon,numLines); vtk_Writer->AddVectorField("epsilon",epsilon);
Delete_N_3DArray(epsilon,numLines); Delete_N_3DArray(epsilon,numLines);
vtk_io->AddVectorField("mue",mue,numLines); vtk_Writer->AddVectorField("mue",mue);
Delete_N_3DArray(mue,numLines); Delete_N_3DArray(mue,numLines);
vtk_io->AddVectorField("kappa",kappa,numLines); vtk_Writer->AddVectorField("kappa",kappa);
Delete_N_3DArray(kappa,numLines); Delete_N_3DArray(kappa,numLines);
vtk_io->AddVectorField("sigma",sigma,numLines); vtk_Writer->AddVectorField("sigma",sigma);
Delete_N_3DArray(sigma,numLines); Delete_N_3DArray(sigma,numLines);
if (vtk_io->Write()==false) if (vtk_Writer->Write()==false)
cerr << "Operator::DumpMaterial2File: Error: Can't write file... skipping!" << endl; cerr << "Operator::DumpMaterial2File: Error: Can't write file... skipping!" << endl;
cout << " done!" << endl; cout << " done!" << endl;

8
openEMS.pro
View File

@ -124,8 +124,8 @@ SOURCES += Common/operator_base.cpp \
tools/array_ops.cpp \ tools/array_ops.cpp \
tools/ErrorMsg.cpp \ tools/ErrorMsg.cpp \
tools/AdrOp.cpp \ tools/AdrOp.cpp \
tools/vtk_file_io.cpp \ tools/vtk_file_writer.cpp \
tools/base_file_io.cpp tools/hdf5_file_writer.cpp
#### HEADERS ################################################################ #### HEADERS ################################################################
HEADERS += openems.h HEADERS += openems.h
@ -188,8 +188,8 @@ HEADERS += tools/ErrorMsg.h \
tools/global.h \ tools/global.h \
tools/useful.h \ tools/useful.h \
tools/aligned_allocator.h \ tools/aligned_allocator.h \
tools/vtk_file_io.h \ tools/vtk_file_writer.h \
tools/base_file_io.h tools/hdf5_file_writer.h
QMAKE_CXXFLAGS_RELEASE = -O3 \ QMAKE_CXXFLAGS_RELEASE = -O3 \
-g \ -g \

33
tools/base_file_io.cpp
View File

@ -1,33 +0,0 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "base_file_io.h"
Base_File_IO::Base_File_IO(std::string filename, int meshType)
{
SetFilename(filename);
m_MeshType = meshType;
m_NativeDump = false;
m_Binary = true;
m_Compress = true;
m_AppendMode = false;
m_ActiveTS = false;
}
Base_File_IO::~Base_File_IO()
{
}

427
tools/hdf5_file_reader.cpp Normal file
View File

@ -0,0 +1,427 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "hdf5_file_reader.h"
#include "../tools/array_ops.h"
#include <hdf5.h>
#include <sstream>
#include <iostream>
#include <iomanip>
using namespace std;
HDF5_File_Reader::HDF5_File_Reader(string filename)
{
m_filename = filename;
//suppress hdf5 error output
//H5Eset_auto(NULL, NULL);
}
HDF5_File_Reader::~HDF5_File_Reader()
{
}
bool HDF5_File_Reader::IsValid()
{
htri_t val = H5Fis_hdf5(m_filename.c_str());
if (val<0)
{
cerr << "HDF5_File_Reader::IsValid: the given file """ << m_filename << """ is not accessible..." << endl;
return false;
}
if (val==0)
{
cerr << "HDF5_File_Reader::IsValid: the given file """ << m_filename << """ is invalid..." << endl;
return false;
}
if (val==0)
cerr << "HDF5_File_Reader::IsValid: the given file """ << m_filename << """ is valid..." << endl;
return true;
}
bool HDF5_File_Reader::ReadDataSet(string ds_name, hsize_t &nDim, hsize_t* &dims, float* &data)
{
if (IsValid()==false)
return false;
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
if (hdf5_file==-1)
{
cerr << "HDF5_File_Reader::ReadDataSet: opening the given file """ << m_filename << """ failed" << endl;
return false;
}
hid_t dataset = H5Dopen(hdf5_file, ds_name.c_str() );
if (dataset<0)
{
cerr << "HDF5_File_Reader::ReadDataSet: dataset not found" << endl;
return false;
}
hid_t type = H5Dget_type(dataset);
if (type<0)
{
cerr << "HDF5_File_Reader::ReadDataSet: dataset type error" << endl;
return false;
}
if (H5Tget_class(type)!=H5T_FLOAT)
{
cerr << "HDF5_File_Reader::ReadDataSet: dataset type not a float" << endl;
return false;
}
hid_t space = H5Dget_space(dataset);
nDim = H5Sget_simple_extent_ndims(space);
dims = new hsize_t[nDim];
H5Sget_simple_extent_dims(space, dims, NULL );
hsize_t data_size = 1;
for (unsigned int d=0;d<nDim;++d)
data_size*=dims[d];
data = new float[data_size];
if (H5Dread(dataset,type,H5S_ALL,H5S_ALL,H5P_DEFAULT,data)<0)
{
cerr << "HDF5_File_Reader::ReadDataSet: error reading data" << endl;
return false;
}
return true;
}
bool HDF5_File_Reader::ReadMesh(float** lines, unsigned int* numLines, int &meshType)
{
if (IsValid()==false)
return false;
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
if (hdf5_file==-1)
{
cerr << "HDF5_File_Reader::ReadMesh: opening the given file """ << m_filename << """ failed" << endl;
return false;
}
vector<string> Mesh_Names;
if (H5Lexists(hdf5_file, "/Mesh/x", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/y", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/z", H5P_DEFAULT))
{
meshType = 0;
Mesh_Names.push_back("/Mesh/x");
Mesh_Names.push_back("/Mesh/y");
Mesh_Names.push_back("/Mesh/z");
}
else if (H5Lexists(hdf5_file, "/Mesh/rho", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/alpha", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/z", H5P_DEFAULT))
{
meshType = 1;
Mesh_Names.push_back("/Mesh/rho");
Mesh_Names.push_back("/Mesh/alpha");
Mesh_Names.push_back("/Mesh/z");
}
else
{
cerr << "HDF5_File_Reader::ReadMesh: no falid mesh information found" << endl;
return false;
}
for (int n=0;n<3;++n)
{
hsize_t nDim;
hsize_t* dims=NULL;
float* data=NULL;
ReadDataSet(Mesh_Names.at(n), nDim, dims, data);
if (nDim!=1)
{
cerr << "HDF5_File_Reader::ReadMesh: mesh dimension error" << endl;
delete[] dims;
return false;
}
numLines[n]=dims[0];
delete[] dims;
lines[n]=data;
}
return true;
}
unsigned int HDF5_File_Reader::GetNumTimeSteps()
{
if (IsValid()==false)
return false;
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
if (hdf5_file==-1)
{
cerr << "HDF5_File_Reader::GetNumTimeSteps: opening the given file """ << m_filename << """ failed" << endl;
return false;
}
if (H5Lexists(hdf5_file, "/FieldData/TD", H5P_DEFAULT)<0)
return 0;
hid_t TD_grp = H5Gopen(hdf5_file, "/FieldData/TD" );
if (TD_grp<0)
{
cerr << "HDF5_File_Reader::GetNumTimeSteps: can't open group ""/FieldData/TD""" << endl;
return 0;
}
hsize_t numObj;
if (H5Gget_num_objs(TD_grp,&numObj)<0)
{
cerr << "HDF5_File_Reader::GetNumTimeSteps: can't read number of datasets" << endl;
return 0;
}
return numObj;
}
bool HDF5_File_Reader::ReadTimeSteps(vector<unsigned int> &timestep, vector<string> &names)
{
if (IsValid()==false)
return false;
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
if (hdf5_file==-1)
{
cerr << "HDF5_File_Reader::ReadTimeSteps: opening the given file """ << m_filename << """ failed" << endl;
return false;
}
if (H5Lexists(hdf5_file, "/FieldData/TD", H5P_DEFAULT)<0)
{
cerr << "HDF5_File_Reader::ReadTimeSteps: can't open ""/FieldData/TD""" << endl;
return false;
}
hid_t TD_grp = H5Gopen(hdf5_file, "/FieldData/TD" );
if (TD_grp<0)
{
cerr << "HDF5_File_Reader::ReadTimeSteps: can't open ""/FieldData/TD""" << endl;
return false;
}
hsize_t numObj;
if (H5Gget_num_objs(TD_grp,&numObj)<0)
{
cerr << "HDF5_File_Reader::ReadTimeSteps: can't read number of datasets" << endl;
return false;
}
char name[100];
timestep.clear();
timestep.resize(numObj,0);
names.clear();
names.resize(numObj);
for (hsize_t n=0;n<numObj;++n)
{
if (H5Gget_objtype_by_idx(TD_grp, n) != H5G_DATASET)
{
cerr << "HDF5_File_Reader::ReadTimeSteps: invalid timestep found!" << endl;
return false;
}
ssize_t name_size = H5Gget_objname_by_idx(TD_grp, n, name, 100 );
istringstream is(name);
unsigned int num;
if (is >> num)
{
timestep.at(n)=num;
names.at(n)=name;
}
else
{
cerr << "HDF5_File_Reader::ReadTimeSteps: invalid timestep format found!" << endl;
return false;
}
}
return true;
}
float**** HDF5_File_Reader::GetTDVectorData(size_t idx, float &time, unsigned int data_size[])
{
if (IsValid()==false)
return false;
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
if (hdf5_file==-1)
{
cerr << "HDF5_File_Reader::GetTDVectorData: opening the given file """ << m_filename << """ failed" << endl;
return false;
}
time = 0;
if (H5Lexists(hdf5_file, "/FieldData/TD", H5P_DEFAULT)<0)
{
cerr << "HDF5_File_Reader::GetTDVectorData: can't open ""/FieldData/TD""" << endl;
return NULL;
}
hid_t TD_grp = H5Gopen(hdf5_file, "/FieldData/TD" );
if (TD_grp<0)
{
cerr << "HDF5_File_Reader::GetTDVectorData: can't open ""/FieldData/TD""" << endl;
return NULL;
}
hsize_t numObj;
if (H5Gget_num_objs(TD_grp,&numObj)<0)
{
cerr << "HDF5_File_Reader::GetTDVectorData: can't read number of datasets" << endl;
return NULL;
}
if (idx>=numObj)
return NULL;
if (H5Gget_objtype_by_idx(TD_grp, idx) != H5G_DATASET)
{
cerr << "HDF5_File_Reader::GetTDVectorData: invalid timestep found!" << endl;
return NULL;
}
char name[100];
ssize_t name_size = H5Gget_objname_by_idx(TD_grp, idx, name, 100 );
string ds_name = "/FieldData/TD/" + string(name);
hid_t attr = H5Aopen_by_name(hdf5_file, ds_name.c_str(), "time", H5P_DEFAULT, H5P_DEFAULT);
if (attr<0)
{
cerr << "HDF5_File_Reader::GetTDVectorData: time attribute not found!" << endl;
return NULL;
}
if (H5Aread(attr, H5T_NATIVE_FLOAT, &time)<0)
{
cerr << "HDF5_File_Reader::GetTDVectorData: can't read time attribute!" << endl;
return NULL;
}
hsize_t nDim;
hsize_t* dims=NULL;
float* data=NULL;
ReadDataSet(ds_name, nDim, dims, data);
if (nDim!=4)
{
cerr << "HDF5_File_Reader::GetTDVectorData: data dimension invalid" << endl;
delete[] dims;
return NULL;
}
if (dims[0]!=3)
{
cerr << "HDF5_File_Reader::GetTDVectorData: vector data dimension invalid" << endl;
delete[] dims;
return NULL;
}
data_size[0]=dims[3];
data_size[1]=dims[2];
data_size[2]=dims[1];
delete[] dims;
data_size[3]=3;
size_t pos = 0;
float**** field = Create_N_3DArray<float>(data_size);
for (unsigned int d=0;d<3;++d)
for (unsigned int k=0;k<data_size[2];++k)
for (unsigned int j=0;j<data_size[1];++j)
for (unsigned int i=0;i<data_size[0];++i)
{
field[d][i][j][k]=data[pos++];
}
delete[] data;
return field;
}
unsigned int HDF5_File_Reader::GetNumFrequencies()
{
cerr << "HDF5_File_Reader::GetNumFrequencies(): not implemented yet!" << endl;
return 0;
}
bool HDF5_File_Reader::ReadFrequencies(vector<float> &frequencies)
{
cerr << "HDF5_File_Reader::ReadFrequencies(): not implemented yet!" << endl;
return false;
}
complex<float>**** HDF5_File_Reader::GetFDVectorData(size_t idx, float &frequency, unsigned int data_size[])
{
cerr << "HDF5_File_Reader::GetFDVectorData(): not implemented yet!" << endl;
return NULL;
}
bool HDF5_File_Reader::CalcFDVectorData(vector<float> &frequencies, vector<complex<float>****> &FD_data, unsigned int data_size[4])
{
FD_data.clear();
float time;
//read first TD data
float**** field = this->GetTDVectorData(0,time,data_size);
if (field==NULL)
{
cerr << "HDF5_File_Reader::CalcFDVectorData: error, no TD data found..." << endl;
return false;
}
//init
FD_data.resize(frequencies.size(), NULL);
for (size_t fn=0;fn<frequencies.size();++fn)
FD_data.at(fn) = Create_N_3DArray<complex<float> >(data_size);
size_t ts=0;
unsigned int pos[3];
complex<float> PI_2_I(0.0,-2.0*M_PI);
complex<float> exp_jwt_2_dt;
float time_diff=0;
float time_old =0;
complex<float>**** field_fd = NULL;
while (field)
{
if ((ts>1) && abs(time_diff - (time - time_old))>1e15)
{
cerr << "HDF5_File_Reader::CalcFDVectorData: time interval error..." << endl;
for (size_t fn=0;fn<frequencies.size();++fn)
Delete_N_3DArray(FD_data.at(fn),data_size);
FD_data.clear();
return false;
}
time_diff = time - time_old;
for (size_t fn=0;fn<frequencies.size();++fn)
{
exp_jwt_2_dt = exp( (complex<float>)(PI_2_I * frequencies.at(fn) * time) );
field_fd = FD_data.at(fn);
for (pos[0]=0; pos[0]<data_size[0]; ++pos[0])
{
for (pos[1]=0; pos[1]<data_size[1]; ++pos[1])
{
for (pos[2]=0; pos[2]<data_size[2]; ++pos[2])
{
field_fd[0][pos[0]][pos[1]][pos[2]] += field[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
field_fd[1][pos[0]][pos[1]][pos[2]] += field[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
field_fd[2][pos[0]][pos[1]][pos[2]] += field[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
}
}
}
}
++ts;
Delete_N_3DArray(field,data_size);
time_old = time;
field = this->GetTDVectorData(ts,time,data_size);
}
// finalize data
time_diff*=2;
for (size_t fn=0;fn<frequencies.size();++fn)
{
field_fd = FD_data.at(fn);
for (pos[0]=0; pos[0]<data_size[0]; ++pos[0])
{
for (pos[1]=0; pos[1]<data_size[1]; ++pos[1])
{
for (pos[2]=0; pos[2]<data_size[2]; ++pos[2])
{
field_fd[0][pos[0]][pos[1]][pos[2]] *= time_diff;
field_fd[1][pos[0]][pos[1]][pos[2]] *= time_diff;
field_fd[2][pos[0]][pos[1]][pos[2]] *= time_diff;
}
}
}
}
return true;
}

64
tools/hdf5_file_reader.h Normal file
View File

@ -0,0 +1,64 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef HDF5_FILE_READER_H
#define HDF5_FILE_READER_H
#include <string>
#include <vector>
#include <complex>
#include <hdf5.h>
class HDF5_File_Reader
{
public:
HDF5_File_Reader(std::string filename);
virtual ~HDF5_File_Reader();
bool ReadMesh(float** lines, unsigned int* numLines, int &meshType);
//! Get the number of timesteps stored at /FieldData/TD/<NUMBER_OF_TS>
unsigned int GetNumTimeSteps();
bool ReadTimeSteps(std::vector<unsigned int> &timestep, std::vector<std::string> &names);
/*!
Get time-domain data stored at /FieldData/TD/<NUMBER_OF_TS>
\param[in] ids time step index to extract
\param[out] time time attribute for the given timestep
\param[out] data_size data size found
\return field data found in given timestep, caller must delete array, returns NULL if timestep was not found
*/
float**** GetTDVectorData(size_t idx, float &time, unsigned int data_size[4]);
unsigned int GetNumFrequencies();
bool ReadFrequencies(std::vector<float> &frequencies);
std::complex<float>**** GetFDVectorData(size_t idx, float &frequency, unsigned int data_size[4]);
/*!
Calculate
*/
bool CalcFDVectorData(std::vector<float> &frequencies, std::vector<std::complex<float>****> &FD_data, unsigned int data_size[4]);
bool IsValid();
protected:
std::string m_filename;
bool ReadDataSet(std::string ds_name, hsize_t &nDim, hsize_t* &dims, float* &data);
};
#endif // HDF5_FILE_READER_H

335
tools/hdf5_file_writer.cpp Normal file
View File

@ -0,0 +1,335 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
using namespace std;
#include "hdf5_file_writer.h"
#include <boost/algorithm/string.hpp>
#include <hdf5.h>
#include <sstream>
#include <iostream>
#include <iomanip>
HDF5_File_Writer::HDF5_File_Writer(string filename)
{
m_filename = filename;
m_Group = "/";
hid_t hdf5_file = H5Fcreate(m_filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
if (hdf5_file<0)
{
cerr << "HDF5_File_Writer::HDF5_File_Writer: Error, creating the given file """ << m_filename << """ failed" << endl;
}
H5Fclose(hdf5_file);
}
HDF5_File_Writer::~HDF5_File_Writer()
{
}
hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group)
{
if (hdf5_file<0)
{
cerr << "HDF5_File_Writer::CreateGroup: Error, invalid file id" << endl;
return -1;
}
vector<string> results;
boost::split(results, group, boost::is_any_of("/"));
hid_t grp=H5Gopen(hdf5_file,"/");
if (grp<0)
{
cerr << "HDF5_File_Writer::OpenGroup: Error, opening root group " << endl;
return -1;
}
for (size_t n=0;n<results.size();++n)
{
if (!results.at(n).empty())
{
if (H5Lexists(grp, results.at(n).c_str(), H5P_DEFAULT))
{
grp = H5Gopen(grp, results.at(n).c_str());
if (grp<0)
{
cerr << "HDF5_File_Writer::OpenGroup: Error, failed to open existing group" << endl;
return -1;
}
}
else
{
grp = H5Gcreate(grp,results.at(n).c_str(),0);
if (grp<0)
{
cerr << "HDF5_File_Writer::OpenGroup: Error, creating group """ << group << """ failed" << endl;
return -1;
}
}
}
}
return grp;
}
bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, double const* const* discLines, int MeshType, double scaling)
{
float* array[3];
for (int n=0; n<3; ++n)
{
array[n] = new float[numLines[n]];
for (unsigned int i=0; i<numLines[n]; ++i)
array[n][i]=discLines[n][i];
}
bool success = WriteRectMesh(numLines,array,MeshType,scaling);
for (int n=0; n<3; ++n)
delete[] array[n];
return success;
}
bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, float const* const* discLines, int MeshType, float scaling)
{
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
if (hdf5_file<0)
{
cerr << "HDF5_File_Writer::WriteRectMesh: Error, opening the given file """ << m_filename << """ failed" << endl;
return false;
}
if (H5Lexists(hdf5_file, "/Mesh", H5P_DEFAULT))
{
cerr << "HDF5_File_Writer::WriteRectMesh: Error, group ""/Mesh"" already exists" << endl;
return false;
}
hid_t mesh_grp = H5Gcreate(hdf5_file,"/Mesh",0);
if (mesh_grp<0)
{
cerr << "HDF5_File_Writer::WriteRectMesh: Error, creating group ""/Mesh"" failed" << endl;
return false;
}
string names[] = {"x","y","z"};
if (MeshType==1)
{
names[0]="rho";
names[1]="alpha";
}
if (MeshType==2)
{
names[0]="r";
names[1]="theta";
names[2]="phi";
}
for (int n=0; n<3; ++n)
{
hsize_t dims[1]={numLines[n]};
hid_t space = H5Screate_simple(1, dims, NULL);
hid_t dataset = H5Dcreate(mesh_grp, names[n].c_str(), H5T_NATIVE_FLOAT, space, H5P_DEFAULT);
float* array = new float[numLines[n]];
for (unsigned int i=0; i<numLines[n]; ++i)
{
if ((MeshType==1) && (n==1)) //check for alpha-direction
array[i] = discLines[n][i];
else if ((MeshType==2) && (n>0)) //check for theta/phi-direction
array[i] = discLines[n][i];
else
array[i] = discLines[n][i] * scaling;
}
if (H5Dwrite(dataset, H5T_NATIVE_FLOAT, space, H5P_DEFAULT, H5P_DEFAULT, array))
{
cerr << "HDF5_File_Writer::WriteRectMesh: Error, writing to dataset failed" << endl;
delete[] array;
return false;
}
delete[] array;
}
H5Fclose(hdf5_file);
return true;
}
bool HDF5_File_Writer::WriteScalarField(std::string dataSetName, float const* const* const* field, size_t datasize[3])
{
size_t pos = 0;
size_t size = datasize[0]*datasize[1]*datasize[2];
size_t n_size[3]={datasize[2],datasize[1],datasize[0]};
float* buffer = new float[size];
for (size_t k=0;k<datasize[2];++k)
for (size_t j=0;j<datasize[1];++j)
for (size_t i=0;i<datasize[0];++i)
{
buffer[pos++]=field[i][j][k];
}
bool success = WriteData(dataSetName,buffer,3,n_size);
delete[] buffer;
return success;
}
bool HDF5_File_Writer::WriteScalarField(std::string dataSetName, std::complex<float> const* const* const* field, size_t datasize[3])
{
size_t pos = 0;
size_t size = datasize[0]*datasize[1]*datasize[2];
size_t n_size[3]={datasize[2],datasize[1],datasize[0]};
float* buffer = new float[size];
for (size_t k=0;k<datasize[2];++k)
for (size_t j=0;j<datasize[1];++j)
for (size_t i=0;i<datasize[0];++i)
{
buffer[pos++]=real(field[i][j][k]);
}
bool success = WriteData(dataSetName + "_real",buffer,3,n_size);
pos = 0;
for (size_t k=0;k<datasize[2];++k)
for (size_t j=0;j<datasize[1];++j)
for (size_t i=0;i<datasize[0];++i)
{
buffer[pos++]=imag(field[i][j][k]);
}
success &= WriteData(dataSetName + "_imag",buffer,3,n_size);
delete[] buffer;
return success;
}
bool HDF5_File_Writer::WriteVectorField(std::string dataSetName, float const* const* const* const* field, size_t datasize[3])
{
size_t pos = 0;
size_t size = datasize[0]*datasize[1]*datasize[2]*3;
float* buffer = new float[size];
for (int n=0;n<3;++n)
for (size_t k=0;k<datasize[2];++k)
for (size_t j=0;j<datasize[1];++j)
for (size_t i=0;i<datasize[0];++i)
{
buffer[pos++]=field[n][i][j][k];
}
size_t n_size[4]={3,datasize[2],datasize[1],datasize[0]};
bool success = WriteData(dataSetName,buffer,4,n_size);
delete[] buffer;
return success;
}
bool HDF5_File_Writer::WriteVectorField(std::string dataSetName, std::complex<float> const* const* const* const* field, size_t datasize[3])
{
size_t pos = 0;
size_t size = datasize[0]*datasize[1]*datasize[2]*3;
size_t n_size[4]={3,datasize[2],datasize[1],datasize[0]};
float* buffer = new float[size];
for (int n=0;n<3;++n)
for (size_t k=0;k<datasize[2];++k)
for (size_t j=0;j<datasize[1];++j)
for (size_t i=0;i<datasize[0];++i)
{
buffer[pos++]=real(field[n][i][j][k]);
}
bool success = WriteData(dataSetName + "_real",buffer,4,n_size);
pos = 0;
for (int n=0;n<3;++n)
for (size_t k=0;k<datasize[2];++k)
for (size_t j=0;j<datasize[1];++j)
for (size_t i=0;i<datasize[0];++i)
{
buffer[pos++]=imag(field[n][i][j][k]);
}
success &= WriteData(dataSetName + "_imag",buffer,4,n_size);
delete[] buffer;
return success;
}
bool HDF5_File_Writer::WriteData(std::string dataSetName, float const* field_buf, size_t dim, size_t* datasize)
{
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
if (hdf5_file<0)
{
cerr << "HDF5_File_Writer::WriteData: Error, opening the given file """ << m_filename << """ failed" << endl;
return false;
}
hid_t group = OpenGroup(hdf5_file,m_Group);
hsize_t dims[dim];
for (size_t n=0;n<dim;++n)
dims[n]=datasize[n];
hid_t space = H5Screate_simple(dim, dims, NULL);
hid_t dataset = H5Dcreate(group, dataSetName.c_str(), H5T_NATIVE_FLOAT, space, H5P_DEFAULT);
if (H5Dwrite(dataset, H5T_NATIVE_FLOAT, space, H5P_DEFAULT, H5P_DEFAULT, field_buf))
{
cerr << "HDF5_File_Writer::WriteData: Error, writing to dataset failed" << endl;
return false;
}
H5Fclose(hdf5_file);
return true;
}
bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, float const* value, hsize_t size)
{
return WriteAtrribute(locName,attr_name, value, size, H5T_NATIVE_FLOAT);
}
bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, void const* value, hsize_t size, hid_t mem_type)
{
hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
if (hdf5_file<0)
{
cerr << "HDF5_File_Writer::WriteAtrribute: Error, opening the given file """ << m_filename << """ failed" << endl;
return false;
}
if (H5Lexists(hdf5_file, locName.c_str(), H5P_DEFAULT)<0)
{
cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to find location: """ << locName << """" << endl;
return false;
}
hid_t loc = H5Oopen(hdf5_file, locName.c_str(), H5P_DEFAULT);
if (loc<0)
{
cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to open location: """ << locName << """" << endl;
return false;
}
hid_t dataspace_id = H5Screate_simple(1, &size, NULL);
/* Create a dataset attribute. */
hid_t attribute_id = H5Acreate(loc, attr_name.c_str(), mem_type, dataspace_id,H5P_DEFAULT);
if (attribute_id<0)
{
cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to create the attrbute" << endl;
return false;
}
/* Write the attribute data. */
if (H5Awrite(attribute_id, mem_type, value)<0)
{
cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to write the attrbute" << endl;
return false;
}
return true;
}
bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, vector<float> values)
{
float val[values.size()];
for (size_t n=0;n<values.size();++n)
val[n]=values.at(n);
return HDF5_File_Writer::WriteAtrribute(locName, attr_name,val,values.size());
}

55
tools/hdf5_file_writer.h Normal file
View File

@ -0,0 +1,55 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef HDF5_FILE_WRITER_H
#define HDF5_FILE_WRITER_H
#include <string>
#include <vector>
#include <complex>
#include <hdf5.h>
class HDF5_File_Writer
{
public:
HDF5_File_Writer(std::string filename);
~HDF5_File_Writer();
bool WriteRectMesh(unsigned int const* numLines, double const* const* discLines, int MeshType=0, double scaling=1);
bool WriteRectMesh(unsigned int const* numLines, float const* const* discLines, int MeshType=0, float scaling=1);
bool WriteScalarField(std::string dataSetName, float const* const* const* field, size_t datasize[3]);
bool WriteScalarField(std::string dataSetName, std::complex<float> const* const* const* field, size_t datasize[3]);
bool WriteVectorField(std::string dataSetName, float const* const* const* const* field, size_t datasize[3]);
bool WriteVectorField(std::string dataSetName, std::complex<float> const* const* const* const* field, size_t datasize[3]);
bool WriteData(std::string dataSetName, float const* field_buf, size_t dim, size_t* datasize);
bool WriteAtrribute(std::string locName, std::string attr_name, void const* value, hsize_t size, hid_t mem_type);
bool WriteAtrribute(std::string locName, std::string attr_name, float const* value, hsize_t size);
bool WriteAtrribute(std::string locName, std::string attr_name, std::vector<float> values);
void SetCurrentGroup(std::string group) {m_Group = group;}
protected:
std::string m_filename;
std::string m_Group;
hid_t OpenGroup(hid_t hdf5_file, std::string group);
};
#endif // HDF5_FILE_WRITER_H

54
tools/vtk_file_io.h
View File

@ -1,54 +0,0 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef VTK_FILE_IO_H
#define VTK_FILE_IO_H
#include "base_file_io.h"
class vtkDataSet;
class VTK_File_IO : public Base_File_IO
{
public:
VTK_File_IO(std::string filename, int meshType=0);
virtual ~VTK_File_IO();
virtual void SetMeshLines(double const* const* lines, unsigned int const* count, double scaling=1);
virtual void AddScalarField(std::string fieldname, double const* const* const* field, unsigned int const* size);
virtual void AddScalarField(std::string fieldname, float const* const* const* field, unsigned int const* size);
virtual void AddVectorField(std::string fieldname, double const* const* const* const* field, unsigned int const* size);
virtual void AddVectorField(std::string fieldname, float const* const* const* const* field, unsigned int const* size);
virtual int GetNumberOfFields() const;
virtual void ClearAllFields();
virtual bool Write();
virtual bool WriteASCII();
virtual bool WriteXML();
protected:
vtkDataSet* m_GridData;
std::vector<double> m_MeshLines[3];
virtual std::string GetTimestepFilename(int pad_length=10) const;
};
#endif // VTK_FILE_IO_H

82
tools/vtk_file_io.cpp → tools/vtk_file_writer.cpp
View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de) * Copyright (C) 2011,2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
* *
* This program is free software: you can redistribute it and/or modify * This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by * it under the terms of the GNU General Public License as published by
@ -17,7 +17,7 @@
using namespace std; using namespace std;
#include "vtk_file_io.h" #include "vtk_file_writer.h"
#include <vtkRectilinearGrid.h> #include <vtkRectilinearGrid.h>
#include <vtkRectilinearGridWriter.h> #include <vtkRectilinearGridWriter.h>
@ -35,34 +35,42 @@ using namespace std;
#include <iomanip> #include <iomanip>
VTK_File_IO::VTK_File_IO(string filename, int meshType) : Base_File_IO(filename, meshType) VTK_File_Writer::VTK_File_Writer(string filename, int meshType)
{ {
SetFilename(filename);
m_MeshType = meshType;
m_NativeDump = false;
m_Binary = true;
m_Compress = true;
m_AppendMode = false;
m_ActiveTS = false;
if (m_MeshType==0) //cartesian mesh if (m_MeshType==0) //cartesian mesh
m_GridData = vtkRectilinearGrid::New(); m_GridData = vtkRectilinearGrid::New();
else if (m_MeshType==1) //cylindrical mesh else if (m_MeshType==1) //cylindrical mesh
m_GridData = vtkStructuredGrid::New(); m_GridData = vtkStructuredGrid::New();
else else
{ {
cerr << "VTK_File_IO::VTK_File_IO: Error, unknown mesh type: " << m_MeshType << endl; cerr << "VTK_File_Writer::VTK_File_Writer: Error, unknown mesh type: " << m_MeshType << endl;
m_GridData=NULL; m_GridData=NULL;
} }
} }
VTK_File_IO::~VTK_File_IO() VTK_File_Writer::~VTK_File_Writer()
{ {
if (m_GridData) if (m_GridData)
m_GridData->Delete(); m_GridData->Delete();
m_GridData = NULL; m_GridData = NULL;
} }
void VTK_File_IO::SetMeshLines(double const* const* lines, unsigned int const* count, double scaling) void VTK_File_Writer::SetMeshLines(double const* const* lines, unsigned int const* count, double scaling)
{ {
if (m_MeshType==0) //cartesian mesh if (m_MeshType==0) //cartesian mesh
{ {
vtkRectilinearGrid* RectGrid = dynamic_cast<vtkRectilinearGrid*>(m_GridData); vtkRectilinearGrid* RectGrid = dynamic_cast<vtkRectilinearGrid*>(m_GridData);
if (RectGrid==NULL) if (RectGrid==NULL)
{ {
cerr << "VTK_File_IO::SetMeshLines: Error, grid invalid, this should not have happend! " << endl; cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
exit(1); exit(1);
} }
RectGrid->SetDimensions(count[0],count[1],count[2]); RectGrid->SetDimensions(count[0],count[1],count[2]);
@ -89,7 +97,7 @@ void VTK_File_IO::SetMeshLines(double const* const* lines, unsigned int const* c
vtkStructuredGrid* StructGrid = dynamic_cast<vtkStructuredGrid*>(m_GridData); vtkStructuredGrid* StructGrid = dynamic_cast<vtkStructuredGrid*>(m_GridData);
if (StructGrid==NULL) if (StructGrid==NULL)
{ {
cerr << "VTK_File_IO::SetMeshLines: Error, grid invalid, this should not have happend! " << endl; cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
exit(1); exit(1);
} }
@ -123,21 +131,21 @@ void VTK_File_IO::SetMeshLines(double const* const* lines, unsigned int const* c
} }
else else
{ {
cerr << "VTK_File_IO::SetMeshLines: Error, unknown mesh type: " << m_MeshType << endl; cerr << "VTK_File_Writer::SetMeshLines: Error, unknown mesh type: " << m_MeshType << endl;
} }
} }
void VTK_File_IO::AddScalarField(string fieldname, double const* const* const* field, unsigned int const* size) void VTK_File_Writer::AddScalarField(string fieldname, double const* const* const* field)
{ {
vtkDoubleArray* array = vtkDoubleArray::New(); vtkDoubleArray* array = vtkDoubleArray::New();
array->SetNumberOfTuples(size[0]*size[1]*size[2]); array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
array->SetName(fieldname.c_str()); array->SetName(fieldname.c_str());
int id=0; int id=0;
for (unsigned int k=0;k<size[2];++k) for (unsigned int k=0;k<m_MeshLines[2].size();++k)
{ {
for (unsigned int j=0;j<size[1];++j) for (unsigned int j=0;j<m_MeshLines[1].size();++j)
{ {
for (unsigned int i=0;i<size[0];++i) for (unsigned int i=0;i<m_MeshLines[0].size();++i)
{ {
array->SetTuple1(id++,field[i][j][k]); array->SetTuple1(id++,field[i][j][k]);
} }
@ -147,17 +155,17 @@ void VTK_File_IO::AddScalarField(string fieldname, double const* const* const* f
array->Delete(); array->Delete();
} }
void VTK_File_IO::AddScalarField(string fieldname, float const* const* const* field, unsigned int const* size) void VTK_File_Writer::AddScalarField(string fieldname, float const* const* const* field)
{ {
vtkFloatArray* array = vtkFloatArray::New(); vtkFloatArray* array = vtkFloatArray::New();
array->SetNumberOfTuples(size[0]*size[1]*size[2]); array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
array->SetName(fieldname.c_str()); array->SetName(fieldname.c_str());
int id=0; int id=0;
for (unsigned int k=0;k<size[2];++k) for (unsigned int k=0;k<m_MeshLines[2].size();++k)
{ {
for (unsigned int j=0;j<size[1];++j) for (unsigned int j=0;j<m_MeshLines[1].size();++j)
{ {
for (unsigned int i=0;i<size[0];++i) for (unsigned int i=0;i<m_MeshLines[0].size();++i)
{ {
array->SetTuple1(id++,field[i][j][k]); array->SetTuple1(id++,field[i][j][k]);
} }
@ -167,21 +175,21 @@ void VTK_File_IO::AddScalarField(string fieldname, float const* const* const* fi
array->Delete(); array->Delete();
} }
void VTK_File_IO::AddVectorField(string fieldname, double const* const* const* const* field, unsigned int const* size) void VTK_File_Writer::AddVectorField(string fieldname, double const* const* const* const* field)
{ {
vtkDoubleArray* array = vtkDoubleArray::New(); vtkDoubleArray* array = vtkDoubleArray::New();
array->SetNumberOfComponents(3); array->SetNumberOfComponents(3);
array->SetNumberOfTuples(size[0]*size[1]*size[2]); array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
array->SetName(fieldname.c_str()); array->SetName(fieldname.c_str());
int id=0; int id=0;
double out[3]; double out[3];
for (unsigned int k=0;k<size[2];++k) for (unsigned int k=0;k<m_MeshLines[2].size();++k)
{ {
for (unsigned int j=0;j<size[1];++j) for (unsigned int j=0;j<m_MeshLines[1].size();++j)
{ {
double cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh) double cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
double sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh) double sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
for (unsigned int i=0;i<size[0];++i) for (unsigned int i=0;i<m_MeshLines[0].size();++i)
{ {
if ((m_MeshType==0) || (m_NativeDump)) if ((m_MeshType==0) || (m_NativeDump))
array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]); array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
@ -199,21 +207,21 @@ void VTK_File_IO::AddVectorField(string fieldname, double const* const* const* c
array->Delete(); array->Delete();
} }
void VTK_File_IO::AddVectorField(string fieldname, float const* const* const* const* field, unsigned int const* size) void VTK_File_Writer::AddVectorField(string fieldname, float const* const* const* const* field)
{ {
vtkFloatArray* array = vtkFloatArray::New(); vtkFloatArray* array = vtkFloatArray::New();
array->SetNumberOfComponents(3); array->SetNumberOfComponents(3);
array->SetNumberOfTuples(size[0]*size[1]*size[2]); array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
array->SetName(fieldname.c_str()); array->SetName(fieldname.c_str());
int id=0; int id=0;
float out[3]; float out[3];
for (unsigned int k=0;k<size[2];++k) for (unsigned int k=0;k<m_MeshLines[2].size();++k)
{ {
for (unsigned int j=0;j<size[1];++j) for (unsigned int j=0;j<m_MeshLines[1].size();++j)
{ {
float cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh) float cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
float sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh) float sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
for (unsigned int i=0;i<size[0];++i) for (unsigned int i=0;i<m_MeshLines[0].size();++i)
{ {
if ((m_MeshType==0) || (m_NativeDump)) if ((m_MeshType==0) || (m_NativeDump))
array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]); array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
@ -232,12 +240,12 @@ void VTK_File_IO::AddVectorField(string fieldname, float const* const* const* co
} }
int VTK_File_IO::GetNumberOfFields() const int VTK_File_Writer::GetNumberOfFields() const
{ {
return m_GridData->GetPointData()->GetNumberOfArrays(); return m_GridData->GetPointData()->GetNumberOfArrays();
} }
void VTK_File_IO::ClearAllFields() void VTK_File_Writer::ClearAllFields()
{ {
while (m_GridData->GetPointData()->GetNumberOfArrays()>0) while (m_GridData->GetPointData()->GetNumberOfArrays()>0)
{ {
@ -246,12 +254,12 @@ void VTK_File_IO::ClearAllFields()
} }
} }
bool VTK_File_IO::Write() bool VTK_File_Writer::Write()
{ {
return WriteXML(); return WriteXML();
} }
string VTK_File_IO::GetTimestepFilename(int pad_length) const string VTK_File_Writer::GetTimestepFilename(int pad_length) const
{ {
if (m_ActiveTS==false) if (m_ActiveTS==false)
return m_filename; return m_filename;
@ -263,7 +271,7 @@ string VTK_File_IO::GetTimestepFilename(int pad_length) const
} }
bool VTK_File_IO::WriteASCII() bool VTK_File_Writer::WriteASCII()
{ {
vtkDataWriter* writer = NULL; vtkDataWriter* writer = NULL;
if (m_MeshType==0) //cartesian mesh if (m_MeshType==0) //cartesian mesh
@ -272,7 +280,7 @@ bool VTK_File_IO::WriteASCII()
writer = vtkStructuredGridWriter::New(); writer = vtkStructuredGridWriter::New();
else else
{ {
cerr << "VTK_File_IO::WriteASCII: Error, unknown mesh type: " << m_MeshType << endl; cerr << "VTK_File_Writer::WriteASCII: Error, unknown mesh type: " << m_MeshType << endl;
return false; return false;
} }
@ -291,7 +299,7 @@ bool VTK_File_IO::WriteASCII()
return true; return true;
} }
bool VTK_File_IO::WriteXML() bool VTK_File_Writer::WriteXML()
{ {
vtkXMLStructuredDataWriter* writer = NULL; vtkXMLStructuredDataWriter* writer = NULL;
if (m_MeshType==0) //cartesian mesh if (m_MeshType==0) //cartesian mesh
@ -300,7 +308,7 @@ bool VTK_File_IO::WriteXML()
writer = vtkXMLStructuredGridWriter::New(); writer = vtkXMLStructuredGridWriter::New();
else else
{ {
cerr << "VTK_File_IO::WriteXML: Error, unknown mesh type: " << m_MeshType << endl; cerr << "VTK_File_Writer::WriteXML: Error, unknown mesh type: " << m_MeshType << endl;
return false; return false;
} }

54
tools/base_file_io.h → tools/vtk_file_writer.h
View File

@ -1,5 +1,5 @@
/* /*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de) * Copyright (C) 2011,2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
* *
* This program is free software: you can redistribute it and/or modify * This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by * it under the terms of the GNU General Public License as published by
@ -15,20 +15,23 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * along with this program. If not, see <http://www.gnu.org/licenses/>.
*/ */
#ifndef BASE_FILE_IO_H #ifndef VTK_FILE_WRITER_H
#define BASE_FILE_IO_H #define VTK_FILE_WRITER_H
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <iostream> #include <iostream>
#include <string> #include <string>
#include <vector> #include <vector>
#include <complex>
//! Abstract base class for dumping scalar or vector field data class vtkDataSet;
class Base_File_IO
class VTK_File_Writer
{ {
public: public:
virtual ~Base_File_IO(); VTK_File_Writer(std::string filename, int meshType=0);
virtual ~VTK_File_Writer();
//! Set the filename //! Set the filename
virtual void SetFilename(std::string filename) {m_filename=filename;} virtual void SetFilename(std::string filename) {m_filename=filename;}
@ -42,24 +45,17 @@ public:
//! Set compression flag (if the file type supports it) //! Set compression flag (if the file type supports it)
virtual void SetCompress(bool val) {m_Compress=val;} virtual void SetCompress(bool val) {m_Compress=val;}
//! Set the mesh lines for the given mesh type. void SetNativeDump(bool val) {m_NativeDump=val;}
virtual void SetMeshLines(double const* const* lines, unsigned int const* count, double scaling=1) = 0 ;
void SetNativeDump(bool val) {m_NativeDump=val;}; virtual void SetMeshLines(double const* const* lines, unsigned int const* count, double scaling=1);
//! Add a scalar field. \sa GetNumberOfFields \sa ClearAllFields virtual void AddScalarField(std::string fieldname, double const* const* const* field);
virtual void AddScalarField(std::string fieldname, double const* const* const* field, unsigned int const* size) = 0; virtual void AddScalarField(std::string fieldname, float const* const* const* field);
//! Add a scalar field. \sa GetNumberOfFields \sa ClearAllFields virtual void AddVectorField(std::string fieldname, double const* const* const* const* field);
virtual void AddScalarField(std::string fieldname, float const* const* const* field, unsigned int const* size) = 0; virtual void AddVectorField(std::string fieldname, float const* const* const* const* field);
//! Add a vector field. \sa GetNumberOfFields \sa ClearAllFields
virtual void AddVectorField(std::string fieldname, double const* const* const* const* field, unsigned int const* size) = 0;
//! Add a vector field. \sa GetNumberOfFields \sa ClearAllFields
virtual void AddVectorField(std::string fieldname, float const* const* const* const* field, unsigned int const* size) = 0;
//! Get the number of fields. \sa ClearAllFields virtual int GetNumberOfFields() const;
virtual int GetNumberOfFields() const = 0; virtual void ClearAllFields();
//! Clear all included fields. \sa GetNumberOfFields
virtual void ClearAllFields() = 0;
//! Get if timestep file series is active. \sa SetTimestepActive //! Get if timestep file series is active. \sa SetTimestepActive
virtual bool GetTimestepActive() {return m_ActiveTS;} virtual bool GetTimestepActive() {return m_ActiveTS;}
@ -68,23 +64,31 @@ public:
//! Set the current timestep, this will set the timestep flag to true. \sa SetTimestepActive //! Set the current timestep, this will set the timestep flag to true. \sa SetTimestepActive
virtual void SetTimestep(unsigned int ts) {m_timestep=ts;SetTimestepActive(true);} virtual void SetTimestep(unsigned int ts) {m_timestep=ts;SetTimestepActive(true);}
virtual bool Write() = 0; virtual bool Write();
virtual bool WriteASCII();
virtual bool WriteXML();
protected: protected:
Base_File_IO(std::string filename, int meshType=0);
std::string m_filename; std::string m_filename;
std::string m_header; std::string m_header;
//timestep properties
bool m_ActiveTS; bool m_ActiveTS;
unsigned int m_timestep; unsigned int m_timestep;
vtkDataSet* m_GridData;
//mesh information
int m_MeshType; int m_MeshType;
std::vector<double> m_MeshLines[3];
bool m_NativeDump; bool m_NativeDump;
bool m_AppendMode; bool m_AppendMode;
bool m_Binary; bool m_Binary;
bool m_Compress; bool m_Compress;
virtual std::string GetTimestepFilename(int pad_length=10) const;
}; };
#endif // VTK_FILE_Writer_H
#endif // BASE_FILE_IO_H