diff --git a/Common/engine_interface_base.cpp b/Common/engine_interface_base.cpp new file mode 100644 index 0000000..f751fc6 --- /dev/null +++ b/Common/engine_interface_base.cpp @@ -0,0 +1,38 @@ +/* +* 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 . +*/ + +#include "engine_interface_base.h" +#include "string" + +Engine_Interface_Base::Engine_Interface_Base() +{ + m_InterpolType = NO_INTERPOLATION; +} + +std::string Engine_Interface_Base::GetInterpolationNameByType(InterpolationType mode) +{ + switch (mode) + { + case NO_INTERPOLATION: + return std::string("None"); + case NODE_INTERPOLATE: + return std::string("Node"); + case CELL_INTERPOLATE: + return std::string("Cell"); + } + return std::string(); +} diff --git a/Common/engine_interface_base.h b/Common/engine_interface_base.h new file mode 100644 index 0000000..846b608 --- /dev/null +++ b/Common/engine_interface_base.h @@ -0,0 +1,56 @@ +/* +* 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 . +*/ + +#ifndef ENGINE_INTERFACE_BASE_H +#define ENGINE_INTERFACE_BASE_H + +#include "tools/global.h" + +//! This is the abstact base for all Engine Interface classes. +/*! + This is the abstact base for all Engine Interface classes. It will provide unified access to the field information of the corresponding engine. + All processing methods should only access this base class. +*/ +class Engine_Interface_Base +{ +public: + enum InterpolationType { NO_INTERPOLATION, NODE_INTERPOLATE, CELL_INTERPOLATE }; + + //! Set the current interpolation type \sa GetInterpolationType + void SetInterpolationType(InterpolationType type) {m_InterpolType=type;} + //! Set the current interpolation type \sa GetInterpolationType + void SetInterpolationType(int type) {m_InterpolType=(InterpolationType)type;} + //! Get the current interpolation type as string \sa SetInterpolationType GetInterpolationType GetInterpolationNameByType + std::string GetInterpolationTypeString() {return GetInterpolationNameByType(m_InterpolType);} + //! Get the current interpolation type \sa SetInterpolationType + InterpolationType GetInterpolationType() {return m_InterpolType;} + + //! Get the (interpolated) electric field at \p pos. \sa SetInterpolationType + virtual double* GetEField(const unsigned int* pos, double* out) const {UNUSED(pos);return out;} + //! Get the (interpolated) magnetic field at \p pos. \sa SetInterpolationType + virtual double* GetHField(const unsigned int* pos, double* out) const {UNUSED(pos);return out;} + + //! Convert the interpolation type into a string. + static std::string GetInterpolationNameByType(InterpolationType mode); + +protected: + Engine_Interface_Base(); + + InterpolationType m_InterpolType; +}; + +#endif // ENGINE_INTERFACE_BASE_H diff --git a/Common/readme.txt b/Common/readme.txt new file mode 100644 index 0000000..89192f3 --- /dev/null +++ b/Common/readme.txt @@ -0,0 +1,5 @@ +readme for openEMS/Commen + +- This folder contains all classes common for all numerical solver included in openEMS (currently only FDTD) + - Engine-Interface classes + - Commen processing classes diff --git a/FDTD/engine_interface_fdtd.cpp b/FDTD/engine_interface_fdtd.cpp new file mode 100644 index 0000000..e138e1f --- /dev/null +++ b/FDTD/engine_interface_fdtd.cpp @@ -0,0 +1,163 @@ +/* +* 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 . +*/ + +#include "engine_interface_fdtd.h" + +Engine_Interface_FDTD::Engine_Interface_FDTD() : Engine_Interface_Base() +{ + m_Op = NULL; + m_Eng = NULL; +} + +Engine_Interface_FDTD::Engine_Interface_FDTD(Operator* op, Engine* eng) : Engine_Interface_Base() +{ + SetOperator(op); + SetEngine(eng); +} + +Engine_Interface_FDTD::~Engine_Interface_FDTD() +{ +} + +double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) const +{ + unsigned int iPos[] = {pos[0],pos[1],pos[2]}; + int nP,nPP; + double delta; + switch (m_InterpolType) + { + default: + case NO_INTERPOLATION: + for (int n=0;n<3;++n) + { + delta = m_Op->GetEdgeLength(n,pos,false); + if (delta) + out[n] = m_Eng->GetVolt(n,pos) / delta; + else + out[n] = 0.0; + } + break; + case NODE_INTERPOLATE: + for (int n=0;n<3;++n) + { + delta = m_Op->GetMeshDelta(n,iPos); + out[n] = m_Eng->GetVolt(n,iPos); + if (delta==0) + { + out[n]=0; + continue; + } + if (pos[n]==0) + { + out[n] /= (delta * 2.0); //make it consistant with upper PEC boundary + continue; + } + --iPos[n]; + double deltaDown = m_Op->GetMeshDelta(n,iPos); + double deltaRel = delta / (delta+deltaDown); + out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetVolt(n,iPos)/deltaDown*deltaRel; + ++iPos[n]; + } + break; + case CELL_INTERPOLATE: + for (int n=0;n<3;++n) + { + nP = (n+1)%3; + nPP = (n+2)%3; + if ((pos[0]==m_Op->GetNumberOfLines(0)-1) || (pos[1]==m_Op->GetNumberOfLines(1)-1) || (pos[2]==m_Op->GetNumberOfLines(2)-1)) + { + out[n] = 0; //electric field outside the field domain is always zero + continue; + } + delta = m_Op->GetMeshDelta(n,iPos); + if (delta) + out[n]=m_Eng->GetVolt(n,iPos)/delta; + ++iPos[nP]; + delta = m_Op->GetMeshDelta(n,iPos); + if (delta) + out[n]+=m_Eng->GetVolt(n,iPos)/delta; + ++iPos[nPP]; + delta = m_Op->GetMeshDelta(n,iPos); + if (delta) + out[n]+=m_Eng->GetVolt(n,iPos)/delta; + --iPos[nP]; + delta = m_Op->GetMeshDelta(n,iPos); + if (delta) + out[n]+=m_Eng->GetVolt(n,iPos)/delta; + --iPos[nPP]; + out[n]/=4; + } + break; + } + return out; +} + +double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) const +{ + unsigned int iPos[] = {pos[0],pos[1],pos[2]}; + int nP,nPP; + double delta; + switch (m_InterpolType) + { + default: + case NO_INTERPOLATION: + out[0] = m_Eng->GetCurr(0,pos) / m_Op->GetEdgeLength(0,pos,true); + out[1] = m_Eng->GetCurr(1,pos) / m_Op->GetEdgeLength(1,pos,true); + out[2] = m_Eng->GetCurr(2,pos) / m_Op->GetEdgeLength(2,pos,true); + break; + case NODE_INTERPOLATE: + for (int n=0;n<3;++n) + { + nP = (n+1)%3; + nPP = (n+2)%3; + if ((pos[0]==m_Op->GetNumberOfLines(0)-1) || (pos[1]==m_Op->GetNumberOfLines(1)-1) || (pos[2]==m_Op->GetNumberOfLines(2)-1) || (pos[nP]==0) || (pos[nPP]==0)) + { + out[n] = 0; + continue; + } + out[n]=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + --iPos[nP]; + out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + --iPos[nPP]; + out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + ++iPos[nP]; + out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + ++iPos[nPP]; + out[n]/=4; + } + break; + case CELL_INTERPOLATE: + for (int n=0;n<3;++n) + { + delta = m_Op->GetMeshDelta(n,iPos,true); + out[n] = m_Eng->GetCurr(n,iPos); + if ((pos[n]>=m_Op->GetNumberOfLines(n)-1)) + { + out[n] = 0; //magnetic field on the outer boundaries is always zero + continue; + } + ++iPos[n]; + double deltaUp = m_Op->GetMeshDelta(n,iPos,true); + double deltaRel = delta / (delta+deltaUp); + out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetCurr(n,iPos)/deltaUp*deltaRel; + --iPos[n]; + } + break; + } + + return out; +} diff --git a/FDTD/engine_interface_fdtd.h b/FDTD/engine_interface_fdtd.h new file mode 100644 index 0000000..c4e0ad1 --- /dev/null +++ b/FDTD/engine_interface_fdtd.h @@ -0,0 +1,52 @@ +/* +* 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 . +*/ + +#ifndef ENGINE_INTERFACE_FDTD_H +#define ENGINE_INTERFACE_FDTD_H + +#include + +#include "Common/engine_interface_base.h" +#include "operator.h" +#include "engine.h" + +class Engine_Interface_FDTD : public Engine_Interface_Base +{ +public: + Engine_Interface_FDTD(); + Engine_Interface_FDTD(Operator* op, Engine* eng); + virtual ~Engine_Interface_FDTD(); + + //! Set the FDTD operator + virtual void SetOperator(Operator* op) {m_Op=op;} + //! Set the FDTD engine + virtual void SetEngine(Engine* eng) {m_Eng=eng;} + + //! Get the FDTD engine in case direct access is needed. Direct access is not recommended! + const Engine* GetFDTDEngine() const {return m_Eng;} + //! Get the FDTD operator in case direct access is needed. Direct access is not recommended! + const Operator* GetFDTDOperator() const {return m_Op;} + + virtual double* GetEField(const unsigned int* pos, double* out) const; + virtual double* GetHField(const unsigned int* pos, double* out) const; + +protected: + Operator* m_Op; + Engine* m_Eng; +}; + +#endif // ENGINE_INTERFACE_FDTD_H diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp index be2c7a9..5a3c790 100644 --- a/FDTD/processfields.cpp +++ b/FDTD/processfields.cpp @@ -22,7 +22,6 @@ ProcessFields::ProcessFields(Operator* op, Engine* eng) : Processing(op, eng) { - m_DumpMode = NO_INTERPOLATION; m_DumpType = E_FIELD_DUMP; // vtk-file is default m_fileType = VTK_FILETYPE; @@ -32,7 +31,6 @@ ProcessFields::ProcessFields(Operator* op, Engine* eng) : Processing(op, eng) for (int n=0;n<3;++n) { numLines[n]=0; - discDLines[n]=NULL; discLines[n]=NULL; } } @@ -41,8 +39,6 @@ ProcessFields::~ProcessFields() { for (int n=0;n<3;++n) { - delete[] discDLines[n]; - discDLines[n]=NULL; delete[] discLines[n]; discLines[n]=NULL; } @@ -55,49 +51,29 @@ void ProcessFields::InitProcess() string names[] = {Op->GetDirName(0),Op->GetDirName(1),Op->GetDirName(2)}; if (m_fileType==HDF5_FILETYPE) { - unsigned int* NrLines; - double** Lines; - - if (m_DumpMode==NO_INTERPOLATION) - { - NrLines = numLines; - Lines = discLines; - } - else if (m_DumpMode==NODE_INTERPOLATE) - { - NrLines = numLines; - Lines = discLines; - } - else if (m_DumpMode==CELL_INTERPOLATE) - { - NrLines = numDLines; - Lines = discDLines; - } - else - return; - m_filename+= ".h5"; + H5::H5File* file = new H5::H5File( m_filename , H5F_ACC_TRUNC ); H5::Group* group = new H5::Group( file->createGroup( "/Mesh" )); for (int n=0;n<3;++n) { hsize_t dimsf[1]; // dataset dimensions - dimsf[0] = NrLines[n]; + 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[NrLines[n]]; - for (unsigned int i=0;iGetGridDelta(); + array[i] = discLines[n][i] * Op->GetGridDelta(); #endif } //write to dataset @@ -123,21 +99,6 @@ string ProcessFields::GetFieldNameByType(DumpType type) return "unknown field"; } -string ProcessFields::GetInterpolationNameByType(DumpMode mode) -{ - switch (mode) - { - case NO_INTERPOLATION: - return string("Interpolation: None"); - case NODE_INTERPOLATE: - return string("Interpolation: Node"); - case CELL_INTERPOLATE: - return string("Interpolation: Cell"); - } - return string(); -} - - void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) { vector lines; @@ -147,7 +108,8 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) if (m_DumpType == H_FIELD_DUMP) dualMesh = true; - if (m_DumpMode==NO_INTERPOLATION) + Engine_Interface_Base::InterpolationType m_DumpMode = m_Eng_Interface->GetInterpolationType(); + if (m_DumpMode==Engine_Interface_Base::NO_INTERPOLATION) { if (!Op->SnapToMesh(dstart,start,dualMesh)) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl; @@ -177,7 +139,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) discLines[n][i] = lines.at(i); } } - else if (m_DumpMode==NODE_INTERPOLATE) + else if (m_DumpMode==Engine_Interface_Base::NODE_INTERPOLATE) { if (Op->SnapToMesh(dstart,start)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl; if (Op->SnapToMesh(dstop,stop)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl; @@ -191,8 +153,8 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) start[n]=stop[n]; stop[n]=help; } - if (stop[n]==Op->GetNumberOfLines(n)-1) - --stop[n]; +// 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]) @@ -206,7 +168,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) discLines[n][i] = lines.at(i); } } - else if (m_DumpMode==CELL_INTERPOLATE) + else if (m_DumpMode==Engine_Interface_Base::CELL_INTERPOLATE) { if (Op->SnapToMesh(dstart,start,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl; if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl; @@ -227,11 +189,11 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) { lines.push_back(Op->GetDiscLine(n,i,true));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i])); } - numDLines[n] = lines.size(); - delete[] discDLines[n]; - discDLines[n] = new double[numDLines[n]]; - for (unsigned int i=0;iSetInterpolationType(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_DumpMode=NODE_INTERPOLATE;} + void SetDumpMode2Node() {m_Eng_Interface->SetInterpolationType(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_DumpMode=CELL_INTERPOLATE;} + void SetDumpMode2Cell() {m_Eng_Interface->SetInterpolationType(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;} @@ -79,11 +78,9 @@ public: void SetFileType(FileType fileType) {m_fileType=fileType;} static string GetFieldNameByType(DumpType type); - static string GetInterpolationNameByType(DumpMode mode); // virtual void Process(); protected: - DumpMode m_DumpMode; DumpType m_DumpType; string filePattern; FileType m_fileType; @@ -94,9 +91,6 @@ protected: //! dump mesh unsigned int numLines[3]; double* discLines[3]; - //! dual dump mesh - unsigned int numDLines[3]; - double* discDLines[3]; }; #endif // PROCESSFIELDS_H diff --git a/FDTD/processfields_td.cpp b/FDTD/processfields_td.cpp index e867f1b..95968be 100644 --- a/FDTD/processfields_td.cpp +++ b/FDTD/processfields_td.cpp @@ -29,270 +29,23 @@ ProcessFieldsTD::~ProcessFieldsTD() { } -void ProcessFieldsTD::DumpNodeInterpol(string filename) +int ProcessFieldsTD::Process() { -#ifdef OUTPUT_IN_DRAWINGUNITS - double discLines_scaling = 1; -#else - double discLines_scaling = Op->GetGridDelta(); -#endif + if (Enabled==false) return -1; + if (filePattern.empty()) return -1; + if (CheckTimestep()==false) return GetNextInterval(); - if (m_DumpType==H_FIELD_DUMP) + string filename; + + if (m_fileType==VTK_FILETYPE) { - //create array - FDTD_FLOAT**** H_T = Create_N_3DArray(numLines); - unsigned int pos[3] = {start[0],start[1],start[2]}; - unsigned int OpPos[3]; - for (pos[0]=0;pos[0]GetCurr(0,OpPos) / Op->GetMeshDelta(0,OpPos,true); - OpPos[1]++; - H_T[0][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(0,OpPos) / Op->GetMeshDelta(0,OpPos,true); - OpPos[2]++; - H_T[0][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(0,OpPos) / Op->GetMeshDelta(0,OpPos,true); - OpPos[1]--; - H_T[0][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(0,OpPos) / Op->GetMeshDelta(0,OpPos,true); - OpPos[2]--; - H_T[0][pos[0]][pos[1]][pos[2]] /= 4.0; - - //in y - H_T[1][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(1,OpPos) / Op->GetMeshDelta(1,OpPos,true); - OpPos[0]++; - H_T[1][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(1,OpPos) / Op->GetMeshDelta(1,OpPos,true); - OpPos[2]++; - H_T[1][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(1,OpPos) / Op->GetMeshDelta(1,OpPos,true); - OpPos[0]--; - H_T[1][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(1,OpPos) / Op->GetMeshDelta(1,OpPos,true); - OpPos[2]--; - H_T[1][pos[0]][pos[1]][pos[2]] /= 4.0; - - //in z - H_T[2][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(2,OpPos) / Op->GetMeshDelta(2,OpPos,true); - OpPos[1]++; - H_T[2][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(2,OpPos) / Op->GetMeshDelta(2,OpPos,true); - OpPos[0]++; - H_T[2][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(2,OpPos) / Op->GetMeshDelta(2,OpPos,true); - OpPos[1]--; - H_T[2][pos[0]][pos[1]][pos[2]] += Eng->GetCurr(2,OpPos) / Op->GetMeshDelta(2,OpPos,true); - OpPos[0]--; - H_T[2][pos[0]][pos[1]][pos[2]] /= 4.0; - } - } - } - - if (m_fileType==VTK_FILETYPE) - { - ofstream file(filename.c_str()); - if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;}; - DumpVectorArray2VTK(file,string("H-Field"),H_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type, discLines_scaling); - file.close(); - } - else if (m_fileType==HDF5_FILETYPE) - { - stringstream ss; - ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); - DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),H_T,numLines,Eng->GetNumberOfTimesteps()*Op->GetTimestep()); - } - else - cerr << "ProcessFieldsTD::DumpNodeInterpol: unknown File-Type" << endl; - Delete_N_3DArray(H_T,numLines); - H_T = NULL; + stringstream ss; + ss << filePattern << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps() << ".vtk"; + filename = ss.str(); } + else + filename = m_filename; - if (m_DumpType==E_FIELD_DUMP) - { - //create array - FDTD_FLOAT**** E_T = Create_N_3DArray(numLines); - unsigned int pos[3] = {start[0],start[1],start[2]}; - unsigned int OpPos[3]; - unsigned int OpPosUp[3]; - double delta, deltaUp, deltaRel; -// cerr << "processing h-fields... " << endl; - for (pos[0]=0;pos[0]GetMeshDelta(n,OpPos); - ++OpPosUp[n]; - deltaUp = Op->GetMeshDelta(n,OpPos); - deltaRel = delta / (delta+deltaUp); - if (delta*deltaUp) - { - E_T[n][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(n,OpPos)*(1-deltaRel)/delta + Eng->GetVolt(n,OpPosUp)/deltaUp*deltaRel; - } - --OpPosUp[n]; - } - } - } - } - if (m_fileType==VTK_FILETYPE) - { - ofstream file(filename.c_str()); - if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;}; - DumpVectorArray2VTK(file,string("E-Field"),E_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type, discLines_scaling); - file.close(); - } - else if (m_fileType==HDF5_FILETYPE) - { - stringstream ss; - ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); - DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),E_T,numLines,(0.5+Eng->GetNumberOfTimesteps())*Op->GetTimestep()); - } - else - cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl; - Delete_N_3DArray(E_T,numLines); - E_T = NULL; - } -} - -void ProcessFieldsTD::DumpCellInterpol(string filename) -{ -#ifdef OUTPUT_IN_DRAWINGUNITS - double discLines_scaling = 1; -#else - double discLines_scaling = Op->GetGridDelta(); -#endif - - if (m_DumpType==E_FIELD_DUMP) - { - //create array - FDTD_FLOAT**** E_T = Create_N_3DArray(numDLines); - unsigned int pos[3] = {start[0],start[1],start[2]}; - unsigned int OpPos[3]; - double delta; -// cerr << "processing e-fields... " << endl; - for (pos[0]=0;pos[0]GetMeshDelta(0,OpPos); //Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; - if (delta) - { - E_T[0][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(0,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetVolt(0,OpPos[0],OpPos[1]+1,OpPos[2]) + Eng->GetVolt(0,OpPos[0],OpPos[1],OpPos[2]+1) + Eng->GetVolt(0,OpPos[0],OpPos[1]+1,OpPos[2]+1); - E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); - } - //in y - delta = Op->GetMeshDelta(1,OpPos); //Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; - if (delta) - { - E_T[1][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(1,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetVolt(1,OpPos[0]+1,OpPos[1],OpPos[2]) + Eng->GetVolt(1,OpPos[0],OpPos[1],OpPos[2]+1) + Eng->GetVolt(1,OpPos[0]+1,OpPos[1],OpPos[2]+1); - E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); - } - //in z - delta = Op->GetMeshDelta(2,OpPos); //Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; - if (delta) - { - E_T[2][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(2,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetVolt(2,OpPos[0],OpPos[1]+1,OpPos[2]) + Eng->GetVolt(2,OpPos[0]+1,OpPos[1],OpPos[2]) + Eng->GetVolt(2,OpPos[0]+1,OpPos[1]+1,OpPos[2]); - E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); - } - } - } - } - - if (m_fileType==VTK_FILETYPE) - { - ofstream file(filename.c_str()); - if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;}; - DumpVectorArray2VTK(file,string("E-Field"),E_T,discDLines,numDLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type, discLines_scaling); - file.close(); - } - else if (m_fileType==HDF5_FILETYPE) - { - stringstream ss; - ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); - DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),E_T,numDLines,Eng->GetNumberOfTimesteps()*Op->GetTimestep()); - } - else - cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl; - Delete_N_3DArray(E_T,numDLines); - E_T = NULL; - } - - if (m_DumpType==1) - { - //create array - FDTD_FLOAT**** H_T = Create_N_3DArray(numDLines); - unsigned int pos[3] = {start[0],start[1],start[2]}; - unsigned int OpPos[3]; - unsigned int OpPosUp[3]; - double delta, deltaUp, deltaRel; -// cerr << "processing h-fields... " << endl; - for (pos[0]=0;pos[0]GetMeshDelta(n,OpPos,true); - ++OpPosUp[n]; - deltaUp = Op->GetMeshDelta(n,OpPos,true); - deltaRel = delta / (delta+deltaUp); - if (delta*deltaUp) - { - H_T[n][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(n,OpPos)*(1-deltaRel)/delta + Eng->GetCurr(n,OpPosUp)/deltaUp*deltaRel; - } - --OpPosUp[n]; - } - } - } - } - if (m_fileType==VTK_FILETYPE) - { - ofstream file(filename.c_str()); - if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;}; - DumpVectorArray2VTK(file,string("H-Field"),H_T,discDLines,numDLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type, discLines_scaling); - file.close(); - } - else if (m_fileType==HDF5_FILETYPE) - { - stringstream ss; - ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); - DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),H_T,numDLines,(0.5+Eng->GetNumberOfTimesteps())*Op->GetTimestep()); - } - else - cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl; - Delete_N_3DArray(H_T,numDLines); - H_T = NULL; - } -} - -void ProcessFieldsTD::DumpNoInterpol(string filename) -{ #ifdef OUTPUT_IN_DRAWINGUNITS double discLines_scaling = 1; #else @@ -301,125 +54,67 @@ void ProcessFieldsTD::DumpNoInterpol(string filename) unsigned int pos[3]; unsigned int OpPos[3]; - double delta[3]; + double out[3]; + //create array + FDTD_FLOAT**** field = Create_N_3DArray(numLines); if (m_DumpType==E_FIELD_DUMP) { - //create array - FDTD_FLOAT**** E_T = Create_N_3DArray(numLines); for (pos[0]=0;pos[0]GetMeshDelta(0,OpPos);//fabs(Op->MainOp->GetIndexDelta(0,OpPos[0])); for (pos[1]=0;pos[1]GetMeshDelta(1,OpPos);//fabs(Op->MainOp->GetIndexDelta(1,OpPos[1])); for (pos[2]=0;pos[2]GetMeshDelta(2,OpPos);//fabs(Op->MainOp->GetIndexDelta(2,OpPos[2])); - if (delta[0]) - E_T[0][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(0,OpPos[0],OpPos[1],OpPos[2])/delta[0];// /Op->gridDelta; - if (delta[1]) - E_T[1][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(1,OpPos[0],OpPos[1],OpPos[2])/delta[1];// /Op->gridDelta; - if (delta[2]) - E_T[2][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(2,OpPos[0],OpPos[1],OpPos[2])/delta[2];// /Op->gridDelta; + m_Eng_Interface->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]; } } } - if (m_fileType==VTK_FILETYPE) - { - ofstream file(filename.c_str()); - if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;}; - DumpVectorArray2VTK(file,string("E-Field"),E_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type, discLines_scaling); - file.close(); - } - else if (m_fileType==HDF5_FILETYPE) - { - stringstream ss; - ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); - DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),E_T,numLines,Eng->GetNumberOfTimesteps()*Op->GetTimestep()); - } - else - cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl; - - Delete_N_3DArray(E_T,numLines); - E_T = NULL; } if (m_DumpType==H_FIELD_DUMP) { - //create array - FDTD_FLOAT**** H_T = Create_N_3DArray(numLines); for (pos[0]=0;pos[0]GetMeshDelta(0,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(0,OpPos[0])); for (pos[1]=0;pos[1]GetMeshDelta(1,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(1,OpPos[1])); for (pos[2]=0;pos[2]GetMeshDelta(2,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(2,OpPos[2])); - //in x - if (delta[0]) - H_T[0][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(0,OpPos[0],OpPos[1],OpPos[2])/delta[0];// /Op->gridDelta; - if (delta[1]) - H_T[1][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(1,OpPos[0],OpPos[1],OpPos[2])/delta[1];// /Op->gridDelta; - if (delta[2]) - H_T[2][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(2,OpPos[0],OpPos[1],OpPos[2])/delta[2];// /Op->gridDelta; + m_Eng_Interface->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]; } } } - if (m_fileType==VTK_FILETYPE) - { - ofstream file(filename.c_str()); - if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;}; - DumpVectorArray2VTK(file,string("H-Field"),H_T,discLines,numLines,m_precision,GetInterpolationNameByType(m_DumpMode), m_Mesh_Type, discLines_scaling); - file.close(); - } - else if (m_fileType==HDF5_FILETYPE) - { - stringstream ss; - ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); - DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),H_T,numLines,(0.5+Eng->GetNumberOfTimesteps())*Op->GetTimestep()); - } - else - cerr << "ProcessFieldsTD::DumpCellInterpol: unknown File-Type" << endl; - - Delete_N_3DArray(H_T,numLines); - H_T = NULL; } -} - -int ProcessFieldsTD::Process() -{ - if (Enabled==false) return -1; - if (filePattern.empty()) return -1; - if (CheckTimestep()==false) return GetNextInterval(); - stringstream ss; - ss << filePattern << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); if (m_fileType==VTK_FILETYPE) { - ss << ".vtk"; - if (m_DumpMode==NO_INTERPOLATION) - DumpNoInterpol(ss.str()); - if (m_DumpMode==NODE_INTERPOLATE) - DumpNodeInterpol(ss.str()); - if (m_DumpMode==CELL_INTERPOLATE) - DumpCellInterpol(ss.str()); + ofstream file(filename.c_str()); + if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' 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(); } else if (m_fileType==HDF5_FILETYPE) { - if (m_DumpMode==NO_INTERPOLATION) - DumpNoInterpol(m_filename); - if (m_DumpMode==NODE_INTERPOLATE) - DumpNodeInterpol(m_filename); - if (m_DumpMode==CELL_INTERPOLATE) - DumpCellInterpol(m_filename); + stringstream ss; + ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); + DumpVectorArray2HDF5(filename.c_str(),string( ss.str() ),field,numLines,(0.5+Eng->GetNumberOfTimesteps())*Op->GetTimestep()); } + else + cerr << "ProcessFieldsTD::Process: unknown File-Type" << endl; + + Delete_N_3DArray(field,numLines); + + return GetNextInterval(); } diff --git a/FDTD/processfields_td.h b/FDTD/processfields_td.h index 6233762..99f4616 100644 --- a/FDTD/processfields_td.h +++ b/FDTD/processfields_td.h @@ -33,10 +33,6 @@ public: protected: int pad_length; - - void DumpNoInterpol(string filename); - void DumpNodeInterpol(string filename); - void DumpCellInterpol(string filename); }; #endif // PROCESSFIELDS_TD_H diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp index 604c98c..95e6033 100644 --- a/FDTD/processing.cpp +++ b/FDTD/processing.cpp @@ -35,10 +35,12 @@ Processing::Processing(Operator* op, Engine* eng) m_Flush = false; m_dualMesh = false; m_Mesh_Type = CARTESIAN_MESH; + m_Eng_Interface = NULL; } Processing::~Processing() { + SetEngineInterface(NULL); file.close(); } @@ -47,6 +49,12 @@ void Processing::Reset() m_PS_pos=0; } +void Processing::SetEngineInterface(Engine_Interface_Base* eng_if) +{ + delete m_Eng_Interface; + m_Eng_Interface = eng_if; +} + bool Processing::CheckTimestep() { if (m_ProcessSteps.size()>m_PS_pos) diff --git a/FDTD/processing.h b/FDTD/processing.h index b6fe330..cee5c25 100644 --- a/FDTD/processing.h +++ b/FDTD/processing.h @@ -25,15 +25,20 @@ typedef std::complex double_complex; #include #include #include +#include "Common/engine_interface_base.h" #include "operator.h" #include "engine.h" + class Processing { public: Processing(Operator* op, Engine* eng); virtual ~Processing(); + //! Set the interface to the engine. Each processing needs its own engine interface. This class will take ownership and cleanup the interface on deletion! + void SetEngineInterface(Engine_Interface_Base* eng_if); + enum MeshType { CARTESIAN_MESH, CYLINDRICAL_MESH}; virtual void SetName(string val) {m_Name=val;} @@ -78,6 +83,7 @@ public: virtual void DumpBox2File( string vtkfilenameprefix, bool dualMesh) const; protected: + Engine_Interface_Base* m_Eng_Interface; Operator* Op; Engine* Eng; MeshType m_Mesh_Type; diff --git a/FDTD/processmodematch.cpp b/FDTD/processmodematch.cpp index 5d73c1d..9c2bd7b 100644 --- a/FDTD/processmodematch.cpp +++ b/FDTD/processmodematch.cpp @@ -51,6 +51,14 @@ void ProcessModeMatch::InitProcess() m_TimeShift = Op->GetTimestep()/2.0; } + if (m_Eng_Interface==NULL) + { + cerr << "ProcessModeMatch::InitProcess: Error, Engine_Interface is NULL, abort mode mathcing..." << endl; + Enabled=false; + return; + } + m_Eng_Interface->SetInterpolationType(Engine_Interface_Base::NODE_INTERPOLATE); + int Dump_Dim=0; for (int n=0;n<3;++n) { @@ -191,56 +199,6 @@ void ProcessModeMatch::SetFieldType(int type) cerr << "ProcessModeMatch::SetFieldType: Warning, unknown field type..." << endl; } -double ProcessModeMatch::GetField(int ny, const unsigned int pos[3]) -{ - if (m_ModeFieldType==0) - return GetEField(ny,pos); - if (m_ModeFieldType==1) - return GetHField(ny,pos); - return 0; -} - -double ProcessModeMatch::GetEField(int ny, const unsigned int pos[3]) -{ - if ((pos[ny]==0) || (pos[ny]==Op->GetNumberOfLines(ny)-1)) - return 0.0; - unsigned int DownPos[] = {pos[0],pos[1],pos[2]}; - --DownPos[ny]; - double delta = Op->GetMeshDelta(ny,pos); - double deltaDown = Op->GetMeshDelta(ny,DownPos); - double deltaRel = delta / (delta+deltaDown); - if (delta*deltaDown) - { - return (double)Eng->GetVolt(ny,pos)*(1.0-deltaRel)/delta + (double)Eng->GetVolt(ny,DownPos)/deltaDown*deltaRel; - } - return 0.0; -} - -double ProcessModeMatch::GetHField(int ny, const unsigned int pos[3]) -{ - if ((pos[ny]==0) || (pos[ny]>=Op->GetNumberOfLines(ny)-1)) - return 0.0; - - unsigned int EngPos[] = {pos[0],pos[1],pos[2]}; - - int nyP = (ny+1)%3; - if (pos[nyP] == 0) - return 0.0; - int nyPP = (ny+2)%3; - if (pos[nyPP] == 0) - return 0.0; - - double hfield = Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); - EngPos[nyP]--; - hfield += Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); - EngPos[nyPP]--; - hfield += Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); - EngPos[nyP]++; - hfield += Eng->GetCurr(ny,EngPos) / Op->GetMeshDelta(ny,EngPos,true); - return hfield/4.0; -} - - double* ProcessModeMatch::CalcMultipleIntegrals() { double value = 0; @@ -254,6 +212,8 @@ double* ProcessModeMatch::CalcMultipleIntegrals() unsigned int pos[3] = {0,0,0}; pos[m_ny] = start[m_ny]; + double out[3]={0,0,0}; + for (unsigned int posP = 0;posPGetNodeArea(m_ny,pos,m_dualMesh); + if (m_ModeFieldType==0) + m_Eng_Interface->GetEField(pos,out); + if (m_ModeFieldType==1) + m_Eng_Interface->GetHField(pos,out); for (int n=0;n<2;++n) { - field = GetField((m_ny+n+1)%3,pos); + field = out[(m_ny+n+1)%3]; value += field * m_ModeDist[n][posP][posPP] * area; purity += field*field * area; } diff --git a/openEMS.pro b/openEMS.pro index 76b8124..c556a8b 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -86,7 +86,9 @@ SOURCES += main.cpp \ FDTD/engine_cylindermultigrid.cpp \ FDTD/engine_ext_cylindermultigrid.cpp \ FDTD/operator_ext_upml.cpp \ - FDTD/engine_ext_upml.cpp + FDTD/engine_ext_upml.cpp \ + Common/engine_interface_base.cpp \ + FDTD/engine_interface_fdtd.cpp HEADERS += tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ @@ -130,7 +132,9 @@ HEADERS += tools/ErrorMsg.h \ FDTD/engine_ext_cylindermultigrid.h \ tools/aligned_allocator.h \ FDTD/operator_ext_upml.h \ - FDTD/engine_ext_upml.h + FDTD/engine_ext_upml.h \ + Common/engine_interface_base.h \ + FDTD/engine_interface_fdtd.h QMAKE_CXXFLAGS_RELEASE = -O3 \ -g \ -march=native diff --git a/openems.cpp b/openems.cpp index 3f18fbf..c3fed9a 100644 --- a/openems.cpp +++ b/openems.cpp @@ -25,6 +25,7 @@ #include "FDTD/operator_ext_mur_abc.h" #include "FDTD/operator_ext_pml_sf.h" #include "FDTD/operator_ext_upml.h" +#include "FDTD/engine_interface_fdtd.h" #include "FDTD/processvoltage.h" #include "FDTD/processcurrent.h" #include "FDTD/process_efield.h" @@ -422,6 +423,7 @@ int openEMS::SetupFDTD(const char* file) //*************** setup processing ************// cout << "Setting up processing..." << endl; + unsigned int Nyquist = FDTD_Op->Exc->GetNyquistNum(); PA = new ProcessingArray(Nyquist); @@ -473,6 +475,7 @@ int openEMS::SetupFDTD(const char* file) } if (CylinderCoords) proc->SetMeshType(Processing::CYLINDRICAL_MESH); + proc->SetEngineInterface(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng)); proc->SetProcessInterval(Nyquist/m_OverSampling); proc->AddFrequency(pb->GetFDSamples()); proc->SetName(pb->GetName()); @@ -508,8 +511,9 @@ int openEMS::SetupFDTD(const char* file) CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox(); if (db) { + ProcTD->SetEngineInterface(new Engine_Interface_FDTD(FDTD_Op,FDTD_Eng)); ProcTD->SetDumpType((ProcessFields::DumpType)db->GetDumpType()); - ProcTD->SetDumpMode((ProcessFields::DumpMode)db->GetDumpMode()); + ProcTD->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode()); ProcTD->SetFileType((ProcessFields::FileType)db->GetFileType()); if (CylinderCoords) ProcTD->SetMeshType(Processing::CYLINDRICAL_MESH); diff --git a/openems.h b/openems.h index bb89df4..20e7c07 100644 --- a/openems.h +++ b/openems.h @@ -24,6 +24,7 @@ using namespace std; class Operator; class Engine; +class Engine_Interface_FDTD; class ProcessingArray; class TiXmlElement; diff --git a/tools/global.h b/tools/global.h index 98b5ccb..55ae94e 100644 --- a/tools/global.h +++ b/tools/global.h @@ -48,4 +48,7 @@ protected: extern Global g_settings; +// declare a parameter as unused +#define UNUSED(x) (void)(x); + #endif // GLOBAL_H