caiyu/openEMS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

operator: speed improvement using a primitive pre-selection

Browse Source 
Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
This commit is contained in:
Thorsten Liebig 2014-10-19 21:59:39 +02:00
parent 650037c874
commit ac3dff8818
13 changed files with 129 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Common/operator_base.h
View File

@ -114,6 +114,8 @@ public:
//! Get the cell center coordinate usable for material averaging (Warning, may not be the yee cell center) //! 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 bool GetCellCenterMaterialAvgCoord(const int pos[3], double coord[3]) const = 0;
virtual vector<CSPrimitives*> GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type=CSProperties::ANY) const = 0;
//! Set the background material (default is vacuum) //! 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); virtual void SetBackgroundMaterial(double epsR=0, double mueR=0, double kappa=0, double sigma=0, double density=0);

4
Common/processfields_sar.cpp
View File

@ -203,6 +203,7 @@ void ProcessFieldsSAR::DumpFDData()
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
orig_pos[1] = posLines[1][pos[1]]; orig_pos[1] = posLines[1][pos[1]];
vector<CSPrimitives*> vPrims = Op->GetPrimitivesBoundBox(orig_pos[0], orig_pos[1], -1, CSProperties::MATERIAL);
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
orig_pos[2] = posLines[2][pos[2]]; orig_pos[2] = posLines[2][pos[2]];
@ -211,7 +212,8 @@ void ProcessFieldsSAR::DumpFDData()
cell_density[pos[0]][pos[1]][pos[2]] = 0.0; cell_density[pos[0]][pos[1]][pos[2]] = 0.0;
Op->GetCellCenterMaterialAvgCoord(orig_pos, coord); Op->GetCellCenterMaterialAvgCoord(orig_pos, coord);
prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); prop = CSX->GetPropertyByCoordPriority(coord, vPrims);
// prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
if (prop!=0) if (prop!=0)
{ {
matProp = dynamic_cast<CSPropMaterial*>(prop); matProp = dynamic_cast<CSPropMaterial*>(prop);

4
FDTD/extensions/operator_ext_conductingsheet.cpp
View File

@ -61,6 +61,7 @@ bool Operator_Ext_ConductingSheet::BuildExtension()
{ {
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
b_pos_on = false; b_pos_on = false;
@ -86,7 +87,8 @@ bool Operator_Ext_ConductingSheet::BuildExtension()
if (m_CC_R0_included && (n==2) && (pos[0]==0)) if (m_CC_R0_included && (n==2) && (pos[0]==0))
disable_pos = true; disable_pos = true;
CSProperties* prop = m_Op->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<CSPropConductingSheet*>(prop); CSPropConductingSheet* cs_prop = dynamic_cast<CSPropConductingSheet*>(prop);
if (cs_prop) if (cs_prop)
{ {

6
FDTD/extensions/operator_ext_cylinder.cpp
View File

@ -58,13 +58,15 @@ bool Operator_Ext_Cylinder::BuildExtension()
double inEC[4]; double inEC[4];
double dT = m_Op->GetTimestep(); double dT = m_Op->GetTimestep();
pos[0]=0; pos[0]=0;
vector<CSPrimitives*> vPrims_metal = m_Op->GetPrimitivesBoundBox(pos[0], -1, -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
for (pos[2]=0; pos[2]<m_Op->GetNumberOfLines(2,true); ++pos[2]) for (pos[2]=0; pos[2]<m_Op->GetNumberOfLines(2,true); ++pos[2])
{ {
double C=0; double C=0;
double G=0; double G=0;
vector<CSPrimitives*> vPrims_mat = m_Op->GetPrimitivesBoundBox(pos[0], -1, pos[2], CSProperties::MATERIAL);
for (pos[1]=0; pos[1]<m_Op->GetNumberOfLines(1,true)-2; ++pos[1]) for (pos[1]=0; pos[1]<m_Op->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]; C+=inEC[0];
G+=inEC[1]; G+=inEC[1];
} }
@ -80,7 +82,7 @@ bool Operator_Ext_Cylinder::BuildExtension()
//search for metal on z-axis //search for metal on z-axis
m_Op_Cyl->GetYeeCoords(2,pos,coord,false); 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)
{ {
if (prop->GetType()==CSProperties::METAL) //set to PEC if (prop->GetType()==CSProperties::METAL) //set to PEC

7
FDTD/extensions/operator_ext_lorentzmaterial.cpp
View File

@ -204,6 +204,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
{ {
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
unsigned int index = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]); unsigned int index = m_Op->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) if (m_Op->GetVI(n,pos[0],pos[1],pos[2])==0)
continue; 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())) if ((mat = prop->ToLorentzMaterial()))
{ {
w_plasma = mat->GetEpsPlasmaFreqWeighted(order,n,coord) * 2 * PI; 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) if (m_Op->GetIV(n,pos[0],pos[1],pos[2])==0)
continue; 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())) if ((mat = prop->ToLorentzMaterial()))
{ {
w_plasma = mat->GetMuePlasmaFreqWeighted(order,n,coord) * 2 * PI; w_plasma = mat->GetMuePlasmaFreqWeighted(order,n,coord) * 2 * PI;

9
FDTD/extensions/operator_ext_mur_abc.cpp
View File

@ -140,13 +140,20 @@ bool Operator_Ext_Mur_ABC::BuildExtension()
else else
coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) - delta/2 / m_Op->GetGridDelta(); 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]<m_numLines[0]; ++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_numLines[0]; ++pos[m_nyP])
{ {
posBB[m_nyP]=pos[m_nyP];
vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(posBB[0], posBB[1], posBB[2], CSProperties::MATERIAL);
coord[m_nyP] = m_Op->GetDiscLine(m_nyP,pos[m_nyP]); coord[m_nyP] = m_Op->GetDiscLine(m_nyP,pos[m_nyP]);
for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
coord[m_nyPP] = m_Op->GetDiscLine(m_nyPP,pos[m_nyPP]); coord[m_nyPP] = m_Op->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) if (prop)
{ {
CSPropMaterial* mat = prop->ToMaterial(); CSPropMaterial* mat = prop->ToMaterial();

3
FDTD/extensions/operator_ext_upml.cpp
View File

@ -391,12 +391,13 @@ bool Operator_Ext_UPML::BuildExtension()
for (loc_pos[1]=0; loc_pos[1]<m_numLines[1]; ++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_StartPos[1]; pos[1] = loc_pos[1] + m_StartPos[1];
vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
for (loc_pos[2]=0; loc_pos[2]<m_numLines[2]; ++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_StartPos[2]; pos[2] = loc_pos[2] + m_StartPos[2];
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
{ {
m_Op->Calc_EffMatPos(n,pos,eff_Mat); m_Op->Calc_EffMatPos(n,pos,eff_Mat,vPrims);
CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i); CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i);
nP = (n+1)%3; nP = (n+1)%3;
nPP = (n+2)%3; nPP = (n+2)%3;

131
FDTD/operator.cpp
View File

@ -753,12 +753,13 @@ void Operator::DumpMaterial2File(string filename)
{ {
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
{ {
double inMat[4]; double inMat[4];
Calc_EffMatPos(n, pos, inMat); Calc_EffMatPos(n, pos, inMat, vPrims);
epsilon[n][pos[0]][pos[1]][pos[2]] = inMat[0]/__EPS0__; epsilon[n][pos[0]][pos[1]][pos[2]] = inMat[0]/__EPS0__;
mue[n][pos[0]][pos[1]][pos[2]] = inMat[2]/__MUE0__; mue[n][pos[0]][pos[1]][pos[2]] = inMat[2]/__MUE0__;
kappa[n][pos[0]][pos[1]][pos[2]] = inMat[1]; kappa[n][pos[0]][pos[1]][pos[2]] = inMat[1];
@ -1165,10 +1166,10 @@ void Operator::ApplyMagneticBC(bool* dirs)
} }
} }
bool Operator::Calc_ECPos(int ny, const unsigned int* pos, double* EC) const bool Operator::Calc_ECPos(int ny, const unsigned int* pos, double* EC, vector<CSPrimitives*> vPrims) const
{ {
double EffMat[4]; double EffMat[4];
Calc_EffMatPos(ny,pos,EffMat); Calc_EffMatPos(ny,pos,EffMat, vPrims);
if (m_epsR) if (m_epsR)
m_epsR[ny][pos[0]][pos[1]][pos[2]] = EffMat[0]; 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; 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<CSPrimitives*> 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<CSPropMaterial*>(prop); CSPropMaterial* mat = dynamic_cast<CSPropMaterial*>(prop);
if (mat) 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<CSPrimitives *> vPrims) const
{ {
int n=ny; int n=ny;
double coord[3]; double coord[3];
@ -1302,8 +1310,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff
if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{ {
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0)*A_n; EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1)*A_n; EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=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)) if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{ {
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0)*A_n; EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1)*A_n; EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=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)) if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{ {
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0)*A_n; EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1)*A_n; EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=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)) if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{ {
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, coord, 0)*A_n; EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, coord, 1)*A_n; EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
area+=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)) if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{ {
delta_ny = GetNodeWidth(n,loc_pos,true); delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] += delta_ny / GetMaterial(n, coord, 2); EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims);
sigma = GetMaterial(n, coord, 3); sigma = GetMaterial(n, coord, 3, vPrims);
if (sigma) if (sigma)
EffMat[3] += delta_ny / sigma; EffMat[3] += delta_ny / sigma;
else else
@ -1366,8 +1374,8 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff
if (GetCellCenterMaterialAvgCoord(loc_pos,coord)) if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
{ {
delta_ny = GetNodeWidth(n,loc_pos,true); delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] += delta_ny / GetMaterial(n, coord, 2); EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims);
sigma = GetMaterial(n, coord, 3); sigma = GetMaterial(n, coord, 3, vPrims);
if (sigma) if (sigma)
EffMat[3] += delta_ny / sigma; EffMat[3] += delta_ny / sigma;
else else
@ -1388,7 +1396,7 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff
return true; 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<CSPrimitives*> vPrims) const
{ {
int n=ny; int n=ny;
double coord[3]; 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[nP] = coord[nP]+deltaP*0.25;
shiftCoord[nPP] = coord[nPP]+deltaPP*0.25; shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] = GetMaterial(n, shiftCoord, 0)*A_n; EffMat[0] = GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] = GetMaterial(n, shiftCoord, 1)*A_n; EffMat[1] = GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n; area+=A_n;
//shift up-left //shift up-left
@ -1426,8 +1434,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef
--loc_pos[nP]; --loc_pos[nP];
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n; area+=A_n;
//shift down-right //shift down-right
@ -1437,8 +1445,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef
++loc_pos[nP]; ++loc_pos[nP];
--loc_pos[nPP]; --loc_pos[nPP];
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n; area+=A_n;
//shift down-left //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; shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
--loc_pos[nP]; --loc_pos[nP];
A_n = GetNodeArea(ny,loc_pos,true); A_n = GetNodeArea(ny,loc_pos,true);
EffMat[0] += GetMaterial(n, shiftCoord, 0)*A_n; EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
EffMat[1] += GetMaterial(n, shiftCoord, 1)*A_n; EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
area+=A_n; area+=A_n;
EffMat[0]*=__EPS0__/area; 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; shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
--loc_pos[n]; --loc_pos[n];
double delta_ny = GetNodeWidth(n,loc_pos,true); double delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2); EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2, vPrims);
double sigma = GetMaterial(n, shiftCoord, 3); double sigma = GetMaterial(n, shiftCoord, 3, vPrims);
if (sigma) if (sigma)
EffMat[3] = delta_ny / sigma; EffMat[3] = delta_ny / sigma;
else else
@ -1480,8 +1488,8 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef
shiftCoord[nPP] = coord[nPP]+deltaPP*0.5; shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
++loc_pos[n]; ++loc_pos[n];
delta_ny = GetNodeWidth(n,loc_pos,true); delta_ny = GetNodeWidth(n,loc_pos,true);
EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2); EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2, vPrims);
sigma = GetMaterial(n, shiftCoord, 3); sigma = GetMaterial(n, shiftCoord, 3, vPrims);
if (sigma) if (sigma)
EffMat[3] += delta_ny / sigma; EffMat[3] += delta_ny / sigma;
else else
@ -1502,14 +1510,14 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef
return true; 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<CSPrimitives *> vPrims) const
{ {
switch (m_MatAverageMethod) switch (m_MatAverageMethod)
{ {
case QuarterCell: case QuarterCell:
return AverageMatQuarterCell(ny, pos, EffMat); return AverageMatQuarterCell(ny, pos, EffMat, vPrims);
case CentralCell: case CentralCell:
return AverageMatCellCenter(ny, pos, EffMat); return AverageMatCellCenter(ny, pos, EffMat, vPrims);
default: default:
cerr << "Operator:: " << __func__ << ": Error, unknown material averaging method... exit" << endl; cerr << "Operator:: " << __func__ << ": Error, unknown material averaging method... exit" << endl;
exit(1); exit(1);
@ -1717,19 +1725,51 @@ bool Operator::Calc_EC()
cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl; cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl;
return false; return false;
} }
MainOp->SetPos(0,0,0);
Calc_EC_Range(0,numLines[0]-1);
return true;
}
vector<CSPrimitives*> 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<CSPrimitives*> vPrim = this->CSX->GetPrimitivesByBoundBox(boundBox, true, type);
return vPrim;
}
void Operator::Calc_EC_Range(unsigned int xStart, unsigned int xStop)
{
// vector<CSPrimitives*> vPrims = this->CSX->GetAllPrimitives(true, CSProperties::MATERIAL);
unsigned int ipos; unsigned int ipos;
unsigned int pos[3]; unsigned int pos[3];
double inEC[4]; 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]<numLines[2]; ++pos[2]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0]) ipos = MainOp->GetPos(pos[0],pos[1],pos[2]);
for (int n=0; n<3; ++n)
{ {
Calc_ECPos(n,pos,inEC); Calc_ECPos(n,pos,inEC,vPrims);
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
EC_C[n][ipos]=inEC[0]; EC_C[n][ipos]=inEC[0];
EC_G[n][ipos]=inEC[1]; EC_G[n][ipos]=inEC[1];
EC_L[n][ipos]=inEC[2]; EC_L[n][ipos]=inEC[2];
@ -1738,8 +1778,6 @@ bool Operator::Calc_EC()
} }
} }
} }
return true;
} }
void Operator::SetTimestepFactor(double factor) 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]<numLines[1]; ++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (int n=0; n<3; ++n) for (int n=0; n<3; ++n)
{ {
GetYeeCoords(n,pos,coord,false); GetYeeCoords(n,pos,coord,false);
CSProperties* prop = CSX->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)
{ {
if (prop->GetType()==CSProperties::METAL) //set to PEC if (prop->GetType()==CSProperties::METAL) //set to PEC

13
FDTD/operator.h
View File

@ -179,6 +179,8 @@ public:
virtual double CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const; virtual double CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const;
virtual vector<CSPrimitives*> GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type=CSProperties::ANY) const;
protected: protected:
//! use New() for creating a new Operator //! use New() for creating a new Operator
Operator(); Operator();
@ -215,7 +217,7 @@ protected:
double CalcTimestep_Var3(); double CalcTimestep_Var3();
//! Calculate the FDTD equivalent circuit parameter for the given position and direction ny. \sa Calc_EffMat_Pos //! 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<CSPrimitives *> vPrims) const;
//! Get the FDTD raw disc delta, needed by Calc_EffMatPos() \sa Calc_EffMatPos //! 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; virtual double GetRawDiscDelta(int ny, const int pos) const;
//! Get the material at a given coordinate, direction and type from CSX (internal use only) //! 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<CSPrimitives*> vPrims, bool markAsUsed=true) const;
MatAverageMethods m_MatAverageMethod; MatAverageMethods m_MatAverageMethod;
//! Calculate the effective/averaged material properties at the given position and direction ny. //! 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<CSPrimitives*> vPrims) const;
virtual bool AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat) const; virtual bool AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives*> vPrims) const;
virtual bool AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat) const; virtual bool AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives*> vPrims) const;
//! Calc operator at certain \a pos //! Calc operator at certain \a pos
virtual void Calc_ECOperatorPos(int n, unsigned int* pos); virtual void Calc_ECOperatorPos(int n, unsigned int* pos);
@ -254,6 +256,7 @@ protected:
//EC elements, internal only! //EC elements, internal only!
virtual void Init_EC(); virtual void Init_EC();
virtual bool Calc_EC(); virtual bool Calc_EC();
virtual void Calc_EC_Range(unsigned int xStart, unsigned int xStop);
FDTD_FLOAT* EC_C[3]; FDTD_FLOAT* EC_C[3];
FDTD_FLOAT* EC_G[3]; FDTD_FLOAT* EC_G[3];
FDTD_FLOAT* EC_L[3]; FDTD_FLOAT* EC_L[3];

4
FDTD/operator_cylinder.cpp
View File

@ -67,14 +67,14 @@ double Operator_Cylinder::GetRawDiscDelta(int ny, const int pos) const
return Operator_Multithread::GetRawDiscDelta(ny,pos); 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<CSPrimitives*> vPrims, bool markAsUsed) const
{ {
double l_coords[] = {coords[0],coords[1],coords[2]}; double l_coords[] = {coords[0],coords[1],coords[2]};
if (CC_closedAlpha && (coords[1]>GetDiscLine(1,0,false)+2*PI)) if (CC_closedAlpha && (coords[1]>GetDiscLine(1,0,false)+2*PI))
l_coords[1]-=2*PI; l_coords[1]-=2*PI;
if (CC_closedAlpha && (coords[1]<GetDiscLine(1,0,false))) if (CC_closedAlpha && (coords[1]<GetDiscLine(1,0,false)))
l_coords[1] += 2*PI; l_coords[1] += 2*PI;
return Operator_Multithread::GetMaterial(ny,l_coords,MatType,markAsUsed); return Operator_Multithread::GetMaterial(ny,l_coords,MatType,vPrims,markAsUsed);
} }
int Operator_Cylinder::CalcECOperator( DebugFlags debugFlags ) int Operator_Cylinder::CalcECOperator( DebugFlags debugFlags )

2
FDTD/operator_cylinder.h
View File

@ -102,7 +102,7 @@ protected:
virtual double GetRawDiscDelta(int ny, const int pos) const; virtual double GetRawDiscDelta(int ny, const int pos) const;
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<CSPrimitives*> vPrims, bool markAsUsed=true) const;
virtual int CalcECOperator( DebugFlags debugFlags = None ); virtual int CalcECOperator( DebugFlags debugFlags = None );

5
FDTD/operator_cylindermultigrid.cpp
View File

@ -212,18 +212,17 @@ void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, ve
void Operator_CylinderMultiGrid::FillMissingDataStorage() void Operator_CylinderMultiGrid::FillMissingDataStorage()
{ {
unsigned int pos[3]; unsigned int pos[3];
double EffMat[4]; double EffMat[4];
for (int ny=0; ny<3; ++ny) for (int ny=0; ny<3; ++ny)
{ {
for (pos[0]=0; pos[0]<m_Split_Pos-1; ++pos[0]) for (pos[0]=0; pos[0]<m_Split_Pos-1; ++pos[0])
{ {
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
Calc_EffMatPos(ny,pos,EffMat); Calc_EffMatPos(ny,pos,EffMat,vPrims);
if (m_epsR) if (m_epsR)
m_epsR[ny][pos[0]][pos[1]][pos[2]] = EffMat[0]; m_epsR[ny][pos[0]][pos[1]][pos[2]] = EffMat[0];

22
FDTD/operator_multithread.cpp
View File

@ -189,27 +189,7 @@ void Operator_Thread::operator()()
{ {
//************** calculate EC (Calc_EC) ***********************// //************** calculate EC (Calc_EC) ***********************//
m_OpPtr->m_CalcEC_Start->wait(); m_OpPtr->m_CalcEC_Start->wait();
unsigned int ipos; m_OpPtr->Calc_EC_Range(m_start,m_stop);
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]<m_OpPtr->numLines[1]; ++pos[1])
{
for (pos[2]=0; pos[2]<m_OpPtr->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->m_CalcEC_Stop->wait(); m_OpPtr->m_CalcEC_Stop->wait();
//************** calculate EC (Calc_EC) ***********************// //************** calculate EC (Calc_EC) ***********************//