diff --git a/Common/operator_base.h b/Common/operator_base.h index 523cb14..f8a4825 100644 --- a/Common/operator_base.h +++ b/Common/operator_base.h @@ -44,9 +44,6 @@ public: //! Get the grid drawing unit in m virtual double GetGridDelta() const =0; - //! Get the mesh delta times the grid delta for a 3D position (unit is meter) - virtual double GetMeshDelta(int n, const unsigned int* pos, bool dualMesh=false) const =0; - //! Get the disc line in \a n direction (in drawing units) virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const =0; diff --git a/Common/process_efield.cpp b/Common/process_efield.cpp index c55a10f..4e7713e 100644 --- a/Common/process_efield.cpp +++ b/Common/process_efield.cpp @@ -91,7 +91,7 @@ int ProcessEField::Process() file << setprecision(m_precision) << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep(); for (int n=0; n<3; n++) { - FDTD_FLOAT field = Eng->GetVolt(n,start) / Op->GetMeshDelta(n,start); + FDTD_FLOAT field = Eng->GetVolt(n,start) / Op->GetEdgeLength(n,start); field *= m_weight; // TD_Values.push_back(voltage); file << "\t" << field; @@ -108,7 +108,7 @@ int ProcessEField::Process() double T = (double)Eng->GetNumberOfTimesteps() * Op->GetTimestep(); for (int pol=0; pol<3; pol++) { - FDTD_FLOAT field = Eng->GetVolt(pol,start) / Op->GetMeshDelta(pol,start); + FDTD_FLOAT field = Eng->GetVolt(pol,start) / Op->GetEdgeLength(pol,start); field *= m_weight; for (size_t n=0; nGetNumberOfTimesteps()+0.5)*Op->GetTimestep(); for (int n=0; n<3; n++) { - FDTD_FLOAT field = Eng->GetCurr(n,start) / Op->GetMeshDelta(n,start,true); + FDTD_FLOAT field = Eng->GetCurr(n,start) / Op->GetEdgeLength(n,start,true); field *= m_weight; // TD_Values.push_back(voltage); file << "\t" << field; @@ -92,7 +92,7 @@ int ProcessHField::Process() double T = (double)Eng->GetNumberOfTimesteps() * Op->GetTimestep(); for (int pol=0; pol<3; pol++) { - FDTD_FLOAT field = Eng->GetCurr(pol,start) / Op->GetMeshDelta(pol,start,true); + FDTD_FLOAT field = Eng->GetCurr(pol,start) / Op->GetEdgeLength(pol,start,true); field *= m_weight; for (size_t n=0; n enlarge for (int i=0; i<3; i++) { - double delta = min( Op->GetMeshDelta( i, start,dualMesh ), Op->GetMeshDelta( i, stop,dualMesh ) ) / Op->GetGridDelta() / 4.0; + double delta = min( Op->GetEdgeLength( i, start,dualMesh ), Op->GetEdgeLength( i, stop,dualMesh ) ) / Op->GetGridDelta() / 4.0; s1[i] -= delta; s2[i] += delta; } diff --git a/FDTD/engine_interface_fdtd.cpp b/FDTD/engine_interface_fdtd.cpp index d787fad..da87055 100644 --- a/FDTD/engine_interface_fdtd.cpp +++ b/FDTD/engine_interface_fdtd.cpp @@ -48,7 +48,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c case NODE_INTERPOLATE: for (int n=0; n<3; ++n) { - delta = m_Op->GetMeshDelta(n,iPos); + delta = m_Op->GetEdgeLength(n,iPos); out[n] = m_Eng->GetVolt(n,iPos); if (delta==0) { @@ -61,7 +61,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c continue; } --iPos[n]; - double deltaDown = m_Op->GetMeshDelta(n,iPos); + double deltaDown = m_Op->GetEdgeLength(n,iPos); double deltaRel = delta / (delta+deltaDown); out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetVolt(n,iPos)/deltaDown*deltaRel; ++iPos[n]; @@ -77,19 +77,19 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c out[n] = 0; //electric field outside the field domain is always zero continue; } - delta = m_Op->GetMeshDelta(n,iPos); + delta = m_Op->GetEdgeLength(n,iPos); if (delta) out[n]=m_Eng->GetVolt(n,iPos)/delta; ++iPos[nP]; - delta = m_Op->GetMeshDelta(n,iPos); + delta = m_Op->GetEdgeLength(n,iPos); if (delta) out[n]+=m_Eng->GetVolt(n,iPos)/delta; ++iPos[nPP]; - delta = m_Op->GetMeshDelta(n,iPos); + delta = m_Op->GetEdgeLength(n,iPos); if (delta) out[n]+=m_Eng->GetVolt(n,iPos)/delta; --iPos[nP]; - delta = m_Op->GetMeshDelta(n,iPos); + delta = m_Op->GetEdgeLength(n,iPos); if (delta) out[n]+=m_Eng->GetVolt(n,iPos)/delta; --iPos[nPP]; @@ -123,13 +123,13 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c out[n] = 0; continue; } - out[n]=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + out[n]=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true); --iPos[nP]; - out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true); --iPos[nPP]; - out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true); ++iPos[nP]; - out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetMeshDelta(n,iPos,true); + out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true); ++iPos[nPP]; out[n]/=4; } @@ -137,7 +137,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c case CELL_INTERPOLATE: for (int n=0; n<3; ++n) { - delta = m_Op->GetMeshDelta(n,iPos,true); + delta = m_Op->GetEdgeLength(n,iPos,true); out[n] = m_Eng->GetCurr(n,iPos); if ((pos[n]>=m_Op->GetNumberOfLines(n)-1)) { @@ -145,7 +145,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c continue; } ++iPos[n]; - double deltaUp = m_Op->GetMeshDelta(n,iPos,true); + double deltaUp = m_Op->GetEdgeLength(n,iPos,true); double deltaRel = delta / (delta+deltaUp); out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetCurr(n,iPos)/deltaUp*deltaRel; --iPos[n]; diff --git a/FDTD/extensions/operator_ext_cylinder.cpp b/FDTD/extensions/operator_ext_cylinder.cpp index 307024c..9cf8108 100644 --- a/FDTD/extensions/operator_ext_cylinder.cpp +++ b/FDTD/extensions/operator_ext_cylinder.cpp @@ -64,8 +64,8 @@ bool Operator_Ext_Cylinder::BuildExtension() for (pos[1]=0; pos[1]GetOriginalNumLines(1)-1; ++pos[1]) { m_Op_Cyl->Calc_ECPos(2,pos,inEC); - C+=inEC[0]*0.5; - G+=inEC[1]*0.5; + C+=inEC[0]; + G+=inEC[1]; } m_Op->SetVV(2,0,0,pos[2], 1); vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C); diff --git a/FDTD/extensions/operator_ext_mur_abc.cpp b/FDTD/extensions/operator_ext_mur_abc.cpp index ee271b5..26a209c 100644 --- a/FDTD/extensions/operator_ext_mur_abc.cpp +++ b/FDTD/extensions/operator_ext_mur_abc.cpp @@ -124,7 +124,7 @@ bool Operator_Ext_Mur_ABC::BuildExtension() double dT = m_Op->GetTimestep(); unsigned int pos[] = {0,0,0}; pos[m_ny] = m_LineNr; - double delta = fabs(m_Op->GetMeshDelta(m_ny,pos)); + double delta = fabs(m_Op->GetEdgeLength(m_ny,pos)); double coord[] = {0,0,0}; coord[0] = m_Op->GetDiscLine(0,pos[0]); coord[1] = m_Op->GetDiscLine(1,pos[1]); diff --git a/FDTD/extensions/operator_ext_pml_sf.cpp b/FDTD/extensions/operator_ext_pml_sf.cpp index 41be114..9841801 100644 --- a/FDTD/extensions/operator_ext_pml_sf.cpp +++ b/FDTD/extensions/operator_ext_pml_sf.cpp @@ -227,7 +227,7 @@ void Operator_Ext_PML_SF_Plane::SetDirection(int ny, bool top_ny) m_LineNr = (unsigned int)((int)m_top * (int)(m_Op->GetNumberOfLines(m_ny)-1)); pos[m_ny] = m_LineNr; - m_pml_delta = m_Op->GetMeshDelta(m_ny,pos); + m_pml_delta = m_Op->GetEdgeLength(m_ny,pos); } void Operator_Ext_PML_SF_Plane::SetPMLLength(int width) diff --git a/FDTD/extensions/operator_ext_upml.cpp b/FDTD/extensions/operator_ext_upml.cpp index b461e49..8672630 100644 --- a/FDTD/extensions/operator_ext_upml.cpp +++ b/FDTD/extensions/operator_ext_upml.cpp @@ -277,17 +277,17 @@ void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm, depth = width - (m_Op->GetDiscLine(n,pos[n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta(); if (n==ny) - depth-=m_Op->GetMeshDelta(n,pos)/2; + depth-=m_Op->GetEdgeLength(n,pos)/2; double vars[5] = {depth, width/m_Size[2*n], width, Zm, m_Size[2*n]}; if (depth>0) kappa_v[n] = m_GradingFunction->Eval(vars); else kappa_v[n]=0; if (n==ny) - depth+=m_Op->GetMeshDelta(n,pos)/2; + depth+=m_Op->GetEdgeLength(n,pos)/2; if (n!=ny) - depth-=m_Op->GetMeshDelta(n,pos)/2; + depth-=m_Op->GetEdgeLength(n,pos)/2; if (depth<0) depth=0; vars[0]=depth; @@ -302,17 +302,17 @@ void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm, depth = width - (m_Op->GetDiscLine(n,m_Op->GetOriginalNumLines(n)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta(); if (n==ny) - depth+=m_Op->GetMeshDelta(n,pos)/2; + depth+=m_Op->GetEdgeLength(n,pos)/2; double vars[5] = {depth, width/(m_Size[2*n]), width, Zm, m_Size[2*n]}; if (depth>0) kappa_v[n] = m_GradingFunction->Eval(vars); else kappa_v[n]=0; if (n==ny) - depth-=m_Op->GetMeshDelta(n,pos)/2; + depth-=m_Op->GetEdgeLength(n,pos)/2; if (n!=ny) - depth+=m_Op->GetMeshDelta(n,pos)/2; + depth+=m_Op->GetEdgeLength(n,pos)/2; if (depth>width) depth=0; vars[0]=depth; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index ab666f9..5e4a000 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -97,24 +97,43 @@ void Operator::Reset() Operator_Base::Reset(); } -double Operator::GetMeshDelta(int n, const unsigned int* pos, bool dualMesh) const -{ - if ((n<0) || (n>2)) return 0.0; - if (pos[n]>=numLines[n]) return 0.0; - if (dualMesh==false) - return fabs(MainOp->GetIndexDelta(n,pos[n]))*gridDelta; - else - return fabs(MainOp->GetIndexWidth(n,pos[n]))*gridDelta; -} - double Operator::GetDiscLine(int n, unsigned int pos, bool dualMesh) const { if ((n<0) || (n>2)) return 0.0; if (pos>=numLines[n]) return 0.0; if (dualMesh==false) return discLines[n][pos]; + + // return dual mesh node + if (pos2)) return 0.0; + if (pos[n]>=numLines[n]) return 0.0; + double delta=0; + if (dualMesh==false) + { + if (pos[n]GetIndexDelta(n,pos))); + { + if (pos[n]>0) + delta = GetDiscLine(n,pos[n],true) - GetDiscLine(n,pos[n]-1,true); + else + delta = GetDiscLine(n,1,false) - GetDiscLine(n,0,false); + return delta*gridDelta; + } } double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const @@ -1242,7 +1261,7 @@ bool Operator::CalcFieldExcitation() { if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))//&& (pos[n]GetWeightedExcitation(n,volt_coord)*GetMeshDelta(n,pos);// delta[n]*gridDelta; + amp = elec->GetWeightedExcitation(n,volt_coord)*GetEdgeLength(n,pos);// delta[n]*gridDelta; if (amp!=0) { volt_vExcit.push_back(amp); @@ -1284,7 +1303,7 @@ bool Operator::CalcFieldExcitation() { if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==2) || (elec->GetExcitType()==3) )) { - amp = elec->GetWeightedExcitation(n,curr_coord)*GetMeshDelta(n,pos,true);// delta[n]*gridDelta; + amp = elec->GetWeightedExcitation(n,curr_coord)*GetEdgeLength(n,pos,true);// delta[n]*gridDelta; if (amp!=0) { curr_vExcit.push_back(amp); diff --git a/FDTD/operator.h b/FDTD/operator.h index 06edb41..dc7bd85 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -93,20 +93,18 @@ public: virtual void ShowExtStat() const; virtual double GetGridDelta() const {return gridDelta;} - //! Get the mesh delta times the grid delta for a 3D position (unit is meter) - virtual double GetMeshDelta(int n, const unsigned int* pos, bool dualMesh=false) const; //! Get the disc line in \a n direction (in drawing units) virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const; //! 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 GetMeshDelta(ny,pos,!dualMesh);} + virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetEdgeLength(ny,pos,!dualMesh);} //! 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 length of an FDTD edge (unit is meter). - virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetMeshDelta(ny,pos,dualMesh);} + virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const; //! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos /*! diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index 9e0d122..cf76d31 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -85,68 +85,48 @@ string Operator_Cylinder::GetDirName(int ny) const return ""; } -double Operator_Cylinder::GetMeshDelta(int n, const unsigned int* pos, bool dualMesh) const -{ - double delta = Operator_Multithread::GetMeshDelta(n,pos,dualMesh); - if (delta==0) return delta; - if (n==1) - { - return delta * GetDiscLine(0,pos[0],dualMesh); - } - return delta; -} - double Operator_Cylinder::GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh) const { if ((ny<0) || (ny>2)) return 0.0; if (pos[ny]>=numLines[ny]) return 0.0; - double width = 0; - if (dualMesh) - width = fabs(MainOp->GetIndexDelta(ny,pos[ny]))*gridDelta; - else - width = fabs(MainOp->GetIndexWidth(ny,pos[ny]))*gridDelta; + double width = Operator_Multithread::GetEdgeLength(ny,pos,!dualMesh); if (ny==1) width *= GetDiscLine(0,pos[0],dualMesh); return width; } - double Operator_Cylinder::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const { if (pos[ny]>=numLines[ny]) return 0.0; if (pos[0]>=numLines[0]) return 0.0; if (ny==2) { - double da = Operator_Multithread::GetMeshDelta(1,pos,dualMesh)/gridDelta; + double da = Operator_Multithread::GetEdgeLength(1,pos,dualMesh)/gridDelta; double r1,r2; - if (!dualMesh) + + if (dualMesh) { - r1 = (discLines[0][pos[0]] - fabs(MainOp->GetIndexDelta(0,pos[0]-1))/2.0)*gridDelta; - r2 = (discLines[0][pos[0]] + fabs(MainOp->GetIndexDelta(0,pos[0] ))/2.0)*gridDelta; + r1 = GetDiscLine(0,pos[0],false)*gridDelta; + r2 = r1 + GetEdgeLength(0,pos,false); } else { - r1 = discLines[0][pos[0]]*gridDelta; - r2 = (discLines[0][pos[0]] + fabs(MainOp->GetIndexDelta(0,pos[0])))*gridDelta; + r2 = GetDiscLine(0,pos[0],!dualMesh)*gridDelta; + r1 = r2 - GetEdgeLength(0,pos,true); } - if (r1<0) - return da * pow(r2,2); - return da/2* (pow(r2,2) - pow(r1,2)); + + if (r1<=0) + return da/2 * pow(r2,2); + else + return da/2* (pow(r2,2) - pow(r1,2)); } - if (ny==0) - { - if (dualMesh) - return fabs(MainOp->GetIndexDelta(1,pos[1]) * MainOp->GetIndexDelta(2,pos[2]) * GetDiscLine(0,pos[0],true) * gridDelta * gridDelta); - else - return fabs(MainOp->GetIndexDelta(1,pos[1]) * MainOp->GetIndexDelta(2,pos[2]) * GetDiscLine(0,pos[0],false) * gridDelta * gridDelta); - } return Operator_Multithread::GetNodeArea(ny,pos,dualMesh); } double Operator_Cylinder::GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh) const { - double length = Operator_Multithread::GetMeshDelta(ny,pos,dualMesh); + double length = Operator_Multithread::GetEdgeLength(ny,pos,dualMesh); if (ny!=1) return length; return length * GetDiscLine(0,pos[0],dualMesh); diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h index d14a590..311ed9d 100644 --- a/FDTD/operator_cylinder.h +++ b/FDTD/operator_cylinder.h @@ -44,16 +44,13 @@ public: //! Get the name for the given direction: 0 -> rho, 1 -> alpha, 2 -> z virtual string GetDirName(int ny) const; - //! Get the mesh delta times the grid delta for a 3D position, including radius corrected alpha-mesh width - virtual double GetMeshDelta(int n, const unsigned int* pos, bool dualMesh=false) const; - //! Get the node width for a given direction \a n and a given mesh posisition \a pos virtual double GetNodeWidth(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 posisition \a pos virtual double GetNodeArea(int n, const unsigned int* pos, bool dualMesh=false) const; - //! Get the length of an FDTD edge. + //! Get the length of an FDTD edge, including radius corrected alpha-mesh width. virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const; //! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos