diff --git a/FDTD/extensions/engine_ext_cylinder.cpp b/FDTD/extensions/engine_ext_cylinder.cpp index d880fd4..bf93b3b 100644 --- a/FDTD/extensions/engine_ext_cylinder.cpp +++ b/FDTD/extensions/engine_ext_cylinder.cpp @@ -61,14 +61,13 @@ void Engine_Ext_Cylinder::Apply2Voltages() //close alpha unsigned int pos[3]; - // copy voltages from last alpha-plane to first - unsigned int last_A_Line = numLines[1]-1; + // copy tangential voltages from last alpha-plane to first + unsigned int last_A_Line = numLines[1]-2; for (pos[0]=0; pos[0]SetVolt(0,pos[0],0,pos[2], m_Eng->GetVolt(0,pos[0],last_A_Line,pos[2]) ); - m_Eng->SetVolt(1,pos[0],0,pos[2], m_Eng->GetVolt(1,pos[0],last_A_Line,pos[2]) ); m_Eng->SetVolt(2,pos[0],0,pos[2], m_Eng->GetVolt(2,pos[0],last_A_Line,pos[2]) ); } } @@ -80,14 +79,13 @@ void Engine_Ext_Cylinder::Apply2Current() //close alpha unsigned int pos[3]; - // copy currents from first alpha-plane to last + // copy tangential currents from first alpha-plane to last for (pos[0]=0; pos[0]SetCurr(0,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(0,pos[0],0,pos[2]) ); - m_Eng->SetCurr(1,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(1,pos[0],0,pos[2]) ); m_Eng->SetCurr(2,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(2,pos[0],0,pos[2]) ); } } diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 7f0f09c..275292d 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -56,6 +56,8 @@ Engine* Operator::CreateEngine() const void Operator::Init() { + CSX = NULL; + Operator_Base::Init(); vv=NULL; @@ -177,6 +179,15 @@ double Operator::GetEdgeLength(int n, const unsigned int* pos, bool dualMesh) co } } +double Operator::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const +{ + if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) ) + return 0.0; + + unsigned int uiPos[]={pos[0],pos[1],pos[2]}; + return GetNodeWidth(ny, uiPos, dualMesh); +} + double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const { int nyP = (ny+1)%3; @@ -184,6 +195,15 @@ double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) c return GetNodeWidth(nyP,pos,dualMesh) * GetNodeWidth(nyPP,pos,dualMesh); } +double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const +{ + if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) ) + return 0.0; + + unsigned int uiPos[]={pos[0],pos[1],pos[2]}; + return GetNodeArea(ny, uiPos, dualMesh); +} + bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool lower, bool* inside) const { bool ok=true; @@ -577,6 +597,35 @@ void Operator::DumpMaterial2File(string filename) cout << " done!" << endl; } + bool Operator::SetupCSXGrid(CSRectGrid* grid) + { + for (int n=0; n<3; ++n) + { + discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true); + if (numLines[n]<3) + { + cerr << "CartOperator::SetupCSXGrid: you need at least 3 disc-lines in every direction (3D!)!!!" << endl; + Reset(); + return false; + } + } + MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]); + MainOp->SetGrid(discLines[0],discLines[1],discLines[2]); + if (grid->GetDeltaUnit()<=0) + { + cerr << "CartOperator::SetupCSXGrid: grid delta unit must not be <=0 !!!" << endl; + Reset(); + return false; + } + else gridDelta=grid->GetDeltaUnit(); + MainOp->SetGridDelta(1); + MainOp->AddCellAdrOp(); + + //delete the grid clone... + delete grid; + return true; + } + bool Operator::SetGeometryCSX(ContinuousStructure* geo) { if (geo==NULL) return false; @@ -584,29 +633,8 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo) CSX = geo; CSRectGrid* grid=CSX->GetGrid(); - for (int n=0; n<3; ++n) - { - discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true); - if (n==1) - if (numLines[n]<3) - { - cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl; - Reset(); - return false; - } - } - MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]); - MainOp->SetGrid(discLines[0],discLines[1],discLines[2]); - if (grid->GetDeltaUnit()<=0) - { - cerr << "CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!!" << endl; - Reset(); - return false; - } - else gridDelta=grid->GetDeltaUnit(); - MainOp->SetGridDelta(1); - MainOp->AddCellAdrOp(); - return true; + + return SetupCSXGrid(CSRectGrid::Clone(grid)); } void Operator::InitOperator() @@ -982,6 +1010,44 @@ double Operator::GetRawDiscDelta(int ny, const int pos) const return (discLines[ny][pos+1] - discLines[ny][pos]); } +double Operator::GetMaterial(int ny, const double* coords, int MatType, bool markAsUsed) const +{ + CSProperties* prop = CSX->GetPropertyByCoordPriority(coords,CSProperties::MATERIAL,markAsUsed); + CSPropMaterial* mat = dynamic_cast(prop); + if (mat) + { + switch (MatType) + { + case 0: + return mat->GetEpsilonWeighted(ny,coords); + case 1: + return mat->GetKappaWeighted(ny,coords); + case 2: + return mat->GetMueWeighted(ny,coords); + case 3: + return mat->GetSigmaWeighted(ny,coords); + default: + cerr << "Operator::GetMaterial: Error: unknown material type" << endl; + return 0; + } + } + + switch (MatType) + { + case 0: + return 1; + case 1: + return 0; + case 2: + return 1; + case 3: + return 0; + default: + cerr << "Operator::GetMaterial: Error: unknown material type" << endl; + return 0; + } +} + bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const { int n=ny; @@ -1008,23 +1074,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c shiftCoord[n] = coord[n]+delta*0.5; shiftCoord[nP] = coord[nP]+deltaP*0.25; shiftCoord[nPP] = coord[nPP]+deltaPP*0.25; - A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true); -// { -// cerr << ny << " " << pos[0] << " " << pos[1] << " " << pos[2] << ": " << A_n << endl; -// exit(0); -// } - CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true); - if (prop) - { - CSPropMaterial* mat = prop->ToMaterial(); - EffMat[0] = mat->GetEpsilonWeighted(n,shiftCoord)*A_n; - EffMat[1] = mat->GetKappaWeighted(n,shiftCoord)*A_n; - } - else - { - EffMat[0] = 1*A_n; - EffMat[1] = 0; - } + A_n = GetNodeArea(ny,loc_pos,true); + EffMat[0] = GetMaterial(n, shiftCoord, 0)*A_n; + EffMat[1] = GetMaterial(n, shiftCoord, 1)*A_n; area+=A_n; //shift up-left @@ -1033,20 +1085,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c shiftCoord[nPP] = coord[nPP]+deltaPP*0.25; --loc_pos[nP]; - A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true); -// cerr << A_n << endl; - prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true); - if (prop) - { - CSPropMaterial* mat = prop->ToMaterial(); - EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n; - EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n; - } - else - { - EffMat[0] += 1*A_n; - EffMat[1] += 0; - } + A_n = GetNodeArea(ny,loc_pos,true); + EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; + EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; area+=A_n; //shift down-right @@ -1055,19 +1096,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25; ++loc_pos[nP]; --loc_pos[nPP]; - A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true); - prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true); - if (prop) - { - CSPropMaterial* mat = prop->ToMaterial(); - EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n; - EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n; - } - else - { - EffMat[0] += 1*A_n; - EffMat[1] += 0; - } + A_n = GetNodeArea(ny,loc_pos,true); + EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; + EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; area+=A_n; //shift down-left @@ -1075,19 +1106,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c shiftCoord[nP] = coord[nP]-deltaP_M*0.25; shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25; --loc_pos[nP]; - A_n = GetNodeArea(ny,(unsigned int*)loc_pos,true); - prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true); - if (prop) - { - CSPropMaterial* mat = prop->ToMaterial(); - EffMat[0] += mat->GetEpsilonWeighted(n,shiftCoord)*A_n; - EffMat[1] += mat->GetKappaWeighted(n,shiftCoord)*A_n; - } - else - { - EffMat[0] += 1*A_n; - EffMat[1] += 0; - } + A_n = GetNodeArea(ny,loc_pos,true); + EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; + EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; area+=A_n; EffMat[0]*=__EPS0__/area; @@ -1098,27 +1119,19 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c loc_pos[1]=pos[1]; loc_pos[2]=pos[2]; double length=0; + //shift down shiftCoord[n] = coord[n]-delta_M*0.25; shiftCoord[nP] = coord[nP]+deltaP*0.5; shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; --loc_pos[n]; - double delta_ny = GetNodeWidth(n,(unsigned int*)loc_pos,true); - prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true); - if (prop) - { - CSPropMaterial* mat = prop->ToMaterial(); - EffMat[2] = delta_ny / mat->GetMueWeighted(n,shiftCoord); - if (mat->GetSigmaWeighted(n,shiftCoord)) - EffMat[3] = delta_ny / mat->GetSigmaWeighted(n,shiftCoord); - else - EffMat[3] = 0; - } + double delta_ny = GetNodeWidth(n,loc_pos,true); + EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2); + double sigma = GetMaterial(n, shiftCoord, 3); + if (sigma) + EffMat[3] = delta_ny / sigma; else - { - EffMat[2] = delta_ny; EffMat[3] = 0; - } length=delta_ny; //shift up @@ -1126,22 +1139,13 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c shiftCoord[nP] = coord[nP]+deltaP*0.5; shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; ++loc_pos[n]; - delta_ny = GetNodeWidth(n,(unsigned int*)loc_pos,true); - prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL,true); - if (prop) - { - CSPropMaterial* mat = prop->ToMaterial(); - EffMat[2] += delta_ny / mat->GetMueWeighted(n,shiftCoord); - if (mat->GetSigmaWeighted(n,shiftCoord)) - EffMat[3] += delta_ny/mat->GetSigmaWeighted(n,shiftCoord); - else - EffMat[3] = 0; - } + delta_ny = GetNodeWidth(n,loc_pos,true); + EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2); + sigma = GetMaterial(n, shiftCoord, 3); + if (sigma) + EffMat[3] += delta_ny / sigma; else - { - EffMat[2] += 1*delta_ny; EffMat[3] = 0; - } length+=delta_ny; EffMat[2] = length * __MUE0__ / EffMat[2]; diff --git a/FDTD/operator.h b/FDTD/operator.h index fcb9d7f..20c4df2 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -35,6 +35,8 @@ class Operator : public Operator_Base 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; friend class Operator_Ext_Excitation; + friend class Operator_Ext_UPML; + friend class Operator_Ext_Cylinder; public: enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4}; @@ -48,20 +50,6 @@ public: virtual int CalcECOperator( DebugFlags debugFlags = None ); - //! Calculate the FDTD equivalent circuit parameter for the given position and direction ny. \sa Calc_EffMat_Pos - virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC) const; - - //! Get the FDTD raw disc delta, needed by Calc_EffMatPos() \sa Calc_EffMatPos - /*! - Get the raw disc delta for a given position and direction. - The result will be positive if a disc delta inside the simulation domain is requested. - The result will be the negative value of the first or last disc delta respectivly if the position is outside the field domain. - */ - virtual double GetRawDiscDelta(int ny, const int pos) const; - - //! 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; - virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) {return Exc->setupExcitation(Excite,maxTS);}; virtual void DumpExciationSignals(); @@ -112,9 +100,13 @@ public: //! Get the node width for a given direction \a n and a given mesh position \a pos virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeLength(ny,pos,!dualMesh);} + //! Get the node width for a given direction \a n and a given mesh position \a pos + virtual double GetNodeWidth(int ny, const int pos[3], bool dualMesh = false) const; //! Get the node area for a given direction \a n and a given mesh position \a pos virtual double GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const; + //! Get the node area for a given direction \a n and a given mesh position \a pos + virtual double GetNodeArea(int ny, const int pos[3], bool dualMesh = false) const; //! Get the length of an FDTD edge (unit is meter). virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const; @@ -147,6 +139,8 @@ protected: virtual void InitDataStorage(); virtual void InitExcitation(); + virtual bool SetupCSXGrid(CSRectGrid* grid); + struct Grid_Path { vector posPath[3]; @@ -174,6 +168,23 @@ protected: double CalcTimestep_Var1(); double CalcTimestep_Var3(); + //! Calculate the FDTD equivalent circuit parameter for the given position and direction ny. \sa Calc_EffMat_Pos + virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC) const; + + //! Get the FDTD raw disc delta, needed by Calc_EffMatPos() \sa Calc_EffMatPos + /*! + Get the raw disc delta for a given position and direction. + The result will be positive if a disc delta inside the simulation domain is requested. + The result will be the negative value of the first or last disc delta respectivly if the position is outside the field domain. + */ + virtual double GetRawDiscDelta(int ny, const int pos) const; + + //! 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; + + //! 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; + //! Calc operator at certain \a pos virtual void Calc_ECOperatorPos(int n, unsigned int* pos); diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index c812e70..392c391 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -46,6 +46,27 @@ void Operator_Cylinder::Init() Operator_Multithread::Init(); } +double Operator_Cylinder::GetRawDiscDelta(int ny, const int pos) const +{ + if (CC_closedAlpha && ny==1 && pos==-1) + { +// cerr << (discLines[1][numLines[1]-2] - discLines[1][numLines[1]-3]) << " vs " << Operator_Multithread::GetRawDiscDelta(ny,pos) << endl; + return (discLines[1][numLines[1]-2] - discLines[1][numLines[1]-3]); + } + + return Operator_Multithread::GetRawDiscDelta(ny,pos); +} + +double Operator_Cylinder::GetMaterial(int ny, const double* coords, int MatType, bool markAsUsed) const +{ + double l_coords[] = {coords[0],coords[1],coords[2]}; + if (CC_closedAlpha && (coords[1]>GetDiscLine(1,0,false)+2*PI)) + l_coords[1]-=2*PI; + if (CC_closedAlpha && (coords[1]GetDiscLine(1,0,false)+2*PI) - { + if (CC_closedAlpha && (coords[1]>GetDiscLine(1,0,false)+2*PI)) coords[1]-=2*PI; - } + if (CC_closedAlpha && (coords[1]=numLines[ny]) return 0.0; @@ -127,6 +156,18 @@ double Operator_Cylinder::GetNodeArea(int ny, const unsigned int pos[3], bool du return Operator_Multithread::GetNodeArea(ny,pos,dualMesh); } +double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) const +{ + if ( (pos[0]<0) || (pos[1]<0 && CC_closedAlpha==false) || (pos[2]<0) ) + return 0.0; + + unsigned int uiPos[]={pos[0],pos[1],pos[2]}; + if (pos[1]<0 && CC_closedAlpha==true) + uiPos[1]+=numLines[1]-2; + + return GetNodeArea(ny, uiPos, dualMesh); +} + double Operator_Cylinder::GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh) const { double length = Operator_Multithread::GetEdgeLength(ny,pos,dualMesh); @@ -143,31 +184,23 @@ double Operator_Cylinder::GetEdgeArea(int ny, const unsigned int pos[3], bool du return GetEdgeLength(1,pos,!dualMesh) * GetEdgeLength(2,pos,!dualMesh); } -bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) +bool Operator_Cylinder::SetupCSXGrid(CSRectGrid* grid) { - if (Operator_Multithread::SetGeometryCSX(geo)==false) return false; + unsigned int alphaNum; + double* alphaLines = NULL; + alphaLines = grid->GetLines(1,alphaLines,alphaNum,true); - double minmaxA = fabs(discLines[1][numLines[1]-1]-discLines[1][0]); - if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width + double minmaxA = fabs(alphaLines[alphaNum-1]-alphaLines[0]); + if (fabs(minmaxA-2*PI) < OPERATOR_CYLINDER_CLOSED_ALPHA_THRESHOLD) { - cout << "Operator_Cylinder::SetGeometryCSX: Alpha is a full 2*PI => closed Cylinder..." << endl; + cout << "Operator_Cylinder::SetupCSXGrid: Alpha is a full 2*PI => closed Cylinder..." << endl; CC_closedAlpha = true; - discLines[1][numLines[1]-1] = discLines[1][0] + 2*PI; - cerr << "Operator_Cylinder::SetGeometryCSX: Warning, not handling the disc-line width and material averaging correctly yet for a closed cylinder..." << endl; - if (MainOp->GetIndexDelta(1,0)-MainOp->GetIndexDelta(1,numLines[1]-2) > (2*PI)/10/numLines[1]) - { - cerr << "Operator_Cylinder::SetGeometryCSX: first and last angle delta must be the same... deviation to large..." << MainOp->GetIndexDelta(1,0) - MainOp->GetIndexDelta(1,numLines[1]-2) << endl; - exit(1); - } - if (MainOp->GetIndexDelta(1,0)-MainOp->GetIndexDelta(1,numLines[1]-2) > 0) - { - cerr << "Operator_Cylinder::SetGeometryCSX: first and last angle delta must be the same... auto correction of deviation: " << MainOp->GetIndexDelta(1,0) - MainOp->GetIndexDelta(1,numLines[1]-2) << endl; - discLines[1][numLines[1]-2] = discLines[1][numLines[1]-1]-MainOp->GetIndexDelta(1,0); - } + grid->SetLine(1,alphaNum-1,2*PI+alphaLines[0]); + grid->AddDiscLine(1,2*PI+alphaLines[1]); } else if (minmaxA>2*PI) { - cerr << "Operator_Cylinder::SetGeometryCSX: Alpha Max-Min must not be larger than 2*PI!!!" << endl; + cerr << "Operator_Cylinder::SetupCSXGrid: Alpha Max-Min must not be larger than 2*PI!!!" << endl; Reset(); return false; } @@ -176,48 +209,27 @@ bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) CC_closedAlpha=false; } - if (discLines[0][0]<0) + if (grid->GetLine(0,0)<0) { - cerr << "Operator_Cylinder::SetGeometryCSX: r<0 not allowed in Cylinder Coordinates!!!" << endl; + cerr << "Operator_Cylinder::SetupCSXGrid: r<0 not allowed in Cylinder Coordinates!!!" << endl; Reset(); return false; } - else if (discLines[0][0]==0.0) + else if (grid->GetLine(0,0)==0.0) { - cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl; + cout << "Operator_Cylinder::SetupCSXGrid: r=0 included..." << endl; CC_R0_included= true; //also needed for correct ec-calculation } + if (Operator_Multithread::SetupCSXGrid(grid)==false) + return false; + if (CC_closedAlpha || CC_R0_included) this->AddExtension(new Operator_Ext_Cylinder(this)); return true; } -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) - { - // E in alpha direction ( aka volt[1][x][y][z] ) is not defined for r==0 --> always zero... - unsigned int pos[3] = {0,0,0}; - for (pos[1]=0; pos[1] 1e-10) { - cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, mesh has to be homogenous in alpha direction for multi grid engine, violation found at: " << n << endl; + cerr << "Operator_CylinderMultiGrid::SetupCSXGrid: Error, mesh has to be homogenous in alpha direction for multi grid engine, violation found at: " << n << endl; exit(0); } } @@ -88,16 +88,23 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo) if (m_Split_Rad < discLines[0][n]) { m_Split_Pos = n; - cout << "Operator_CylinderMultiGrid::SetGeometryCSX: Found mesh split position @" << m_Split_Pos << endl; + cout << "Operator_CylinderMultiGrid::SetupCSXGrid: Found mesh split position @" << m_Split_Pos << endl; m_Split_Rad = discLines[0][n]; break; } } if ((m_Split_Pos<4) || (m_Split_Pos>numLines[0]-4)) { - cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, split invalid..." << endl; + cerr << "Operator_CylinderMultiGrid::SetupCSXGrid: Error, split invalid..." << endl; return false; } + return true; +} + +bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo) +{ + if (Operator_Cylinder::SetGeometryCSX(geo)==false) + return false; CSRectGrid* grid = geo->GetGrid(); diff --git a/FDTD/operator_cylindermultigrid.h b/FDTD/operator_cylindermultigrid.h index ade2c82..50ce94a 100644 --- a/FDTD/operator_cylindermultigrid.h +++ b/FDTD/operator_cylindermultigrid.h @@ -43,8 +43,6 @@ public: virtual unsigned int GetSplitPos() const {return m_Split_Pos;} - virtual int CalcECOperator( DebugFlags debugFlags = None ); - virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS); virtual void SetBoundaryCondition(int* BCs); @@ -67,6 +65,10 @@ protected: void Delete(); virtual void Reset(); + virtual bool SetupCSXGrid(CSRectGrid* grid); + + virtual int CalcECOperator( DebugFlags debugFlags = None ); + //! The material data storage in the sub-grid area's will not be filled by the base-operator. Check and do this here! void FillMissingDataStorage(); diff --git a/FDTD/operator_multithread.h b/FDTD/operator_multithread.h index eefe188..32103b1 100644 --- a/FDTD/operator_multithread.h +++ b/FDTD/operator_multithread.h @@ -38,8 +38,6 @@ public: static Operator_Multithread* New(unsigned int numThreads = 0); virtual ~Operator_Multithread(); - virtual int CalcECOperator( DebugFlags debugFlags = None ); - virtual void setNumThreads( unsigned int numThreads ); virtual Engine* CreateEngine() const; @@ -55,6 +53,8 @@ protected: unsigned int (*m_Nr_PEC_thread)[3]; //count PEC edges per thread virtual bool CalcPEC(); //this method is using multi-threading + virtual int CalcECOperator( DebugFlags debugFlags = None ); + //Calc_EC barrier boost::barrier* m_CalcEC_Start; boost::barrier* m_CalcEC_Stop; diff --git a/FDTD/operator_sse_compressed.h b/FDTD/operator_sse_compressed.h index 66643ad..8bd5bba 100644 --- a/FDTD/operator_sse_compressed.h +++ b/FDTD/operator_sse_compressed.h @@ -45,8 +45,6 @@ public: virtual Engine* CreateEngine() const; - virtual int CalcECOperator( DebugFlags debugFlags = None ); - inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVV(n,x,y,z);} inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVI(n,x,y,z);} inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetII(n,x,y,z);} @@ -73,6 +71,8 @@ protected: bool CompareOperators(unsigned int pos1[3], unsigned int pos2[3]); + virtual int CalcECOperator( DebugFlags debugFlags = None ); + // engine needs access public: unsigned int*** m_Op_index;