From f62da05c127443cc80e7a749ede2ce23d4606163 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Fri, 27 Apr 2012 16:31:36 +0200 Subject: [PATCH] rewritten SAR calculation --- Common/operator_base.h | 3 + Common/processfields.cpp | 6 +- Common/processfields_sar.cpp | 198 ++++++++++++++++++++++++++--------- Common/processfields_sar.h | 11 +- FDTD/operator.cpp | 9 ++ FDTD/operator.h | 3 + FDTD/operator_cylinder.cpp | 5 + FDTD/operator_cylinder.h | 3 + 8 files changed, 182 insertions(+), 56 deletions(-) diff --git a/Common/operator_base.h b/Common/operator_base.h index ebd8379..114adcf 100644 --- a/Common/operator_base.h +++ b/Common/operator_base.h @@ -69,6 +69,9 @@ public: */ virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const =0; + //! Get the volume of an FDTD cell + virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const =0; + //! Snap the given coodinates to mesh indices, return box dimension virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool* inside=NULL) const =0; diff --git a/Common/processfields.cpp b/Common/processfields.cpp index 2edeeda..5a6c91e 100644 --- a/Common/processfields.cpp +++ b/Common/processfields.cpp @@ -241,7 +241,6 @@ FDTD_FLOAT**** ProcessFields::CalcField() switch (m_DumpType) { case E_FIELD_DUMP: - case SAR_LOCAL_DUMP: for (unsigned int i=0; iGetInterpolationType()!=Engine_Interface_Base::NODE_INTERPOLATE) + if (m_Eng_Interface->GetInterpolationType()!=Engine_Interface_Base::CELL_INTERPOLATE) { - cerr << "ProcessFieldsSAR::InitProcess(): Warning, interpolation type is not supported, resetting to NODE!" << endl; - SetDumpMode2Node(); + cerr << "ProcessFieldsSAR::InitProcess(): Warning, interpolation type is not supported, resetting to CELL!" << endl; + SetDumpMode2Cell(); } ProcessFieldsFD::InitProcess(); -} -double ProcessFieldsSAR::GetKappaDensityRatio(const unsigned int* pos) -{ - - double coord[3] = {discLines[0][pos[0]],discLines[1][pos[1]],discLines[2][pos[2]]}; - unsigned int OpPos[] = {posLines[0][pos[0]],posLines[1][pos[1]],posLines[2][pos[2]]}; - ContinuousStructure* CSX = Op->GetGeometryCSX(); - CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); - - if (prop==0) - return 0.0; - CSPropMaterial* matProp = dynamic_cast(prop); - double density = matProp->GetDensityWeighted(coord); - - if (density==0) - return 0.0; - - double kappa = 0; - double max_kappa = 0; - - for (int n=0;n<3;++n) - { - kappa = Op->GetDiscMaterial(1,n,OpPos); - if (kappa>max_kappa) - max_kappa = kappa; - if (OpPos[n]>0) - { - --OpPos[n]; - kappa = Op->GetDiscMaterial(1,n,OpPos); - if (kappa>max_kappa) - max_kappa = kappa; - ++OpPos[n]; - } - } - - return max_kappa/density; -} - -void ProcessFieldsSAR::DumpFDData() -{ - unsigned int pos[3]; - FDTD_FLOAT*** SAR = Create3DArray(numLines); - std::complex**** field_fd = NULL; - double kdRatio = 0; + if (Enabled==false) return; + //create data structures... for (size_t n = 0; n >(numLines)); + m_J_FD_Fields.push_back(Create_N_3DArray >(numLines)); + } +} + +int ProcessFieldsSAR::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(); + + std::complex**** field_fd = NULL; + unsigned int pos[3]; + double T; + FDTD_FLOAT**** field_td=NULL; + + // calc E-field + m_DumpType = E_FIELD_DUMP; + field_td = CalcField(); + T = m_Eng_Interface->GetTime(m_dualTime); + for (size_t n = 0; n exp_jwt_2_dt = std::exp( (std::complex)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) ); + exp_jwt_2_dt *= 2; // *2 for single-sided spectrum + exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval + field_fd = m_E_FD_Fields.at(n); for (pos[0]=0; pos[0](field_td,numLines); + + // calc J-field + m_DumpType = J_FIELD_DUMP; + field_td = CalcField(); + T = m_Eng_Interface->GetTime(m_dualTime); + for (size_t n = 0; n exp_jwt_2_dt = std::exp( (std::complex)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) ); + exp_jwt_2_dt *= 2; // *2 for single-sided spectrum + exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval + field_fd = m_J_FD_Fields.at(n); + for (pos[0]=0; pos[0](field_td,numLines); + + //reset dump type + m_DumpType = SAR_LOCAL_DUMP; + + ++m_FD_SampleCount; + return GetNextInterval(); +} + +void ProcessFieldsSAR::DumpFDData() +{ + unsigned int pos[3]; + unsigned int orig_pos[3]; + FDTD_FLOAT*** SAR = Create3DArray(numLines); + std::complex**** E_field_fd = NULL; + std::complex**** J_field_fd = NULL; + double coord[3]; + double density; + ContinuousStructure* CSX = Op->GetGeometryCSX(); + CSProperties* prop = NULL; + CSPropMaterial* matProp = NULL; + + double power; + + for (size_t n = 0; nGetDiscLine(0,orig_pos[0],true); + for (pos[1]=0; pos[1]GetDiscLine(1,orig_pos[1],true); + for (pos[2]=0; pos[2]GetDiscLine(2,orig_pos[2],true); + prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); + SAR[pos[0]][pos[1]][pos[2]] = 0.0; + density=0.0; + if (prop!=0) + { + matProp = dynamic_cast(prop); + density = matProp->GetDensityWeighted(coord); + if (density>0) + { + SAR[pos[0]][pos[1]][pos[2]] = abs(E_field_fd[0][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[0][pos[0]][pos[1]][pos[2]]); + SAR[pos[0]][pos[1]][pos[2]] += abs(E_field_fd[1][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[1][pos[0]][pos[1]][pos[2]]); + SAR[pos[0]][pos[1]][pos[2]] += abs(E_field_fd[2][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[2][pos[0]][pos[1]][pos[2]]); + SAR[pos[0]][pos[1]][pos[2]] *= 0.5/density; + } + else + density=0.0; + } + power+=SAR[pos[0]][pos[1]][pos[2]]*Op->GetCellVolume(orig_pos)*density; } } } @@ -126,6 +217,9 @@ void ProcessFieldsSAR::DumpFDData() 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; + float pow[1]={power}; + if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",pow,1)==false) + cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl; } else cerr << "ProcessFieldsSAR::Process: unknown File-Type" << endl; diff --git a/Common/processfields_sar.h b/Common/processfields_sar.h index adc4975..cfa69ef 100644 --- a/Common/processfields_sar.h +++ b/Common/processfields_sar.h @@ -30,10 +30,19 @@ public: virtual void InitProcess(); + virtual int Process(); + + virtual void SetSubSampling(unsigned int subSampleRate, int dir=-1); + + virtual void SetOptResolution(double optRes, int dir=-1); + protected: virtual void DumpFDData(); - double GetKappaDensityRatio(const unsigned int* pos); + //! frequency domain electric field storage + vector****> m_E_FD_Fields; + //! frequency domain current density storage + vector****> m_J_FD_Fields; }; #endif // PROCESSFIELDS_SAR_H diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index ec57b19..fd3d467 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -181,6 +181,14 @@ double Operator::GetEdgeLength(int n, const unsigned int* pos, bool dualMesh) co } } +double Operator::GetCellVolume(const unsigned int pos[3], bool dualMesh) const +{ + double vol=1; + for (int n=0;n<3;++n) + vol*=GetEdgeLength(n,pos,dualMesh); + return vol; +} + double Operator::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const { if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) ) @@ -1338,6 +1346,7 @@ bool Operator::Calc_LumpedElements() EC_C[ny][ipos] = epsilon * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos); if (R>0) EC_G[ny][ipos] = kappa * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos); + if (R==0) //make lumped element a PEC if resistance is zero { SetVV(ny,pos[0],pos[1],pos[2], 0 ); diff --git a/FDTD/operator.h b/FDTD/operator.h index da817eb..278010f 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -111,6 +111,9 @@ public: //! Get the length of an FDTD edge (unit is meter). virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const; + //! Get the volume of an FDTD cell + virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const; + //! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos /*! This will return the area around an edge with a given direction, measured at the middle of the edge. diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index 546b6f1..9cd160a 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -176,6 +176,11 @@ double Operator_Cylinder::GetEdgeLength(int ny, const unsigned int pos[3], bool return length * GetDiscLine(0,pos[0],dualMesh); } +double Operator_Cylinder::GetCellVolume(const unsigned int pos[3], bool dualMesh) const +{ + return GetEdgeArea(2,pos,dualMesh)*GetEdgeLength(2,pos,dualMesh); +} + double Operator_Cylinder::GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh) const { if (ny!=0) diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h index e4c1824..88005c1 100644 --- a/FDTD/operator_cylinder.h +++ b/FDTD/operator_cylinder.h @@ -58,6 +58,9 @@ public: //! Get the length of an FDTD edge, including radius corrected alpha-mesh width. virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const; + //! Get the volume of an FDTD cell + virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const; + //! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos /*! This will return the area around an edge with a given direction, measured at the middle of the edge.