new ProcessFieldsFD class + fixes

changes:
- new field processing class producing frequency domain field-dumps
- Processing needs dualMesh and dualTime information
- fix in TD field dumps using always dual-Time

Signed-off-by: Thorsten Liebig <thorsten.liebig@gmx.de>
pull/1/head
Thorsten Liebig 2010-12-19 20:41:08 +01:00
parent 4254337ea0
commit 0973f80680
13 changed files with 384 additions and 39 deletions

View File

@ -22,8 +22,6 @@
ProcessCurrent::ProcessCurrent(Engine_Interface_Base* eng_if) : ProcessIntegral(eng_if)
{
m_TimeShift = Op->GetTimestep()/2.0;
m_dualMesh = true;
}
ProcessCurrent::~ProcessCurrent()

View File

@ -28,6 +28,7 @@ ProcessFields::ProcessFields(Engine_Interface_Base* eng_if) : Processing(eng_if)
m_fileType = VTK_FILETYPE;
SetSubSampling(1);
SetPrecision(6);
m_dualTime = false;
for (int n=0; n<3; ++n)
{
@ -84,6 +85,8 @@ void ProcessFields::InitProcess()
group = new H5::Group( file->createGroup( "/FieldData" ));
delete group;
group = new H5::Group( file->createGroup( "/FieldData/FD" ));
delete group;
delete file;
}
}
@ -461,6 +464,76 @@ bool ProcessFields::DumpVectorArray2HDF5(string filename, string name, FDTD_FLOA
return true;
}
bool ProcessFields::DumpVectorArray2HDF5(string filename, 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( "/FieldData/FD" ));
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???
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].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???
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].imag() * weight;
}
}
}
}
dataset.write( hdf5array, H5::PredType::NATIVE_FLOAT );
return true;
}
FDTD_FLOAT**** ProcessFields::CalcField()
{
unsigned int pos[3];
@ -486,6 +559,7 @@ FDTD_FLOAT**** ProcessFields::CalcField()
}
}
}
return field;
}
if (m_DumpType==H_FIELD_DUMP)
@ -506,7 +580,10 @@ FDTD_FLOAT**** ProcessFields::CalcField()
}
}
}
}
return field;
}
cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl;
return field;
}

View File

@ -71,8 +71,12 @@ public:
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, double discLines_scaling = 1);
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, double discLines_scaling = 1);
//! Dump a time-domain vector dump to an HDF5 file
static bool DumpVectorArray2HDF5(string filename, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, float time=0);
//! Dump a frequency-domain complex-vector dump to an HDF5 file
static bool DumpVectorArray2HDF5(string filename, string name, std::complex<float> const* const* const* const* array, unsigned int const* numLines, float weight, float frequency);
double CalcTotalEnergy() const;
void SetFileType(FileType fileType) {m_fileType=fileType;}

205
Common/processfields_fd.cpp Normal file
View File

@ -0,0 +1,205 @@
/*
* Copyright (C) 2010 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 "processfields_fd.h"
#include "Common/operator_base.h"
#include <iomanip>
#include <sstream>
#include <string>
ProcessFieldsFD::ProcessFieldsFD(Engine_Interface_Base* eng_if) : ProcessFields(eng_if)
{
}
ProcessFieldsFD::~ProcessFieldsFD()
{
for (size_t n = 0; n<m_FD_Fields.size(); ++n)
{
Delete_N_3DArray(m_FD_Fields.at(n),numLines);
}
m_FD_Fields.clear();
}
void ProcessFieldsFD::InitProcess()
{
if (Enabled==false) return;
if (m_FD_Samples.size()==0)
{
cerr << "ProcessFieldsFD::InitProcess: No frequencies found... skipping this dump!" << endl;
Enabled=false;
return;
}
//setup the hdf5 file
ProcessFields::InitProcess();
//create data structures...
for (size_t n = 0; n<m_FD_Samples.size(); ++n)
{
std::complex<float>**** field_fd = Create_N_3DArray<std::complex<float> >(numLines);
m_FD_Fields.push_back(field_fd);
}
}
int ProcessFieldsFD::Process()
{
if (Enabled==false) return -1;
if (CheckTimestep()==false) return GetNextInterval();
if ((m_FD_Interval==0) || (m_Eng_Interface->GetNumberOfTimesteps()%m_FD_Interval!=0))
return GetNextInterval();
FDTD_FLOAT**** field_td = CalcField();
std::complex<float>**** field_fd = NULL;
double T = m_Eng_Interface->GetTime(m_dualTime);
unsigned int pos[3];
for (size_t n = 0; n<m_FD_Samples.size(); ++n)
{
std::complex<float> exp_jwt = std::exp( (std::complex<float>)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) );
field_fd = m_FD_Fields.at(n);
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
field_fd[0][pos[0]][pos[1]][pos[2]] += field_td[0][pos[0]][pos[1]][pos[2]] * exp_jwt;
field_fd[1][pos[0]][pos[1]][pos[2]] += field_td[1][pos[0]][pos[1]][pos[2]] * exp_jwt;
field_fd[2][pos[0]][pos[1]][pos[2]] += field_td[2][pos[0]][pos[1]][pos[2]] * exp_jwt;
}
}
}
}
++m_FD_SampleCount;
return GetNextInterval();
}
void ProcessFieldsFD::PostProcess()
{
DumpFDData();
}
void ProcessFieldsFD::DumpFDData()
{
#ifdef OUTPUT_IN_DRAWINGUNITS
double discLines_scaling = 1;
#else
double discLines_scaling = Op->GetGridDelta();
#endif
if (m_fileType==VTK_FILETYPE)
{
unsigned int pos[3];
FDTD_FLOAT**** field = Create_N_3DArray<float>(numLines);
std::complex<float>**** field_fd = NULL;
double angle=0;
int Nr_Ph = 21;
for (size_t n = 0; n<m_FD_Samples.size(); ++n)
{
//dump multiple phase to vtk-files
for (int p=0; p<Nr_Ph; ++p)
{
angle = 2.0 * M_PI * p / Nr_Ph;
std::complex<float> exp_jwt = std::exp( (std::complex<float>)( _I * angle) );
field_fd = m_FD_Fields.at(n);
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
field[0][pos[0]][pos[1]][pos[2]] = real(field_fd[0][pos[0]][pos[1]][pos[2]] * exp_jwt)/(float)m_FD_SampleCount;
field[1][pos[0]][pos[1]][pos[2]] = real(field_fd[1][pos[0]][pos[1]][pos[2]] * exp_jwt)/(float)m_FD_SampleCount;
field[2][pos[0]][pos[1]][pos[2]] = real(field_fd[2][pos[0]][pos[1]][pos[2]] * exp_jwt)/(float)m_FD_SampleCount;
}
}
}
stringstream ss;
ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_p=" << std::setw( 3 ) << std::setfill( '0' ) <<(int)(angle * 180 / M_PI) << ".vtk";
ofstream file(ss.str().c_str());
if (file.is_open()==false)
cerr << "ProcessFieldsFD::DumpFDData: can't open file '" << ss.str() << "' for writing... abort! " << endl;
DumpVectorArray2VTK(file,GetFieldNameByType(m_DumpType),field,discLines,numLines,m_precision,string("Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString(), m_Mesh_Type, discLines_scaling);
file.close();
}
{
//dump magnitude to vtk-files
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
field[0][pos[0]][pos[1]][pos[2]] = abs(field_fd[0][pos[0]][pos[1]][pos[2]])/(float)m_FD_SampleCount;
field[1][pos[0]][pos[1]][pos[2]] = abs(field_fd[1][pos[0]][pos[1]][pos[2]])/(float)m_FD_SampleCount;
field[2][pos[0]][pos[1]][pos[2]] = abs(field_fd[2][pos[0]][pos[1]][pos[2]])/(float)m_FD_SampleCount;
}
}
}
stringstream ss;
ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_abs" << ".vtk";
ofstream file(ss.str().c_str());
if (file.is_open()==false)
cerr << "ProcessFieldsFD::DumpFDData: can't open file '" << ss.str() << "' for writing... abort! " << endl;
DumpVectorArray2VTK(file,GetFieldNameByType(m_DumpType),field,discLines,numLines,m_precision,string("Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString(), m_Mesh_Type, discLines_scaling);
file.close();
}
{
//dump phase to vtk-files
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
field[0][pos[0]][pos[1]][pos[2]] = arg(field_fd[0][pos[0]][pos[1]][pos[2]]);
field[1][pos[0]][pos[1]][pos[2]] = arg(field_fd[1][pos[0]][pos[1]][pos[2]]);
field[2][pos[0]][pos[1]][pos[2]] = arg(field_fd[2][pos[0]][pos[1]][pos[2]]);
}
}
}
stringstream ss;
ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n) << "_arg" << ".vtk";
ofstream file(ss.str().c_str());
if (file.is_open()==false)
cerr << "ProcessFieldsFD::DumpFDData: can't open file '" << ss.str() << "' for writing... abort! " << endl;
DumpVectorArray2VTK(file,GetFieldNameByType(m_DumpType),field,discLines,numLines,m_precision,string("Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString(), m_Mesh_Type, discLines_scaling);
file.close();
}
}
Delete_N_3DArray(field,numLines);
return;
}
if (m_fileType==HDF5_FILETYPE)
{
for (size_t n = 0; n<m_FD_Samples.size(); ++n)
{
stringstream ss;
ss << "f" << n;
DumpVectorArray2HDF5(m_filename.c_str(), ss.str(), m_FD_Fields.at(n),numLines,1.0/(float)m_FD_SampleCount,m_FD_Samples.at(n));
}
return;
}
cerr << "ProcessFieldsTD::Process: unknown File-Type" << endl;
}

41
Common/processfields_fd.h Normal file
View File

@ -0,0 +1,41 @@
/*
* Copyright (C) 2010 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 PROCESSFIELDS_FD_H
#define PROCESSFIELDS_FD_H
#include "processfields.h"
class ProcessFieldsFD : public ProcessFields
{
public:
ProcessFieldsFD(Engine_Interface_Base* eng_if);
virtual ~ProcessFieldsFD();
virtual void InitProcess();
virtual int Process();
virtual void PostProcess();
protected:
void DumpFDData();
//! frequency domain field storage
vector<std::complex<float>****> m_FD_Fields;
};
#endif // PROCESSFIELDS_FD_H

View File

@ -69,7 +69,7 @@ int ProcessFieldsTD::Process()
{
stringstream ss;
ss << std::setw( pad_length ) << std::setfill( '0' ) << m_Eng_Interface->GetNumberOfTimesteps();
DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),field,numLines,(0.5+m_Eng_Interface->GetNumberOfTimesteps())*Op->GetTimestep());
DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),field,numLines,m_Eng_Interface->GetTime(m_dualTime));
}
else
cerr << "ProcessFieldsTD::Process: unknown File-Type" << endl;

View File

@ -37,6 +37,7 @@ Processing::Processing(Engine_Interface_Base* eng_if)
m_weight=1;
m_Flush = false;
m_dualMesh = false;
m_dualTime = false;
m_Mesh_Type = CARTESIAN_MESH;
}

View File

@ -97,6 +97,9 @@ public:
//! Dump probe geometry to file
virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh) const;
virtual void SetDualMesh(bool val) {m_dualMesh=val;}
virtual void SetDualTime(bool val) {m_dualTime=val;}
protected:
Engine_Interface_Base* m_Eng_Interface;
const Operator_Base* Op;
@ -130,6 +133,9 @@ protected:
//! define if given coords are on main or dualMesh (default is false)
bool m_dualMesh;
//! define if given processing uses the dual time concept (default is false);
bool m_dualTime;
//! define/store snapped start/stop coords as mesh index
unsigned int start[3];
unsigned int stop[3];

View File

@ -20,7 +20,6 @@
ProcessIntegral::ProcessIntegral(Engine_Interface_Base* eng_if) : Processing(eng_if)
{
m_TimeShift = 0.0;
m_Results=NULL;
}
@ -56,7 +55,7 @@ int ProcessIntegral::Process()
int NrInt = GetNumberOfIntegrals();
double integral = m_Results[0] * m_weight;
double time = m_Eng_Interface->GetTime(false) + m_TimeShift;
double time = m_Eng_Interface->GetTime(m_dualTime);
if (ProcessInterval)
{

View File

@ -53,9 +53,6 @@ public:
protected:
ProcessIntegral(Engine_Interface_Base* eng_if);
//! timeshift to be used in TD and FD data, e.g. 0.5*dT in case of current based parameter
double m_TimeShift;
vector<double> TD_Values;
vector<double_complex> FD_Values;

View File

@ -27,8 +27,6 @@ ProcessModeMatch::ProcessModeMatch(Engine_Interface_Base* eng_if) : ProcessInteg
m_ModeParser[n] = new CSFunctionParser();
m_ModeDist[n] = NULL;
}
m_dualMesh = false;
delete[] m_Results;
m_Results = new double[2];
}
@ -47,11 +45,6 @@ void ProcessModeMatch::InitProcess()
{
if (!Enabled) return;
if (m_ModeFieldType==1)
{
m_TimeShift = Op->GetTimestep()/2.0;
}
if (m_Eng_Interface==NULL)
{
cerr << "ProcessModeMatch::InitProcess: Error, Engine_Interface is NULL, abort mode mathcing..." << endl;

View File

@ -64,7 +64,8 @@ SOURCES += main.cpp \
FDTD/excitation.cpp \
FDTD/operator_cylindermultigrid.cpp \
FDTD/engine_cylindermultigrid.cpp \
FDTD/engine_interface_fdtd.cpp
FDTD/engine_interface_fdtd.cpp \
Common/processfields_fd.cpp
# FDTD/extensions source files
SOURCES += FDTD/extensions/engine_extension.cpp \
@ -116,7 +117,8 @@ HEADERS += tools/ErrorMsg.h \
FDTD/operator_cylindermultigrid.h \
FDTD/engine_cylindermultigrid.h \
tools/aligned_allocator.h \
FDTD/engine_interface_fdtd.h
FDTD/engine_interface_fdtd.h \
Common/processfields_fd.h
# FDTD/extensions header files
HEADERS += FDTD/operator_extension.h \

View File

@ -33,6 +33,7 @@
#include "Common/process_hfield.h"
#include "Common/processmodematch.h"
#include "Common/processfields_td.h"
#include "Common/processfields_fd.h"
#include <sys/time.h>
#include <time.h>
#include <H5Cpp.h> // only for H5get_libversion()
@ -488,6 +489,11 @@ int openEMS::SetupFDTD(const char* file)
}
if (CylinderCoords)
proc->SetMeshType(Processing::CYLINDRICAL_MESH);
if ((pb->GetProbeType()==1) || (pb->GetProbeType()==3) || (pb->GetProbeType()==11))
{
proc->SetDualTime(true);
proc->SetDualMesh(true);
}
proc->SetProcessInterval(Nyquist/m_OverSampling);
proc->AddFrequency(pb->GetFDSamples());
proc->SetName(pb->GetName());
@ -505,15 +511,11 @@ int openEMS::SetupFDTD(const char* file)
vector<CSProperties*> DumpProps = CSX.GetPropertyByType(CSProperties::DUMPBOX);
for (size_t i=0; i<DumpProps.size(); ++i)
{
ProcessFieldsTD* ProcTD = new ProcessFieldsTD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
ProcTD->SetEnable(Enable_Dumps);
ProcTD->SetProcessInterval(Nyquist/m_OverSampling);
ProcessFields* ProcField=NULL;
//only looking for one prim atm
CSPrimitives* prim = DumpProps.at(i)->GetPrimitive(0);
if (prim==NULL)
delete ProcTD;
else
if (prim!=NULL)
{
bool acc;
double bnd[6] = {0,0,0,0,0,0};
@ -527,22 +529,42 @@ int openEMS::SetupFDTD(const char* file)
CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox();
if (db)
{
ProcTD->SetDumpType((ProcessFields::DumpType)db->GetDumpType());
ProcTD->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode());
ProcTD->SetFileType((ProcessFields::FileType)db->GetFileType());
if ((db->GetDumpType()>=0) && (db->GetDumpType()<2))
ProcField = new ProcessFieldsTD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
else if ((db->GetDumpType()>=10) && (db->GetDumpType()<12))
ProcField = new ProcessFieldsFD(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng));
else
cerr << "openEMS::SetupFDTD: unknown dump box type... skipping!" << endl;
if (ProcField)
{
ProcField->SetEnable(Enable_Dumps);
ProcField->SetProcessInterval(Nyquist/m_OverSampling);
if ((db->GetDumpType()==1) || (db->GetDumpType()==11))
{
ProcField->SetDualTime(true);
ProcField->SetDualMesh(true);
}
if ((db->GetDumpType()==10) || (db->GetDumpType()==11))
ProcField->AddFrequency(db->GetFDSamples());
if (db->GetDumpType()>=10)
ProcField->SetDumpType((ProcessFields::DumpType)(db->GetDumpType()-10));
else
ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType());
ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode());
ProcField->SetFileType((ProcessFields::FileType)db->GetFileType());
if (CylinderCoords)
ProcTD->SetMeshType(Processing::CYLINDRICAL_MESH);
ProcField->SetMeshType(Processing::CYLINDRICAL_MESH);
for (int n=0; n<3; ++n)
ProcTD->SetSubSampling(db->GetSubSampling(n),n);
ProcTD->SetFilePattern(db->GetName());
ProcTD->SetFileName(db->GetName());
ProcTD->DefineStartStopCoord(start,stop);
ProcTD->InitProcess();
PA->AddProcessing(ProcTD);
ProcField->SetSubSampling(db->GetSubSampling(n),n);
ProcField->SetFilePattern(db->GetName());
ProcField->SetFileName(db->GetName());
ProcField->DefineStartStopCoord(start,stop);
ProcField->InitProcess();
PA->AddProcessing(ProcField);
prim->SetPrimitiveUsed(true);
}
else
delete ProcTD;
}
}
}