From ff9d362d74c3c2a15dc744d009b2b23d27ea3f57 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 27 Dec 2010 21:23:51 +0100 Subject: [PATCH] ProcessFields: prepare fields dumps at arbitrary positions --- Common/processfields.cpp | 199 +++++++++++++++------------------------ Common/processfields.h | 21 +++-- openems.cpp | 3 - 3 files changed, 88 insertions(+), 135 deletions(-) diff --git a/Common/processfields.cpp b/Common/processfields.cpp index 94d35f7..169b945 100644 --- a/Common/processfields.cpp +++ b/Common/processfields.cpp @@ -27,12 +27,14 @@ ProcessFields::ProcessFields(Engine_Interface_Base* eng_if) : Processing(eng_if) // vtk-file is default m_fileType = VTK_FILETYPE; SetSubSampling(1); + m_SampleType = NONE; SetPrecision(6); m_dualTime = false; for (int n=0; n<3; ++n) { numLines[n]=0; + posLines[n]=NULL; discLines[n]=NULL; } } @@ -41,6 +43,8 @@ ProcessFields::~ProcessFields() { for (int n=0; n<3; ++n) { + delete[] posLines[n]; + posLines[n]=NULL; delete[] discLines[n]; discLines[n]=NULL; } @@ -58,116 +62,34 @@ string ProcessFields::GetFieldNameByType(DumpType type) return "unknown field"; } +void ProcessFields::InitProcess() +{ + CalcMeshPos(); +} + +void ProcessFields::SetDumpMode(Engine_Interface_Base::InterpolationType mode) +{ + m_Eng_Interface->SetInterpolationType(mode); + if (mode==Engine_Interface_Base::CELL_INTERPOLATE) + m_dualMesh=true; + else + m_dualMesh=false; +} + void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) { - vector lines; + Processing::DefineStartStopCoord(dstart,dstop); - // determine mesh type - bool dualMesh = false; - if (m_DumpType == H_FIELD_DUMP) - dualMesh = true; - - Engine_Interface_Base::InterpolationType m_DumpMode = m_Eng_Interface->GetInterpolationType(); - if (m_DumpMode==Engine_Interface_Base::NO_INTERPOLATION) + // normalize order of start and stop + for (int n=0; n<3; ++n) { - if (!Op->SnapToMesh(dstart,start,dualMesh)) - cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl; - if (!Op->SnapToMesh(dstop,stop,dualMesh)) - cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl; - - for (int n=0; n<3; ++n) + if (start[n]>stop[n]) { - // normalize order of start and stop - if (start[n]>stop[n]) - { - unsigned int help = start[n]; - start[n]=stop[n]; - stop[n]=help; - } - - // construct new discLines - lines.clear(); - for (unsigned int i=start[n]; i<=stop[n]; i+=subSample[n]) - { - lines.push_back(Op->GetDiscLine(n,i,dualMesh)); - } - numLines[n] = lines.size(); - delete[] discLines[n]; - discLines[n] = new double[numLines[n]]; - for (unsigned int i=0; iSnapToMesh(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; - - //create mesh - for (int n=0; n<3; ++n) - { - if (start[n]>stop[n]) - { - unsigned int help = start[n]; - start[n]=stop[n]; - stop[n]=help; - } -// if (stop[n]==Op->GetNumberOfLines(n)-1) -// --stop[n]; -// cerr << "start " << start[n] << "stop " << stop[n]; - lines.clear(); - for (unsigned int i=start[n]; i<=stop[n]; i+=subSample[n]) - { - lines.push_back(Op->GetDiscLine(n,i));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i])); - } - numLines[n] = lines.size(); - delete[] discLines[n]; - discLines[n] = new double[numLines[n]]; - for (unsigned int i=0; iSnapToMesh(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; - - //create dual mesh - for (int n=0; n<3; ++n) - { - // cerr << "start " << start[n] << "stop " << stop[n]; - if (start[n]>stop[n]) - { - unsigned int help = start[n]; - start[n]=stop[n]; - stop[n]=help; - } - ++stop[n]; - lines.clear(); - for (unsigned int i=start[n]; iGetDiscLine(n,i,true));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i])); - } - numLines[n] = lines.size(); - delete[] discLines[n]; - discLines[n] = new double[numLines[n]]; - for (unsigned int i=0; iGetDiscLine( 0, start[0], dualMesh ) << "," - << Op->GetDiscLine( 1, start[1], dualMesh ) << "," << Op->GetDiscLine( 2, start[2], dualMesh ) << ") -> (" - << Op->GetDiscLine( 0, stop[0], dualMesh ) << ","<< Op->GetDiscLine( 1, stop[1], dualMesh ) << "," - << Op->GetDiscLine( 2, stop[2], dualMesh ) << ")"; - cerr << " [" << start[0] << "," << start[1] << "," << start[2] << "] -> [" - << stop[0] << "," << stop[1] << "," << stop[2] << "]" << endl; - } - } double ProcessFields::CalcTotalEnergy() const @@ -210,6 +132,34 @@ void ProcessFields::SetSubSampling(unsigned int subSampleRate, int dir) subSample[2]=subSampleRate; } else subSample[dir]=subSampleRate; + m_SampleType = SUBSAMPLE; +} + +void ProcessFields::CalcMeshPos() +{ + if ((m_SampleType==SUBSAMPLE) || (m_SampleType==NONE)) + { + vector tmp_pos; + + for (int n=0; n<3; ++n) + { + // construct new discLines + tmp_pos.clear(); + for (unsigned int i=start[n]; i<=stop[n]; i+=subSample[n]) + tmp_pos.push_back(i); + + numLines[n] = tmp_pos.size(); + delete[] discLines[n]; + discLines[n] = new double[numLines[n]]; + delete[] posLines[n]; + posLines[n] = new unsigned int[numLines[n]]; + for (unsigned int i=0; iGetDiscLine(n,tmp_pos.at(i),m_dualMesh); + } + } + } } void ProcessFields::WriteVTKHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT, double discLines_scaling) @@ -535,25 +485,25 @@ bool ProcessFields::DumpVectorArray2HDF5(string filename, string groupName, stri FDTD_FLOAT**** ProcessFields::CalcField() { unsigned int pos[3]; - unsigned int OpPos[3]; double out[3]; //create array FDTD_FLOAT**** field = Create_N_3DArray(numLines); if (m_DumpType==E_FIELD_DUMP) { - for (pos[0]=0; pos[0]GetEField(OpPos,out); - field[0][pos[0]][pos[1]][pos[2]] = out[0]; - field[1][pos[0]][pos[1]][pos[2]] = out[1]; - field[2][pos[0]][pos[1]][pos[2]] = out[2]; + pos[2]=posLines[2][k]; + + m_Eng_Interface->GetEField(pos,out); + field[0][i][j][k] = out[0]; + field[1][i][j][k] = out[1]; + field[2][i][j][k] = out[2]; } } } @@ -562,19 +512,20 @@ FDTD_FLOAT**** ProcessFields::CalcField() if (m_DumpType==H_FIELD_DUMP) { - for (pos[0]=0; pos[0]GetHField(OpPos,out); - field[0][pos[0]][pos[1]][pos[2]] = out[0]; - field[1][pos[0]][pos[1]][pos[2]] = out[1]; - field[2][pos[0]][pos[1]][pos[2]] = out[2]; + pos[2]=posLines[2][k]; + + m_Eng_Interface->GetHField(pos,out); + field[0][i][j][k] = out[0]; + field[1][i][j][k] = out[1]; + field[2][i][j][k] = out[2]; } } } diff --git a/Common/processfields.h b/Common/processfields.h index db6127e..ecd5afa 100644 --- a/Common/processfields.h +++ b/Common/processfields.h @@ -32,6 +32,8 @@ public: enum FileType { VTK_FILETYPE, HDF5_FILETYPE}; enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP}; + virtual void InitProcess(); + virtual void DefineStartStopCoord(double* dstart, double* dstop); //! Define a field dump sub sampling rate for a given direction (default: \a dir = -1 means all directions) @@ -44,11 +46,11 @@ public: void SetFileName(string fn) {m_filename=fn;} //! Define the Dump-Mode - void SetDumpMode(Engine_Interface_Base::InterpolationType mode) {m_Eng_Interface->SetInterpolationType(mode);} + void SetDumpMode(Engine_Interface_Base::InterpolationType mode); //! This methode will dump all fields on a main cell node using 2 E-field and 4 H-fields per direction. - void SetDumpMode2Node() {m_Eng_Interface->SetInterpolationType(Engine_Interface_Base::NODE_INTERPOLATE);} + void SetDumpMode2Node() {SetDumpMode(Engine_Interface_Base::NODE_INTERPOLATE);} //! 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. - void SetDumpMode2Cell() {m_Eng_Interface->SetInterpolationType(Engine_Interface_Base::CELL_INTERPOLATE);} + void SetDumpMode2Cell() {SetDumpMode(Engine_Interface_Base::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(DumpType type) {m_DumpType=type;} @@ -84,18 +86,21 @@ public: static string GetFieldNameByType(DumpType type); -// virtual void Process(); protected: DumpType m_DumpType; string filePattern; FileType m_fileType; - //! field dump sub-sampling + enum SampleType {NONE, SUBSAMPLE} m_SampleType; + virtual void CalcMeshPos(); + + //! field dump sub-sampling (if enabled) unsigned int subSample[3]; - //! dump mesh - unsigned int numLines[3]; - double* discLines[3]; + //! dump mesh information + unsigned int numLines[3]; //number of lines to dump + unsigned int* posLines[3]; //grid positions to dump + double* discLines[3]; //mesh disc lines to dump //! Calculate and return the defined field. Caller has to cleanup the array. FDTD_FLOAT**** CalcField(); diff --git a/openems.cpp b/openems.cpp index b880131..cfe589d 100644 --- a/openems.cpp +++ b/openems.cpp @@ -540,10 +540,7 @@ int openEMS::SetupFDTD(const char* file) 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)