From 8d5043bd447af42b82bf566b668c286d7120215d Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 18 Feb 2013 10:38:55 +0100 Subject: [PATCH] operator: alternative material averaging method If constant cell material is activated, material probing is performed only in the center of a primary cell. This should improve and simplify SAR calculation if all materials are assumed as constant within a primary YEE cell. Usage from Matlab/Octave: FDTD = InitFDTD('CellConstantMaterial',1); Signed-off-by: Thorsten Liebig --- FDTD/operator.cpp | 21 +++++++++++---------- FDTD/operator.h | 5 +++++ matlab/InitFDTD.m | 2 ++ openems.cpp | 9 +++++++++ 4 files changed, 27 insertions(+), 10 deletions(-) diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 9af1c2e..7b56921 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -95,6 +95,7 @@ void Operator::Init() m_Exc = 0; m_TimeStepFactor = 1; + m_MatCellShiftFaktor = 0.25; } void Operator::Delete() @@ -1220,8 +1221,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //******************************* epsilon,kappa averaging *****************************// //shift up-right shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]+deltaP*0.25; - shiftCoord[nPP] = coord[nPP]+deltaPP*0.25; + shiftCoord[nP] = coord[nP]+deltaP*m_MatCellShiftFaktor; + shiftCoord[nPP] = coord[nPP]+deltaPP*m_MatCellShiftFaktor; A_n = GetNodeArea(ny,loc_pos,true); EffMat[0] = GetMaterial(n, shiftCoord, 0)*A_n; EffMat[1] = GetMaterial(n, shiftCoord, 1)*A_n; @@ -1229,8 +1230,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //shift up-left shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]-deltaP_M*0.25; - shiftCoord[nPP] = coord[nPP]+deltaPP*0.25; + shiftCoord[nP] = coord[nP]-deltaP_M*m_MatCellShiftFaktor; + shiftCoord[nPP] = coord[nPP]+deltaPP*m_MatCellShiftFaktor; --loc_pos[nP]; A_n = GetNodeArea(ny,loc_pos,true); @@ -1240,8 +1241,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //shift down-right shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]+deltaP*0.25; - shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25; + shiftCoord[nP] = coord[nP]+deltaP*m_MatCellShiftFaktor; + shiftCoord[nPP] = coord[nPP]-deltaPP_M*m_MatCellShiftFaktor; ++loc_pos[nP]; --loc_pos[nPP]; A_n = GetNodeArea(ny,loc_pos,true); @@ -1251,8 +1252,8 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c //shift down-left shiftCoord[n] = coord[n]+delta*0.5; - shiftCoord[nP] = coord[nP]-deltaP_M*0.25; - shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25; + shiftCoord[nP] = coord[nP]-deltaP_M*m_MatCellShiftFaktor; + shiftCoord[nPP] = coord[nPP]-deltaPP_M*m_MatCellShiftFaktor; --loc_pos[nP]; A_n = GetNodeArea(ny,loc_pos,true); EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; @@ -1269,7 +1270,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c double length=0; //shift down - shiftCoord[n] = coord[n]-delta_M*0.25; + shiftCoord[n] = coord[n]-delta_M*m_MatCellShiftFaktor; shiftCoord[nP] = coord[nP]+deltaP*0.5; shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; --loc_pos[n]; @@ -1283,7 +1284,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c length=delta_ny; //shift up - shiftCoord[n] = coord[n]+delta*0.25; + shiftCoord[n] = coord[n]+delta*m_MatCellShiftFaktor; shiftCoord[nP] = coord[nP]+deltaP*0.5; shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; ++loc_pos[n]; diff --git a/FDTD/operator.h b/FDTD/operator.h index 97137f2..396db84 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -81,6 +81,9 @@ public: virtual void SetTimestepFactor(double factor); bool GetTimestepValid() const {return !m_InvaildTimestep;} + //! Set operator to assume a constant material inside a cell (material probing in the cell center) + void SetCellConstantMaterial() {m_MatCellShiftFaktor=0.5;} + virtual double GetNumberCells() const; virtual unsigned int GetNumberOfNyquistTimesteps() const {return m_Exc->GetNyquistNum();} @@ -204,6 +207,8 @@ protected: //! Get the material at a given coordinate, direction and type from CSX (internal use only) virtual double GetMaterial(int ny, const double coords[3], int MatType, bool markAsUsed=true) const; + double m_MatCellShiftFaktor; + //! Calculate the effective/averaged material properties at the given position and direction ny. virtual bool Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const; diff --git a/matlab/InitFDTD.m b/matlab/InitFDTD.m index 8051446..1a0afed 100644 --- a/matlab/InitFDTD.m +++ b/matlab/InitFDTD.m @@ -12,6 +12,8 @@ function FDTD = InitFDTD(varargin) % CoordSystem: choose coordinate system (0 Cartesian, 1 Cylindrical) % TimeStep: force to use a given timestep (dangerous!) % TimeStepFactor: reduce the timestep by a given factor (>0 to <=1) +% CellConstantMaterial: set to 1 to assume a material is constant inside +% a cell (material probing in cell center) % % examples: % %default init with 1e9 max. timesteps and -50dB end-criteria diff --git a/openems.cpp b/openems.cpp index 1971c16..a34e00a 100644 --- a/openems.cpp +++ b/openems.cpp @@ -647,6 +647,15 @@ int openEMS::SetupFDTD(const char* file) if (SetupOperator(FDTD_Opts)==false) return 2; + int cellConstantMaterial=0; + FDTD_Opts->QueryIntAttribute("CellConstantMaterial",&cellConstantMaterial); + if (cellConstantMaterial==1) + { + FDTD_Op->SetCellConstantMaterial(); + if (g_settings.GetVerboseLevel()>0) + cerr << "Enabling constant cell material assumption." << endl; + } + m_Exc = new Excitation(); FDTD_Op->SetExcitationSignal(m_Exc); FDTD_Op->AddExtension(new Operator_Ext_Excitation(FDTD_Op));