diff --git a/Common/operator_base.h b/Common/operator_base.h index e83f862..41d3844 100644 --- a/Common/operator_base.h +++ b/Common/operator_base.h @@ -114,6 +114,8 @@ public: //! Get the cell center coordinate usable for material averaging (Warning, may not be the yee cell center) virtual bool GetCellCenterMaterialAvgCoord(const int pos[3], double coord[3]) const = 0; + virtual vector GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type=CSProperties::ANY) const = 0; + //! Set the background material (default is vacuum) virtual void SetBackgroundMaterial(double epsR=0, double mueR=0, double kappa=0, double sigma=0, double density=0); diff --git a/Common/processfields_sar.cpp b/Common/processfields_sar.cpp index 8275b8e..79a9f6d 100644 --- a/Common/processfields_sar.cpp +++ b/Common/processfields_sar.cpp @@ -203,6 +203,7 @@ void ProcessFieldsSAR::DumpFDData() for (pos[1]=0; pos[1] vPrims = Op->GetPrimitivesBoundBox(orig_pos[0], orig_pos[1], -1, CSProperties::MATERIAL); for (pos[2]=0; pos[2]GetCellCenterMaterialAvgCoord(orig_pos, coord); - prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); + prop = CSX->GetPropertyByCoordPriority(coord, vPrims); +// prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); if (prop!=0) { matProp = dynamic_cast(prop); diff --git a/FDTD/extensions/operator_ext_conductingsheet.cpp b/FDTD/extensions/operator_ext_conductingsheet.cpp index 135d333..28c51bb 100644 --- a/FDTD/extensions/operator_ext_conductingsheet.cpp +++ b/FDTD/extensions/operator_ext_conductingsheet.cpp @@ -61,6 +61,7 @@ bool Operator_Ext_ConductingSheet::BuildExtension() { for (pos[1]=0; pos[1] vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL)); for (pos[2]=0; pos[2]GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), false, &cs_sheet); +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), false, &cs_sheet); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, false, &cs_sheet); CSPropConductingSheet* cs_prop = dynamic_cast(prop); if (cs_prop) { diff --git a/FDTD/extensions/operator_ext_cylinder.cpp b/FDTD/extensions/operator_ext_cylinder.cpp index 277cad2..d400bf8 100644 --- a/FDTD/extensions/operator_ext_cylinder.cpp +++ b/FDTD/extensions/operator_ext_cylinder.cpp @@ -58,13 +58,15 @@ bool Operator_Ext_Cylinder::BuildExtension() double inEC[4]; double dT = m_Op->GetTimestep(); pos[0]=0; + vector vPrims_metal = m_Op->GetPrimitivesBoundBox(pos[0], -1, -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL)); for (pos[2]=0; pos[2]GetNumberOfLines(2,true); ++pos[2]) { double C=0; double G=0; + vector vPrims_mat = m_Op->GetPrimitivesBoundBox(pos[0], -1, pos[2], CSProperties::MATERIAL); for (pos[1]=0; pos[1]GetNumberOfLines(1,true)-2; ++pos[1]) { - m_Op_Cyl->Calc_ECPos(2,pos,inEC); + m_Op_Cyl->Calc_ECPos(2,pos,inEC,vPrims_mat); C+=inEC[0]; G+=inEC[1]; } @@ -80,7 +82,7 @@ bool Operator_Ext_Cylinder::BuildExtension() //search for metal on z-axis m_Op_Cyl->GetYeeCoords(2,pos,coord,false); - CSProperties* prop = m_Op->CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL), true); + CSProperties* prop = m_Op->CSX->GetPropertyByCoordPriority(coord, vPrims_metal, true); if (prop) { if (prop->GetType()==CSProperties::METAL) //set to PEC diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/FDTD/extensions/operator_ext_lorentzmaterial.cpp index 7ef355c..9f386fc 100644 --- a/FDTD/extensions/operator_ext_lorentzmaterial.cpp +++ b/FDTD/extensions/operator_ext_lorentzmaterial.cpp @@ -204,6 +204,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() { for (pos[1]=0; pos[1] vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL)); for (pos[2]=0; pos[2]MainOp->SetPos(pos[0],pos[1],pos[2]); @@ -222,7 +223,8 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() if (m_Op->GetVI(n,pos[0],pos[1],pos[2])==0) continue; - CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, true); if ((mat = prop->ToLorentzMaterial())) { w_plasma = mat->GetEpsPlasmaFreqWeighted(order,n,coord) * 2 * PI; @@ -273,7 +275,8 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() if (m_Op->GetIV(n,pos[0],pos[1],pos[2])==0) continue; - CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, true); if ((mat = prop->ToLorentzMaterial())) { w_plasma = mat->GetMuePlasmaFreqWeighted(order,n,coord) * 2 * PI; diff --git a/FDTD/extensions/operator_ext_mur_abc.cpp b/FDTD/extensions/operator_ext_mur_abc.cpp index 19a6389..1df8583 100644 --- a/FDTD/extensions/operator_ext_mur_abc.cpp +++ b/FDTD/extensions/operator_ext_mur_abc.cpp @@ -140,13 +140,20 @@ bool Operator_Ext_Mur_ABC::BuildExtension() else coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) - delta/2 / m_Op->GetGridDelta(); + int posBB[3]; + posBB[m_ny] =pos[m_ny]; + posBB[m_nyPP]=-1; + for (pos[m_nyP]=0; pos[m_nyP] vPrims = m_Op->GetPrimitivesBoundBox(posBB[0], posBB[1], posBB[2], CSProperties::MATERIAL); coord[m_nyP] = m_Op->GetDiscLine(m_nyP,pos[m_nyP]); for (pos[m_nyPP]=0; pos[m_nyPP]GetDiscLine(m_nyPP,pos[m_nyPP]); - CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, CSProperties::MATERIAL, false); +// CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, CSProperties::MATERIAL, false); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, false); if (prop) { CSPropMaterial* mat = prop->ToMaterial(); diff --git a/FDTD/extensions/operator_ext_upml.cpp b/FDTD/extensions/operator_ext_upml.cpp index 7086ae9..756abcd 100644 --- a/FDTD/extensions/operator_ext_upml.cpp +++ b/FDTD/extensions/operator_ext_upml.cpp @@ -391,12 +391,13 @@ bool Operator_Ext_UPML::BuildExtension() for (loc_pos[1]=0; loc_pos[1] vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL); for (loc_pos[2]=0; loc_pos[2]Calc_EffMatPos(n,pos,eff_Mat); + m_Op->Calc_EffMatPos(n,pos,eff_Mat,vPrims); CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i); nP = (n+1)%3; nPP = (n+2)%3; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index ee9fd4a..42f9f74 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -753,12 +753,13 @@ void Operator::DumpMaterial2File(string filename) { for (pos[1]=0; pos[1] vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL); for (pos[2]=0; pos[2] vPrims) const { double EffMat[4]; - Calc_EffMatPos(ny,pos,EffMat); + Calc_EffMatPos(ny,pos,EffMat, vPrims); if (m_epsR) m_epsR[ny][pos[0]][pos[1]][pos[2]] = EffMat[0]; @@ -1240,9 +1241,16 @@ bool Operator::GetCellCenterMaterialAvgCoord(const int pos[], double coord[3]) c return true; } -double Operator::GetMaterial(int ny, const double* coords, int MatType, bool markAsUsed) const +double Operator::GetMaterial(int ny, const double* coords, int MatType, vector vPrims, bool markAsUsed) const { - CSProperties* prop = CSX->GetPropertyByCoordPriority(coords,CSProperties::MATERIAL,markAsUsed); + CSProperties* prop = CSX->GetPropertyByCoordPriority(coords,vPrims,markAsUsed); +// CSProperties* old_prop = CSX->GetPropertyByCoordPriority(coords,CSProperties::MATERIAL,markAsUsed); +// if (old_prop!=prop) +// { +// cerr << "ERROR: Unequal properties!" << endl; +// exit(-1); +// } + CSPropMaterial* mat = dynamic_cast(prop); if (mat) { @@ -1282,7 +1290,7 @@ double Operator::GetMaterial(int ny, const double* coords, int MatType, bool mar } } -bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat) const +bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat, vector vPrims) const { int n=ny; double coord[3]; @@ -1302,8 +1310,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) { A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] += GetMaterial(n, coord, 0)*A_n; - EffMat[1] += GetMaterial(n, coord, 1)*A_n; + EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n; + EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n; area+=A_n; } @@ -1312,8 +1320,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) { A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] += GetMaterial(n, coord, 0)*A_n; - EffMat[1] += GetMaterial(n, coord, 1)*A_n; + EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n; + EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n; area+=A_n; } @@ -1323,8 +1331,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) { A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] += GetMaterial(n, coord, 0)*A_n; - EffMat[1] += GetMaterial(n, coord, 1)*A_n; + EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n; + EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n; area+=A_n; } @@ -1333,8 +1341,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) { A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] += GetMaterial(n, coord, 0)*A_n; - EffMat[1] += GetMaterial(n, coord, 1)*A_n; + EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n; + EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n; area+=A_n; } @@ -1352,8 +1360,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) { delta_ny = GetNodeWidth(n,loc_pos,true); - EffMat[2] += delta_ny / GetMaterial(n, coord, 2); - sigma = GetMaterial(n, coord, 3); + EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims); + sigma = GetMaterial(n, coord, 3, vPrims); if (sigma) EffMat[3] += delta_ny / sigma; else @@ -1366,8 +1374,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) { delta_ny = GetNodeWidth(n,loc_pos,true); - EffMat[2] += delta_ny / GetMaterial(n, coord, 2); - sigma = GetMaterial(n, coord, 3); + EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims); + sigma = GetMaterial(n, coord, 3, vPrims); if (sigma) EffMat[3] += delta_ny / sigma; else @@ -1388,7 +1396,7 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff return true; } -bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat) const +bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat, vector vPrims) const { int n=ny; double coord[3]; @@ -1415,8 +1423,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef shiftCoord[nP] = coord[nP]+deltaP*0.25; shiftCoord[nPP] = coord[nPP]+deltaPP*0.25; A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] = GetMaterial(n, shiftCoord, 0)*A_n; - EffMat[1] = GetMaterial(n, shiftCoord, 1)*A_n; + EffMat[0] = GetMaterial(n, shiftCoord, 0, vPrims)*A_n; + EffMat[1] = GetMaterial(n, shiftCoord, 1, vPrims)*A_n; area+=A_n; //shift up-left @@ -1426,8 +1434,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef --loc_pos[nP]; A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; - EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; + EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n; + EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n; area+=A_n; //shift down-right @@ -1437,8 +1445,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef ++loc_pos[nP]; --loc_pos[nPP]; A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; - EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; + EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n; + EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n; area+=A_n; //shift down-left @@ -1447,8 +1455,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25; --loc_pos[nP]; A_n = GetNodeArea(ny,loc_pos,true); - EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; - EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; + EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n; + EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n; area+=A_n; EffMat[0]*=__EPS0__/area; @@ -1466,8 +1474,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; --loc_pos[n]; double delta_ny = GetNodeWidth(n,loc_pos,true); - EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2); - double sigma = GetMaterial(n, shiftCoord, 3); + EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2, vPrims); + double sigma = GetMaterial(n, shiftCoord, 3, vPrims); if (sigma) EffMat[3] = delta_ny / sigma; else @@ -1480,8 +1488,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; ++loc_pos[n]; delta_ny = GetNodeWidth(n,loc_pos,true); - EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2); - sigma = GetMaterial(n, shiftCoord, 3); + EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2, vPrims); + sigma = GetMaterial(n, shiftCoord, 3, vPrims); if (sigma) EffMat[3] += delta_ny / sigma; else @@ -1502,14 +1510,14 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef return true; } -bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) const +bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, vector vPrims) const { switch (m_MatAverageMethod) { case QuarterCell: - return AverageMatQuarterCell(ny, pos, EffMat); + return AverageMatQuarterCell(ny, pos, EffMat, vPrims); case CentralCell: - return AverageMatCellCenter(ny, pos, EffMat); + return AverageMatCellCenter(ny, pos, EffMat, vPrims); default: cerr << "Operator:: " << __func__ << ": Error, unknown material averaging method... exit" << endl; exit(1); @@ -1717,19 +1725,51 @@ bool Operator::Calc_EC() cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl; return false; } + + MainOp->SetPos(0,0,0); + Calc_EC_Range(0,numLines[0]-1); + return true; +} + +vector Operator::GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type) const +{ + double boundBox[6]; + int BBpos[3] = {posX, posY, posZ}; + for (int n=0;n<3;++n) + { + if (BBpos[n]<0) + { + boundBox[2*n] = this->GetDiscLine(n,0); + boundBox[2*n+1] = this->GetDiscLine(n,numLines[n]-1); + } + else + { + boundBox[2*n] = this->GetDiscLine(n, max(0, BBpos[n]-1)); + boundBox[2*n+1] = this->GetDiscLine(n, min(int(numLines[n])-1, BBpos[n]+1)); + } + } + + vector vPrim = this->CSX->GetPrimitivesByBoundBox(boundBox, true, type); + return vPrim; +} + +void Operator::Calc_EC_Range(unsigned int xStart, unsigned int xStop) +{ +// vector vPrims = this->CSX->GetAllPrimitives(true, CSProperties::MATERIAL); unsigned int ipos; unsigned int pos[3]; double inEC[4]; - for (int n=0; n<3; ++n) + for (pos[0]=xStart; pos[0]<=xStop; ++pos[0]) { - for (pos[2]=0; pos[2] vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL); + for (pos[2]=0; pos[2]GetPos(pos[0],pos[1],pos[2]); + for (int n=0; n<3; ++n) { - Calc_ECPos(n,pos,inEC); - ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); + Calc_ECPos(n,pos,inEC,vPrims); EC_C[n][ipos]=inEC[0]; EC_G[n][ipos]=inEC[1]; EC_L[n][ipos]=inEC[2]; @@ -1738,8 +1778,6 @@ bool Operator::Calc_EC() } } } - - return true; } void Operator::SetTimestepFactor(double factor) @@ -1925,12 +1963,19 @@ void Operator::CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned i { for (pos[1]=0; pos[1] vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL)); for (pos[2]=0; pos[2]GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL), true); + CSProperties* prop = CSX->GetPropertyByCoordPriority(coord, vPrims, true); +// CSProperties* old_prop = CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL), true); +// if (old_prop!=prop) +// { +// cerr << "CalcPEC_Range: " << old_prop << " vs " << prop << endl; +// exit(-1); +// } if (prop) { if (prop->GetType()==CSProperties::METAL) //set to PEC diff --git a/FDTD/operator.h b/FDTD/operator.h index 65aed18..8fc0409 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -179,6 +179,8 @@ public: virtual double CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const; + virtual vector GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type=CSProperties::ANY) const; + protected: //! use New() for creating a new Operator Operator(); @@ -215,7 +217,7 @@ protected: 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; + virtual bool Calc_ECPos(int ny, const unsigned int* pos, double* EC, vector vPrims) const; //! Get the FDTD raw disc delta, needed by Calc_EffMatPos() \sa Calc_EffMatPos /*! @@ -226,15 +228,15 @@ protected: 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; + virtual double GetMaterial(int ny, const double coords[3], int MatType, vector vPrims, bool markAsUsed=true) const; MatAverageMethods m_MatAverageMethod; //! 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 Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, vector vPrims) const; - virtual bool AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat) const; - virtual bool AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat) const; + virtual bool AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat, vector vPrims) const; + virtual bool AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat, vector vPrims) const; //! Calc operator at certain \a pos virtual void Calc_ECOperatorPos(int n, unsigned int* pos); @@ -254,6 +256,7 @@ protected: //EC elements, internal only! virtual void Init_EC(); virtual bool Calc_EC(); + virtual void Calc_EC_Range(unsigned int xStart, unsigned int xStop); FDTD_FLOAT* EC_C[3]; FDTD_FLOAT* EC_G[3]; FDTD_FLOAT* EC_L[3]; diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index 4ca5134..4a53d9b 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -67,14 +67,14 @@ double Operator_Cylinder::GetRawDiscDelta(int ny, const int pos) const return Operator_Multithread::GetRawDiscDelta(ny,pos); } -double Operator_Cylinder::GetMaterial(int ny, const double* coords, int MatType, bool markAsUsed) const +double Operator_Cylinder::GetMaterial(int ny, const double* coords, int MatType, vector vPrims, 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] vPrims, bool markAsUsed=true) const; virtual int CalcECOperator( DebugFlags debugFlags = None ); diff --git a/FDTD/operator_cylindermultigrid.cpp b/FDTD/operator_cylindermultigrid.cpp index c61dd15..88f9f82 100644 --- a/FDTD/operator_cylindermultigrid.cpp +++ b/FDTD/operator_cylindermultigrid.cpp @@ -212,18 +212,17 @@ void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, ve void Operator_CylinderMultiGrid::FillMissingDataStorage() { unsigned int pos[3]; - double EffMat[4]; - for (int ny=0; ny<3; ++ny) { for (pos[0]=0; pos[0] vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL); for (pos[2]=0; pos[2]m_CalcEC_Start->wait(); - unsigned int ipos; - unsigned int pos[3]; - double inEC[4]; - for (int n=0; n<3; ++n) - { - for (pos[0]=m_start; pos[0]<=m_stop; ++pos[0]) - { - for (pos[1]=0; pos[1]numLines[1]; ++pos[1]) - { - for (pos[2]=0; pos[2]numLines[2]; ++pos[2]) - { - m_OpPtr->Calc_ECPos(n,pos,inEC); - ipos = m_OpPtr->MainOp->GetPos(pos[0],pos[1],pos[2]); - m_OpPtr->EC_C[n][ipos]=inEC[0]; - m_OpPtr->EC_G[n][ipos]=inEC[1]; - m_OpPtr->EC_L[n][ipos]=inEC[2]; - m_OpPtr->EC_R[n][ipos]=inEC[3]; - } - } - } - } + m_OpPtr->Calc_EC_Range(m_start,m_stop); m_OpPtr->m_CalcEC_Stop->wait(); //************** calculate EC (Calc_EC) ***********************//