From 2c3ebe5a7da77a59c700f7b92e31d74393246819 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Fri, 7 Jan 2011 16:12:07 +0100 Subject: [PATCH] new dump type: electric current density (J = kappa * E) --- Common/engine_interface_base.h | 2 ++ Common/processfields.cpp | 29 ++++++++++++++++++++++----- Common/processfields.h | 8 +++++++- FDTD/engine_interface_fdtd.cpp | 36 +++++++++++++++++++++++++++------- FDTD/engine_interface_fdtd.h | 4 ++++ FDTD/operator.h | 1 + openems.cpp | 20 +++++++++++++++---- 7 files changed, 83 insertions(+), 17 deletions(-) diff --git a/Common/engine_interface_base.h b/Common/engine_interface_base.h index c568ce6..b4808a9 100644 --- a/Common/engine_interface_base.h +++ b/Common/engine_interface_base.h @@ -50,6 +50,8 @@ public: virtual double* GetEField(const unsigned int* pos, double* out) const =0; //! Get the (interpolated) magnetic field at \p pos. \sa SetInterpolationType virtual double* GetHField(const unsigned int* pos, double* out) const =0; + //! Get the (interpolated) electric current density field at \p pos. \sa SetInterpolationType + virtual double* GetJField(const unsigned int* pos, double* out) const =0; //! Calculate the electric field integral along a given line virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const =0; diff --git a/Common/processfields.cpp b/Common/processfields.cpp index 88663e5..5b13504 100644 --- a/Common/processfields.cpp +++ b/Common/processfields.cpp @@ -63,6 +63,8 @@ string ProcessFields::GetFieldNameByType(DumpType type) return "E-Field"; case H_FIELD_DUMP: return "H-Field"; + case J_FIELD_DUMP: + return "J-Field"; } return "unknown field"; } @@ -539,8 +541,9 @@ FDTD_FLOAT**** ProcessFields::CalcField() double out[3]; //create array FDTD_FLOAT**** field = Create_N_3DArray(numLines); - if (m_DumpType==E_FIELD_DUMP) + switch (m_DumpType) { + case E_FIELD_DUMP: for (unsigned int i=0; iGetJField(pos,out); + field[0][i][j][k] = out[0]; + field[1][i][j][k] = out[1]; + field[2][i][j][k] = out[2]; + } + } + } + return field; } cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl; diff --git a/Common/processfields.h b/Common/processfields.h index fb7f767..344013a 100644 --- a/Common/processfields.h +++ b/Common/processfields.h @@ -29,8 +29,14 @@ public: ProcessFields(Engine_Interface_Base* eng_if); virtual ~ProcessFields(); + //! File type definition. enum FileType { VTK_FILETYPE, HDF5_FILETYPE}; - enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP}; + + //! Dump type definitions. + /*! + Current dump types are electric field (E_FIELD_DUMP), magnetic field (H_FIELD_DUMP) and (conduction) electric current density (kappa*E) (J_FIELD_DUMP). + */ + enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP, J_FIELD_DUMP}; virtual void InitProcess(); diff --git a/FDTD/engine_interface_fdtd.cpp b/FDTD/engine_interface_fdtd.cpp index da87055..24fd651 100644 --- a/FDTD/engine_interface_fdtd.cpp +++ b/FDTD/engine_interface_fdtd.cpp @@ -28,6 +28,16 @@ Engine_Interface_FDTD::~Engine_Interface_FDTD() } double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) const +{ + return GetRawInterpolatedField(pos, out, 0); +} + +double* Engine_Interface_FDTD::GetJField(const unsigned int* pos, double* out) const +{ + return GetRawInterpolatedField(pos, out, 1); +} + +double* Engine_Interface_FDTD::GetRawInterpolatedField(const unsigned int* pos, double* out, int type) const { unsigned int iPos[] = {pos[0],pos[1],pos[2]}; int nP,nPP; @@ -40,7 +50,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c { delta = m_Op->GetEdgeLength(n,pos,false); if (delta) - out[n] = m_Eng->GetVolt(n,pos) / delta; + out[n] = GetRawField(n,pos,type) / delta; else out[n] = 0.0; } @@ -49,7 +59,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c for (int n=0; n<3; ++n) { delta = m_Op->GetEdgeLength(n,iPos); - out[n] = m_Eng->GetVolt(n,iPos); + out[n] = GetRawField(n,iPos); if (delta==0) { out[n]=0; @@ -63,7 +73,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c --iPos[n]; double deltaDown = m_Op->GetEdgeLength(n,iPos); double deltaRel = delta / (delta+deltaDown); - out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetVolt(n,iPos)/deltaDown*deltaRel; + out[n] = out[n]*(1.0-deltaRel)/delta + (double)GetRawField(n,iPos)/deltaDown*deltaRel; ++iPos[n]; } break; @@ -79,19 +89,19 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c } delta = m_Op->GetEdgeLength(n,iPos); if (delta) - out[n]=m_Eng->GetVolt(n,iPos)/delta; + out[n]=GetRawField(n,iPos)/delta; ++iPos[nP]; delta = m_Op->GetEdgeLength(n,iPos); if (delta) - out[n]+=m_Eng->GetVolt(n,iPos)/delta; + out[n]+=GetRawField(n,iPos)/delta; ++iPos[nPP]; delta = m_Op->GetEdgeLength(n,iPos); if (delta) - out[n]+=m_Eng->GetVolt(n,iPos)/delta; + out[n]+=GetRawField(n,iPos)/delta; --iPos[nP]; delta = m_Op->GetEdgeLength(n,iPos); if (delta) - out[n]+=m_Eng->GetVolt(n,iPos)/delta; + out[n]+=GetRawField(n,iPos)/delta; --iPos[nPP]; out[n]/=4; } @@ -176,3 +186,15 @@ double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, con } return result; } + + +double Engine_Interface_FDTD::GetRawField(unsigned int n, const unsigned int* pos, int type) const +{ + double value = m_Eng->GetVolt(n,pos[0],pos[1],pos[2]); + if (type==0) + return value; + if ((type==1) && (m_Op->m_kappa)) + return value*m_Op->m_kappa[n][pos[0]][pos[1]][pos[2]]; + + return 0.0; +} diff --git a/FDTD/engine_interface_fdtd.h b/FDTD/engine_interface_fdtd.h index 0151ece..a725019 100644 --- a/FDTD/engine_interface_fdtd.h +++ b/FDTD/engine_interface_fdtd.h @@ -42,6 +42,7 @@ public: virtual double* GetEField(const unsigned int* pos, double* out) const; virtual double* GetHField(const unsigned int* pos, double* out) const; + virtual double* GetJField(const unsigned int* pos, double* out) const; virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const; @@ -51,6 +52,9 @@ public: protected: Operator* m_Op; Engine* m_Eng; + + double* GetRawInterpolatedField(const unsigned int* pos, double* out, int type = 0) const; + double GetRawField(unsigned int n, const unsigned int* pos, int type = 0) const; }; #endif // ENGINE_INTERFACE_FDTD_H diff --git a/FDTD/operator.h b/FDTD/operator.h index 6db0430..cf8fd95 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -32,6 +32,7 @@ class TiXmlElement; class Operator : public Operator_Base { friend class Engine; + friend class Engine_Interface_FDTD; friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice friend class Operator_Ext_PML_SF_Plane; public: diff --git a/openems.cpp b/openems.cpp index 49ec04e..39d7869 100644 --- a/openems.cpp +++ b/openems.cpp @@ -398,6 +398,9 @@ int openEMS::SetupFDTD(const char* file) if (timestep) FDTD_Op->SetTimestep(timestep); + //save kappa material properties, maybe used for dump + FDTD_Op->SetMaterialStoreFlags(1,true); + Operator::DebugFlags debugFlags = Operator::None; if (DebugMat) debugFlags |= Operator::debugMaterial; @@ -407,6 +410,9 @@ int openEMS::SetupFDTD(const char* file) debugFlags |= Operator::debugPEC; FDTD_Op->CalcECOperator( debugFlags ); + //reset flag for kappa material properties, if no dump-box resets it to true, it will be cleaned up... + FDTD_Op->SetMaterialStoreFlags(1,false); + unsigned int maxTime_TS = (unsigned int)(maxTime/FDTD_Op->GetTimestep()); if ((maxTime_TS>0) && (maxTime_TSToDumpBox(); if (db) { - if ((db->GetDumpType()>=0) && (db->GetDumpType()<2)) + 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)) + 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; @@ -545,13 +551,17 @@ int openEMS::SetupFDTD(const char* file) //make dualMesh the default mesh for h-field dumps, maybe overwritten by interpolation type (node-interpolation) ProcField->SetDualMesh(true); } - if ((db->GetDumpType()==10) || (db->GetDumpType()==11)) - ProcField->AddFrequency(db->GetFDSamples()); if (db->GetDumpType()>=10) + { + ProcField->AddFrequency(db->GetFDSamples()); ProcField->SetDumpType((ProcessFields::DumpType)(db->GetDumpType()-10)); + } else ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType()); + if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12)) && Enable_Dumps ) + FDTD_Op->SetMaterialStoreFlags(1,true); //keep kappa material data (prevent cleanup) + ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode()); ProcField->SetFileType((ProcessFields::FileType)db->GetFileType()); if (CylinderCoords) @@ -573,6 +583,8 @@ int openEMS::SetupFDTD(const char* file) } } + FDTD_Op->CleanupMaterialStorage(); + CSX.WarnUnusedPrimitves(cerr); // dump all boxes (voltage, current, fields, ...)