1 /* 2 * Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. 16 */ 17 18 #ifndef OPERATOR_H 19 #define OPERATOR_H 20 21 #include "tools/AdrOp.h" 22 #include "tools/constants.h" 23 #include "excitation.h" 24 #include "Common/operator_base.h" 25 26 class Operator_Extension; 27 class Operator_Ext_Excitation; 28 class Engine; 29 class TiXmlElement; 30 31 //! Basic FDTD-operator 32 class Operator : public Operator_Base 33 { 34 friend class Engine; 35 friend class Engine_Interface_FDTD; 36 friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice 37 friend class Operator_Ext_ConductingSheet; //we need to find a way around this... friend class Operator_Extension only would be nice 38 friend class Operator_Ext_PML_SF_Plane; 39 friend class Operator_Ext_Excitation; 40 friend class Operator_Ext_UPML; 41 friend class Operator_Ext_Cylinder; 42 public: 43 enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4}; 44 45 enum MatAverageMethods {QuarterCell=0, CentralCell=1}; 46 47 //! Create a new operator 48 static Operator* New(); 49 virtual ~Operator(); 50 51 virtual Engine* CreateEngine(); GetEngine()52 virtual Engine* GetEngine() const {return m_Engine;} 53 54 virtual bool SetGeometryCSX(ContinuousStructure* geo); 55 56 virtual int CalcECOperator( DebugFlags debugFlags = None ); 57 58 // the next four functions need to be reimplemented in a derived class GetVV(unsigned int n,unsigned int x,unsigned int y,unsigned int z)59 inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; } GetVI(unsigned int n,unsigned int x,unsigned int y,unsigned int z)60 inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; } GetII(unsigned int n,unsigned int x,unsigned int y,unsigned int z)61 inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[n][x][y][z]; } GetIV(unsigned int n,unsigned int x,unsigned int y,unsigned int z)62 inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[n][x][y][z]; } 63 64 // convenient access functions GetVV(unsigned int n,unsigned int pos[3])65 inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int pos[3] ) const { return GetVV(n,pos[0],pos[1],pos[2]); } GetVI(unsigned int n,unsigned int pos[3])66 inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int pos[3] ) const { return GetVI(n,pos[0],pos[1],pos[2]); } GetII(unsigned int n,unsigned int pos[3])67 inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int pos[3] ) const { return GetII(n,pos[0],pos[1],pos[2]); } GetIV(unsigned int n,unsigned int pos[3])68 inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int pos[3] ) const { return GetIV(n,pos[0],pos[1],pos[2]); } 69 70 // the next four functions need to be reimplemented in a derived class SetVV(unsigned int n,unsigned int x,unsigned int y,unsigned int z,FDTD_FLOAT value)71 inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vv[n][x][y][z] = value; } SetVI(unsigned int n,unsigned int x,unsigned int y,unsigned int z,FDTD_FLOAT value)72 inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vi[n][x][y][z] = value; } SetII(unsigned int n,unsigned int x,unsigned int y,unsigned int z,FDTD_FLOAT value)73 inline virtual void SetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { ii[n][x][y][z] = value; } SetIV(unsigned int n,unsigned int x,unsigned int y,unsigned int z,FDTD_FLOAT value)74 inline virtual void SetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { iv[n][x][y][z] = value; } 75 76 virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries 77 virtual void ApplyMagneticBC(bool* dirs); 78 SetBCSize(int dir,int size)79 virtual void SetBCSize(int dir, int size) {m_BC_Size[dir]=size;} GetBCSize(int dir)80 virtual int GetBCSize(int dir) {return m_BC_Size[dir];} 81 82 //! Set a forced timestep to use by the operator SetTimestep(double ts)83 virtual void SetTimestep(double ts) {dT = ts;} 84 virtual void SetTimestepFactor(double factor); GetTimestepValid()85 bool GetTimestepValid() const {return !m_InvaildTimestep;} 86 87 //! Choose a time step method (0=auto, 1=CFL, 3=Rennings) SetTimeStepMethod(int var)88 void SetTimeStepMethod(int var) {m_TimeStepVar=var;} 89 90 //! Set the material averaging method /sa MatAverageMethods 91 void SetMaterialAvgMethod(MatAverageMethods method); 92 93 //! Set material averaging method to the advanced quarter cell material interpolation (default) SetQuarterCellMaterialAvg()94 void SetQuarterCellMaterialAvg() {m_MatAverageMethod=QuarterCell;} 95 96 //! Set operator to assume a constant material inside a cell (material probing in the cell center) SetCellConstantMaterial()97 void SetCellConstantMaterial() {m_MatAverageMethod=CentralCell;} 98 99 virtual double GetNumberCells() const; 100 GetNumberOfNyquistTimesteps()101 virtual unsigned int GetNumberOfNyquistTimesteps() const {return m_Exc->GetNyquistNum();} 102 103 virtual unsigned int GetNumberOfLines(int ny, bool full=false) const {UNUSED(full);return numLines[ny];} 104 105 virtual void ShowStat() const; 106 virtual void ShowExtStat() const; 107 GetGridDelta()108 virtual double GetGridDelta() const {return gridDelta;} 109 110 //! Get the disc line in \a n direction (in drawing units) 111 virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const; 112 113 //! Get the disc line delta in \a n direction (in drawing units) 114 virtual double GetDiscDelta(int n, unsigned int pos, bool dualMesh=false) const; 115 116 //! Get the coordinates for a given node index and component, according to the yee-algorithm. Returns true if inside the FDTD domain. 117 virtual bool GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool dualMesh) const; 118 119 virtual bool GetNodeCoords(const unsigned int pos[3], double* coords, bool dualMesh=false, CoordinateSystem c_system=UNDEFINED_CS) const; 120 121 //! Get the node width for a given direction \a n and a given mesh position \a pos 122 virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeLength(ny,pos,!dualMesh);} 123 //! Get the node width for a given direction \a n and a given mesh position \a pos 124 virtual double GetNodeWidth(int ny, const int pos[3], bool dualMesh = false) const; 125 126 //! Get the node area for a given direction \a n and a given mesh position \a pos 127 virtual double GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const; 128 //! Get the node area for a given direction \a n and a given mesh position \a pos 129 virtual double GetNodeArea(int ny, const int pos[3], bool dualMesh = false) const; 130 131 //! Get the length of an FDTD edge (unit is meter). 132 virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const; 133 134 //! Get the volume of an FDTD cell 135 virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const; 136 137 //! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos 138 /*! 139 This will return the area around an edge with a given direction, measured at the middle of the edge. 140 In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems. 141 */ 142 virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetNodeArea(ny,pos,dualMesh);} 143 144 virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false, bool fullMesh=false) const; 145 146 //! Snap the given coodinates to mesh indices 147 virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool fullMesh=false, bool* inside=NULL) const; 148 149 //! Snap a given box to the FDTD mesh 150 virtual int SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false, int SnapMethod=0, bool* bStartIn=NULL, bool* bStopIn=NULL) const; 151 152 //! Snap a given line to the operator mesh 153 /*! 154 \param[in] start coorindate of the line 155 \param[in] stop coorindate of the line 156 \param[out] uiStart the snapped line-start coorindate index 157 \param[out] uiStop the snapped line-stop coorindate index 158 \param[in] dualMesh snap to main or dual mesh (default is main mesh) 159 \return returns a status, 0 = success, 1 = start outside, 2 = stop outside, 3 = both outside 160 */ 161 virtual int SnapLine2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh=false, bool fullMesh=false) const; 162 163 virtual void AddExtension(Operator_Extension* op_ext); 164 virtual void DeleteExtension(Operator_Extension* op_ext); GetNumberOfExtentions()165 virtual size_t GetNumberOfExtentions() const {return m_Op_exts.size();} GetExtension(size_t index)166 virtual Operator_Extension* GetExtension(size_t index) const {return m_Op_exts.at(index);} 167 168 virtual void CleanupMaterialStorage(); 169 170 virtual double GetDiscMaterial(int type, int ny, const unsigned int pos[3]) const; 171 172 //! Get the cell center coordinate usable for material averaging (Warning, may not be the yee cell center) 173 virtual bool GetCellCenterMaterialAvgCoord(const int pos[3], double coord[3]) const; 174 175 virtual void SetExcitationSignal(Excitation* exc); GetExcitationSignal()176 virtual Excitation* GetExcitationSignal() const {return m_Exc;} 177 178 Operator_Ext_Excitation* GetExcitationExtension() const; 179 180 virtual double CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const; 181 182 virtual vector<CSPrimitives*> GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type=CSProperties::ANY) const; 183 184 protected: 185 //! use New() for creating a new Operator 186 Operator(); 187 188 virtual void Init(); 189 void Delete(); 190 virtual void Reset(); 191 virtual void InitOperator(); 192 virtual void InitDataStorage(); 193 194 virtual bool SetupCSXGrid(CSRectGrid* grid); 195 196 virtual Grid_Path FindPath(double start[], double stop[]); 197 198 // debug 199 virtual void DumpOperator2File(string filename); 200 virtual void DumpMaterial2File(string filename); 201 virtual void DumpPEC2File( string filename, unsigned int *range = NULL ); 202 203 unsigned int m_Nr_PEC[3]; //count PEC edges 204 virtual bool CalcPEC(); 205 virtual void CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned int* counter); //internal to CalcPEC 206 virtual void CalcPEC_Curves(); //internal to CalcPEC 207 208 //Calc timestep only internal use 209 int m_TimeStepVar; 210 double m_TimeStepFactor; 211 virtual double CalcTimestep(); 212 double opt_dT; 213 bool m_InvaildTimestep; 214 string m_Used_TS_Name; 215 216 double CalcTimestep_Var1(); 217 double CalcTimestep_Var3(); 218 219 //! Calculate the FDTD equivalent circuit parameter for the given position and direction ny. \sa Calc_EffMat_Pos 220 virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC, vector<CSPrimitives *> vPrims) const; 221 222 //! Get the FDTD raw disc delta, needed by Calc_EffMatPos() \sa Calc_EffMatPos 223 /*! 224 Get the raw disc delta for a given position and direction. 225 The result will be positive if a disc delta inside the simulation domain is requested. 226 The result will be the negative value of the first or last disc delta respectivly if the position is outside the field domain. 227 */ 228 virtual double GetRawDiscDelta(int ny, const int pos) const; 229 230 //! Get the material at a given coordinate, direction and type from CSX (internal use only) 231 virtual double GetMaterial(int ny, const double coords[3], int MatType, vector<CSPrimitives*> vPrims, bool markAsUsed=true) const; 232 233 MatAverageMethods m_MatAverageMethod; 234 235 //! Calculate the effective/averaged material properties at the given position and direction ny. 236 virtual bool Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives*> vPrims) const; 237 238 virtual bool AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives*> vPrims) const; 239 virtual bool AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives*> vPrims) const; 240 241 //! Calc operator at certain \a pos 242 virtual void Calc_ECOperatorPos(int n, unsigned int* pos); 243 244 //! Calculate and setup lumped elements 245 virtual bool Calc_LumpedElements(); 246 247 //! Store the size of the applied boundary conditions 248 int m_BC_Size[6]; 249 250 //store material properties for post-processing 251 float**** m_epsR; 252 float**** m_kappa; 253 float**** m_mueR; 254 float**** m_sigma; 255 256 //EC elements, internal only! 257 virtual void Init_EC(); 258 virtual bool Calc_EC(); 259 virtual void Calc_EC_Range(unsigned int xStart, unsigned int xStop); 260 FDTD_FLOAT* EC_C[3]; 261 FDTD_FLOAT* EC_G[3]; 262 FDTD_FLOAT* EC_L[3]; 263 FDTD_FLOAT* EC_R[3]; 264 265 AdrOp* MainOp; 266 267 vector<Operator_Extension*> m_Op_exts; 268 269 Engine* m_Engine; 270 271 // excitation classes 272 Excitation* m_Exc; // excitation time signal class 273 // Operator_Ext_Excitation* m_Op_Ext_Exc; // excitation extension 274 275 // engine/post-proc needs access 276 public: 277 //EC operator 278 FDTD_FLOAT**** vv; //calc new voltage from old voltage 279 FDTD_FLOAT**** vi; //calc new voltage from old current 280 FDTD_FLOAT**** ii; //calc new current from old current 281 FDTD_FLOAT**** iv; //calc new current from old voltage 282 }; 283 284 inline Operator::DebugFlags operator|( Operator::DebugFlags a, Operator::DebugFlags b ) { return static_cast<Operator::DebugFlags>(static_cast<int>(a) | static_cast<int>(b)); } 285 inline Operator::DebugFlags& operator|=( Operator::DebugFlags& a, const Operator::DebugFlags& b ) { return a = a | b; } 286 287 #endif // OPERATOR_H 288