diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index 09c011a..5c7561e 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -20,7 +20,7 @@ //! \brief construct an Engine instance //! it's the responsibility of the caller to free the returned pointer -Engine* Engine::createEngine(const Operator* op) +Engine* Engine::New(const Operator* op) { Engine* e = new Engine(op); e->Init(); @@ -52,73 +52,90 @@ void Engine::Reset() curr=NULL; } -bool Engine::IterateTS(unsigned int iterTS) +inline void Engine::UpdateVoltages() { unsigned int pos[3]; - int exc_pos; bool shift[3]; + //voltage updates + for (pos[0]=0;pos[0]numLines[0];++pos[0]) + { + shift[0]=pos[0]; + for (pos[1]=0;pos[1]numLines[1];++pos[1]) + { + shift[1]=pos[1]; + for (pos[2]=0;pos[2]numLines[2];++pos[2]) + { + shift[2]=pos[2]; + //do the updates here + //for x + volt[0][pos[0]][pos[1]][pos[2]] *= Op->vv[0][pos[0]][pos[1]][pos[2]]; + volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-shift[1]][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-shift[2]]); + + //for y + volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]]; + volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-shift[0]][pos[1]][pos[2]]); + + //for z + volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]]; + volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-shift[1]][pos[2]]); + } + } + } +} + +inline void Engine::ApplyVoltageExcite() +{ + int exc_pos; + //soft voltage excitation here (E-field excite) + for (unsigned int n=0;nE_Exc_Count;++n) + { + exc_pos = (int)numTS - (int)Op->E_Exc_delay[n]; + exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->ExciteLength); +// if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl; + volt[Op->E_Exc_dir[n]][Op->E_Exc_index[0][n]][Op->E_Exc_index[1][n]][Op->E_Exc_index[2][n]] += Op->E_Exc_amp[n]*Op->ExciteSignal[exc_pos]; + } +} + +inline void Engine::UpdateCurrents() +{ + unsigned int pos[3]; + for (pos[0]=0;pos[0]numLines[0]-1;++pos[0]) + { + for (pos[1]=0;pos[1]numLines[1]-1;++pos[1]) + { + for (pos[2]=0;pos[2]numLines[2]-1;++pos[2]) + { + //do the updates here + //for x + curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]]; + curr[0][pos[0]][pos[1]][pos[2]] += Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]][pos[1]+1][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] + volt[1][pos[0]][pos[1]][pos[2]+1]); + + //for y + curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]]; + curr[1][pos[0]][pos[1]][pos[2]] += Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]+1] - volt[2][pos[0]][pos[1]][pos[2]] + volt[2][pos[0]+1][pos[1]][pos[2]]); + + //for z + curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]]; + curr[2][pos[0]][pos[1]][pos[2]] += Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]+1][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]] + volt[0][pos[0]][pos[1]+1][pos[2]]); + } + } + } +} + +inline void Engine::ApplyCurrentExcite() +{ + int exc_pos; +} + +bool Engine::IterateTS(unsigned int iterTS) +{ for (unsigned int iter=0;iternumLines[0];++pos[0]) - { - shift[0]=pos[0]; - for (pos[1]=0;pos[1]numLines[1];++pos[1]) - { - shift[1]=pos[1]; - for (pos[2]=0;pos[2]numLines[2];++pos[2]) - { - shift[2]=pos[2]; - //do the updates here - //for x - volt[0][pos[0]][pos[1]][pos[2]] *= Op->vv[0][pos[0]][pos[1]][pos[2]]; - volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-shift[1]][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-shift[2]]); - - //for y - volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]]; - volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-shift[0]][pos[1]][pos[2]]); - - //for x - volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]]; - volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-shift[1]][pos[2]]); - } - } - } - - //soft voltage excitation here (E-field excite) - for (unsigned int n=0;nE_Exc_Count;++n) - { - exc_pos = (int)numTS - (int)Op->E_Exc_delay[n]; - exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->ExciteLength); -// if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl; - volt[Op->E_Exc_dir[n]][Op->E_Exc_index[0][n]][Op->E_Exc_index[1][n]][Op->E_Exc_index[2][n]] += Op->E_Exc_amp[n]*Op->ExciteSignal[exc_pos]; - } - - //current updates - for (pos[0]=0;pos[0]numLines[0]-1;++pos[0]) - { - for (pos[1]=0;pos[1]numLines[1]-1;++pos[1]) - { - for (pos[2]=0;pos[2]numLines[2]-1;++pos[2]) - { - //do the updates here - //for x - curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]]; - curr[0][pos[0]][pos[1]][pos[2]] += Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]][pos[1]+1][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] + volt[1][pos[0]][pos[1]][pos[2]+1]); - - //for y - curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]]; - curr[1][pos[0]][pos[1]][pos[2]] += Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]+1] - volt[2][pos[0]][pos[1]][pos[2]] + volt[2][pos[0]+1][pos[1]][pos[2]]); - - //for x - curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]]; - curr[2][pos[0]][pos[1]][pos[2]] += Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]+1][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]] + volt[0][pos[0]][pos[1]+1][pos[2]]); - } - } - } - - //soft current excitation here (H-field excite) + UpdateVoltages(); + ApplyVoltageExcite(); + UpdateCurrents(); + ApplyCurrentExcite(); ++numTS; } return true; diff --git a/FDTD/engine.h b/FDTD/engine.h index 7068ac2..230b2ad 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -23,7 +23,7 @@ class Engine { public: - static Engine* createEngine(const Operator* op); + static Engine* New(const Operator* op); virtual ~Engine(); virtual void Init(); @@ -41,6 +41,11 @@ protected: Engine(const Operator* op); const Operator* Op; + virtual inline void UpdateVoltages(); + virtual inline void ApplyVoltageExcite(); + virtual inline void UpdateCurrents(); + virtual inline void ApplyCurrentExcite(); + FDTD_FLOAT**** volt; FDTD_FLOAT**** curr; unsigned int numTS; diff --git a/FDTD/engine_cylinder.cpp b/FDTD/engine_cylinder.cpp new file mode 100644 index 0000000..2ccf96b --- /dev/null +++ b/FDTD/engine_cylinder.cpp @@ -0,0 +1,44 @@ +/* +* 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_cylinder.h" + +Engine_Cylinder* Engine_Cylinder::New(const Operator_Cylinder* op) +{ + Engine_Cylinder* e = new Engine_Cylinder(op); + e->Init(); + return e; +} + +Engine_Cylinder::Engine_Cylinder(const Operator_Cylinder* op) : Engine(op) +{ +} + +Engine_Cylinder::~Engine_Cylinder() +{ +} + +void Engine_Cylinder::Init() +{ + cerr << "Engine_Cylinder::Init()" << endl; + Engine::Init(); +} + +void Engine_Cylinder::Reset() +{ + Engine::Reset(); +} diff --git a/FDTD/engine_cylinder.h b/FDTD/engine_cylinder.h new file mode 100644 index 0000000..af4d122 --- /dev/null +++ b/FDTD/engine_cylinder.h @@ -0,0 +1,37 @@ +/* +* 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_CYLINDER_H +#define ENGINE_CYLINDER_H + +#include "engine.h" +#include "operator_cylinder.h" + +class Engine_Cylinder : public Engine +{ +public: + static Engine_Cylinder* New(const Operator_Cylinder* op); + + virtual void Init(); + virtual void Reset(); + +protected: + Engine_Cylinder(const Operator_Cylinder* op); + ~Engine_Cylinder(); +}; + +#endif // ENGINE_CYLINDER_H diff --git a/FDTD/engine_multithread.cpp b/FDTD/engine_multithread.cpp index 42ab92a..8d140a2 100644 --- a/FDTD/engine_multithread.cpp +++ b/FDTD/engine_multithread.cpp @@ -34,7 +34,7 @@ //! \brief construct an Engine_Multithread instance //! it's the responsibility of the caller to free the returned pointer -Engine_Multithread* Engine_Multithread::createEngine(const Operator* op, unsigned int numThreads) +Engine_Multithread* Engine_Multithread::New(const Operator* op, unsigned int numThreads) { Engine_Multithread* e = new Engine_Multithread(op); e->setNumThreads( numThreads ); diff --git a/FDTD/engine_multithread.h b/FDTD/engine_multithread.h index ef3623d..55df164 100644 --- a/FDTD/engine_multithread.h +++ b/FDTD/engine_multithread.h @@ -73,7 +73,7 @@ class Engine_Multithread : public Engine friend class NS_Engine_Multithread::thread; friend class NS_Engine_Multithread::thread_e_excitation; public: - static Engine_Multithread* createEngine(const Operator* op, unsigned int numThreads = 0); + static Engine_Multithread* New(const Operator* op, unsigned int numThreads = 0); virtual ~Engine_Multithread(); virtual void setNumThreads( unsigned int numThreads ); diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index f3686d9..cf1ebac 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -21,14 +21,20 @@ #include "tools/array_ops.h" #include "fparser.hh" +Operator* Operator::New() +{ + Operator* op = new Operator(); + op->Init(); + return op; +} + Operator::Operator() { - Operator::Init(); } Operator::~Operator() { - Operator::Reset(); + Reset(); } void Operator::Init() @@ -86,7 +92,7 @@ void Operator::Reset() delete[] EC_R[n]; } - Operator::Init(); + Init(); } unsigned int Operator::CalcNyquistNum(double fmax) @@ -449,6 +455,16 @@ void Operator::InitOperator() ii = Create_N_3DArray(numLines); } +inline void Operator::Calc_ECOperatorPos(int n, unsigned int* pos) +{ + unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]); + vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); + vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); + + ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); + iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); +} + int Operator::CalcECOperator() { if (Calc_EC()==0) @@ -469,12 +485,7 @@ int Operator::CalcECOperator() { for (pos[2]=0;pos[2]SetPos(pos[0],pos[1],pos[2]); - vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); - vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); - - ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); - iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); + Calc_ECOperatorPos(n,pos); } } } diff --git a/FDTD/operator.h b/FDTD/operator.h index 142a0d1..fd4a6eb 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -28,7 +28,8 @@ class Operator { public: - Operator(); + //! Create a new operator + static Operator* New(); virtual ~Operator(); virtual bool SetGeometryCSX(ContinuousStructure* geo); @@ -69,6 +70,9 @@ public: bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false); protected: + //! use New() for creating a new Operator + Operator(); + virtual void Init(); virtual void InitOperator(); @@ -92,10 +96,13 @@ protected: double dT; //FDTD timestep! unsigned int m_nyquistTS; + //! Calc operator at certain pos + virtual inline void Calc_ECOperatorPos(int n, unsigned int* pos); + //EC elements, internal only! - bool Calc_EC(); - bool Calc_ECPos(int n, unsigned int* pos, double* inEC); - bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat); + virtual bool Calc_EC(); + virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC); + virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat); double* EC_C[3]; double* EC_G[3]; double* EC_L[3]; diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp new file mode 100644 index 0000000..51ad60b --- /dev/null +++ b/FDTD/operator_cylinder.cpp @@ -0,0 +1,376 @@ +/* +* 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 "operator_cylinder.h" + +Operator_Cylinder* Operator_Cylinder::New() +{ + Operator_Cylinder* op = new Operator_Cylinder(); + op->Init(); + return op; +} + +Operator_Cylinder::Operator_Cylinder() +{ +} + +Operator_Cylinder::~Operator_Cylinder() +{ + Operator::Reset(); +} + +void Operator_Cylinder::Init() +{ + CC_closedAlpha = false; + CC_R0_included = false; + Operator::Init(); +} + +void Operator_Cylinder::Reset() +{ + Operator::Reset(); +} + +bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) +{ + if (Operator::SetGeometryCSX(geo)==false) return false; + + double minmaxA = fabs(discLines[1][numLines[1]-1]-discLines[1][0]); +// cerr << minmaxA -2*PI << " < " << (2*PI)/10/numLines[1] << endl; + if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width + { + CC_closedAlpha = true; + --numLines[1]; + cout << "Operator_Cylinder::SetGeometryCSX: Alpha is a full 2*PI => closed Cylinder..." << endl; + cerr << "Operator_Cylinder::SetGeometryCSX: closed cylinder not yet implemented... exit... " << endl; exit(1); + } + else if (minmaxA>2*PI) + {cerr << "Operator_Cylinder::SetGeometryCSX: Alpha Max-Min must not be larger than 2*PI!!!" << endl; Reset(); return false;} + else CC_closedAlpha=false; + + if (discLines[0][0]<0) + {cerr << "Operator_Cylinder::SetGeometryCSX: r<0 not allowed in Cylinder Coordinates!!!" << endl; Reset(); return false;} + else if (discLines[0][0]==0.0) + { + cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl; + cerr << "Operator_Cylinder::SetGeometryCSX: r=0 included not yet implemented... exit... " << endl; exit(1); + CC_R0_included=true; + } + + return true; +} + +inline void Operator_Cylinder::Calc_ECOperatorPos(int n, unsigned int* pos) +{ + unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]); + if (EC_C[n][i]>0) + { + vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); + vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); + } + else + { + vv[n][pos[0]][pos[1]][pos[2]] = 0; + vi[n][pos[0]][pos[1]][pos[2]] = 0; + } + + if (EC_L[n][i]>0) + { + ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); + iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); + } + else + { + ii[n][pos[0]][pos[1]][pos[2]] = 0; + iv[n][pos[0]][pos[1]][pos[2]] = 0; + } +} + +void Operator_Cylinder::ApplyElectricBC(bool* dirs) +{ + if (dirs==NULL) return; + if (CC_closedAlpha) + { + dirs[2]=0;dirs[3]=0; //no PEC in alpha directions... + } + if (CC_R0_included) + { + dirs[2]=0; //no PEC in r_min directions... + } + + Operator::ApplyElectricBC(dirs); +} + +void Operator_Cylinder::ApplyMagneticBC(bool* dirs) +{ + if (dirs==NULL) return; + if (CC_closedAlpha) + { + dirs[2]=0;dirs[3]=0; //no PMC in alpha directions... + } + if (CC_R0_included) + { + dirs[2]=0; //no PMC in r_min directions... + } + Operator::ApplyMagneticBC(dirs); +} + +bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) +{ + double coord[3]; + double shiftCoord[3]; + int nP = (n+1)%3; + int nPP = (n+2)%3; + coord[0] = discLines[0][pos[0]]; + coord[1] = discLines[1][pos[1]]; + coord[2] = discLines[2][pos[2]]; + unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); + double delta=MainOp->GetIndexDelta(n,pos[n]); + double deltaP=MainOp->GetIndexDelta(nP,pos[nP]); + double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]); + double delta_M=MainOp->GetIndexDelta(n,pos[n]-1); + double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1); + double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1); + double geom_factor,A_n; + + //******************************* 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; + CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL); + switch (n) + { + case 0: + geom_factor = fabs(deltaPP*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25; + break; + case 1: + geom_factor = fabs(deltaP*deltaPP/(delta*coord[0]))*0.25; + break; + case 2: + geom_factor = fabs(deltaPP/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0))*0.25; + break; + } + geom_factor*=gridDelta; + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + inEC[0] = mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__; + inEC[1] = mat->GetKappaWeighted(n,shiftCoord)*geom_factor; + } + else + { + inEC[0] = 1*geom_factor*__EPS0__; + inEC[1] = 0; + } + + //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; + prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL); + switch (n) + { + case 0: + geom_factor = fabs(deltaPP*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25; + break; + case 1: + geom_factor = fabs(deltaP_M*deltaPP/(delta*coord[0]))*0.25; + break; + case 2: + geom_factor = fabs(deltaPP/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0))*0.25; + break; + } + geom_factor*=gridDelta; + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__; + inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*geom_factor; + } + else + { + inEC[0] += 1*geom_factor*__EPS0__; + inEC[1] += 0; + } + + //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; + prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL); + switch (n) + { + case 0: + geom_factor = fabs(deltaPP_M*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25; + break; + case 1: + geom_factor = fabs(deltaP*deltaPP_M/(delta*coord[0]))*0.25; + break; + case 2: + geom_factor = fabs(deltaPP_M/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0))*0.25; + break; + } + geom_factor*=gridDelta; + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__; + inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*geom_factor; + } + else + { + inEC[0] += 1*geom_factor*__EPS0__; + inEC[1] += 0; + } + + //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; + prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL); + switch (n) + { + case 0: + geom_factor = fabs(deltaPP_M*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25; + break; + case 1: + geom_factor = fabs(deltaP_M*deltaPP_M/(delta*coord[0]))*0.25; + break; + case 2: + geom_factor = fabs(deltaPP_M/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0))*0.25; + break; + } + geom_factor*=gridDelta; + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*geom_factor*__EPS0__; + inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*geom_factor; + } + else + { + inEC[0] += 1*geom_factor*__EPS0__; + inEC[1] += 0; + } + + if (CC_R0_included && (n>0) && (coord[0]==0)) + { + inEC[0]=0; + inEC[1]=0; + } + +// if ((n==2) && (pos[1]==0) && (pos[2]==0)) +// cerr << n << " -> " << pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[0] << endl; + + //******************************* mu,sigma averaging *****************************// + //shift down + shiftCoord[n] = coord[n]-delta_M*0.25; + shiftCoord[nP] = coord[nP]+deltaP*0.5; + shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; + prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL); + double delta_n = fabs(delta_M); + if (n==1) + { + delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius + } + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + inEC[2] = delta_n / mat->GetMueWeighted(n,shiftCoord); + if (mat->GetSigma(n)) + inEC[3] = delta_n / mat->GetSigmaWeighted(n,shiftCoord); + else + inEC[3] = 0; + } + else + { + inEC[2] = delta_n; + inEC[3] = 0; + } + //shift up + shiftCoord[n] = coord[n]+delta*0.25; + shiftCoord[nP] = coord[nP]+deltaP*0.5; + shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; + prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL); + delta_n = fabs(delta); + if (n==1) + { + delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius + } + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + inEC[2] += mat->GetMue(n)*delta_n; + if (mat->GetSigmaWeighted(n,shiftCoord)) + inEC[3] += delta_n/mat->GetSigmaWeighted(n,shiftCoord); + else + inEC[3] = 0; + } + else + { + inEC[2] += 1*delta_n; + inEC[3] = 0; + } + + A_n = fabs(deltaP*deltaPP); + if (n==0) //z-direction n==0 -> r; nP==1 -> alpha; nPP==2 -> z + { + A_n = A_n * coord[0]; + } + if (n==2) //z-direction n==2 -> z; nP==0 -> r; nPP==1 -> alpha + { + A_n = fabs(deltaPP) * (pow(coord[0]+fabs(deltaP),2.0) - pow(coord[0],2.0))*0.5; + } + + inEC[2] = gridDelta * A_n * 2 * __MUE0__ / inEC[2]; + if (inEC[3]) inEC[3]=gridDelta * A_n * 2 / inEC[3]; + +// if ((n==0) && (pos[1]==0) && (pos[2]==0)) +// cerr << inEC[2]/(coord[0]) << endl; +// cerr << n << " -> " << pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[2] << endl; + + return true; +} + +bool Operator_Cylinder::Calc_EffMatPos(int n, unsigned int* pos, double* inMat) +{ + return false; //fixme + +// int nP = (n+1)%3; +// int nPP = (n+2)%3; +// +// unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); +// double delta=MainOp->GetIndexDelta(n,pos[n]); +// double deltaP=MainOp->GetIndexDelta(nP,pos[nP]); +// double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]); +// +// double delta_M=MainOp->GetIndexDelta(n,pos[n]-1); +// double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1); +// double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1); +// +// this->Calc_ECPos(n,pos,inMat); +// +// inMat[0] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta; +// inMat[1] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta; +// +// inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta; +// inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta; +// +// return true; +} + diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h new file mode 100644 index 0000000..de4c3ef --- /dev/null +++ b/FDTD/operator_cylinder.h @@ -0,0 +1,49 @@ +/* +* 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 OPERATOR_CYLINDER_H +#define OPERATOR_CYLINDER_H + +#include "operator.h" + +class Operator_Cylinder : public Operator +{ +public: + static Operator_Cylinder* New(); + virtual ~Operator_Cylinder(); + + virtual bool SetGeometryCSX(ContinuousStructure* geo); + + virtual void ApplyElectricBC(bool* dirs); + virtual void ApplyMagneticBC(bool* dirs); + + virtual void Reset(); + +protected: + Operator_Cylinder(); + virtual void Init(); + + virtual inline void Calc_ECOperatorPos(int n, unsigned int* pos); + + bool CC_closedAlpha; + bool CC_R0_included; + + virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC); + virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat); +}; + +#endif // OPERATOR_CYLINDER_H diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp index bac63d6..cb0facd 100644 --- a/FDTD/processfields.cpp +++ b/FDTD/processfields.cpp @@ -81,7 +81,7 @@ void ProcessFields::InitProcess() H5::DataSet dataset = group->createDataSet( names[n].c_str(), datatype, dataspace ); //convert to float... float* array = new float[NrLines[n]]; - for (int i=0;i #include "tools/array_ops.h" -#include "FDTD/operator.h" #include "FDTD/engine.h" +#include "FDTD/engine_cylinder.h" #include "FDTD/engine_multithread.h" #include "FDTD/processvoltage.h" #include "FDTD/processcurrent.h" @@ -43,6 +43,7 @@ openEMS::openEMS() FDTD_Op=NULL; FDTD_Eng=NULL; PA=NULL; + CylinderCoords = false; Enable_Dumps = true; DebugMat = false; DebugOp = false; @@ -190,6 +191,13 @@ int openEMS::SetupFDTD(const char* file) else NrTS = help; + FDTD_Opts->QueryIntAttribute("CylinderCoords",&help); + if (help==1) + { + cout << "Using a cylinder coordinate FDTD..." << endl; + CylinderCoords = true; + } + FDTD_Opts->QueryDoubleAttribute("endCriteria",&endCrit); if (endCrit==0) endCrit=1e-6; @@ -226,8 +234,16 @@ int openEMS::SetupFDTD(const char* file) //*************** setup operator ************// cout << "Create Operator..." << endl; - FDTD_Op = new Operator(); - if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1); + if (CylinderCoords) + { + FDTD_Op = Operator_Cylinder::New(); + } + else + { + FDTD_Op = Operator::New(); + } + + if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(2); FDTD_Op->CalcECOperator(); @@ -251,13 +267,21 @@ int openEMS::SetupFDTD(const char* file) cout << "Creation time for operator: " << difftime(OpDoneTime,startTime) << " s" << endl; //create FDTD engine - switch (m_engine) { - case EngineType_Multithreaded: - FDTD_Eng = Engine_Multithread::createEngine(FDTD_Op,m_engine_numThreads); - break; - default: - FDTD_Eng = Engine::createEngine(FDTD_Op); - break; + if (CylinderCoords) + { + cerr << "openEMS: creating cylinder coordinate FDTD engine..." << endl; + FDTD_Eng = Engine_Cylinder::New((Operator_Cylinder*)FDTD_Op); + } + else + { + switch (m_engine) { + case EngineType_Multithreaded: + FDTD_Eng = Engine_Multithread::New(FDTD_Op,m_engine_numThreads); + break; + default: + FDTD_Eng = Engine::New(FDTD_Op); + break; + } } time_t currTime = time(NULL); diff --git a/openems.h b/openems.h index e93ac75..3b144e2 100644 --- a/openems.h +++ b/openems.h @@ -45,6 +45,9 @@ public: protected: void SetupExcitation(TiXmlElement* Excite); + + bool CylinderCoords; + //! Number of Timesteps unsigned int NrTS; bool Enable_Dumps;