diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 72c5cb5..219ef49 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -822,6 +822,8 @@ bool Operator::Calc_EC() return true; } +#if 0 //use the old timestep-calc (1) or the new one (0) +////Berechnung nach Andreas Rennings Dissertation 2008, Seite 66, Formel 4.52 double Operator::CalcTimestep() { dT=1e200; @@ -862,6 +864,87 @@ double Operator::CalcTimestep() return 0; } +#else + +double min(double* val, unsigned int count) +{ + if (count==0) + return 0.0; + double min = val[0]; + for (unsigned int n=1;nSetReflection2Cell(); + for (int n=0;n<3;++n) + { + int nP = (n+1)%3; + int nPP = (n+2)%3; + + for (pos[2]=1;pos[2]ResetShift(); + ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); + wqp = 1/(EC_L[nPP][ipos]*EC_C[n][MainOp->GetShiftedPos(nP ,1)]) + 1/(EC_L[nPP][ipos]*EC_C[n][ipos]); + wqp += 1/(EC_L[nP ][ipos]*EC_C[n][MainOp->GetShiftedPos(nPP,1)]) + 1/(EC_L[nP ][ipos]*EC_C[n][ipos]); + ipos = MainOp->Shift(nP,-1); + wqp += 1/(EC_L[nPP][ipos]*EC_C[n][MainOp->GetShiftedPos(nP ,1)]) + 1/(EC_L[nPP][ipos]*EC_C[n][ipos]); + ipos = MainOp->Shift(nPP,-1); + wqp += 1/(EC_L[nP ][ipos]*EC_C[n][MainOp->GetShiftedPos(nPP,1)]) + 1/(EC_L[nP ][ipos]*EC_C[n][ipos]); + + MainOp->ResetShift(); + ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); + wt_4[0] = 1/(EC_L[nPP][ipos] *EC_C[nP ][ipos]); + wt_4[1] = 1/(EC_L[nPP][MainOp->GetShiftedPos(nP ,-1)] *EC_C[nP ][ipos]); + wt_4[2] = 1/(EC_L[nP ][ipos] *EC_C[nPP][ipos]); + wt_4[3] = 1/(EC_L[nP ][MainOp->GetShiftedPos(nPP,-1)] *EC_C[nPP][ipos]); + + wt1 = wt_4[0]+wt_4[1]+wt_4[2]+wt_4[3] - 2*min(wt_4,4); + + MainOp->ResetShift(); + ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); + wt_4[0] = 1/(EC_L[nPP][ipos] *EC_C[nP ][MainOp->GetShiftedPos(n,1)]); + wt_4[1] = 1/(EC_L[nPP][MainOp->GetShiftedPos(nP ,-1)] *EC_C[nP ][MainOp->GetShiftedPos(n,1)]); + wt_4[2] = 1/(EC_L[nP ][ipos] *EC_C[nPP][MainOp->GetShiftedPos(n,1)]); + wt_4[3] = 1/(EC_L[nP ][MainOp->GetShiftedPos(nPP,-1)] *EC_C[nPP][MainOp->GetShiftedPos(n,1)]); + + wt2 = wt_4[0]+wt_4[1]+wt_4[2]+wt_4[3] - 2*min(wt_4,4); + + w_total = wqp + wt1 + wt2; + newT = 2/sqrt( w_total ); + if ((newT0.0)) + dT=newT; + } + } + } + } + if (dT==0) + { + cerr << "Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit!" << endl; + exit(3); + } +// cerr << "Operator Timestep: " << dT << endl; + return 0; +} + +#endif + bool Operator::CalcFieldExcitation() { if (dT==0) diff --git a/tools/AdrOp.cpp b/tools/AdrOp.cpp index b8201fc..6c2c5fd 100644 --- a/tools/AdrOp.cpp +++ b/tools/AdrOp.cpp @@ -347,9 +347,13 @@ bool AdrOp::CheckShift(int ny, int step) else return false; } -unsigned int AdrOp::GetShiftedPos() +unsigned int AdrOp::GetShiftedPos(int ny, int step) { - return GetPos(iIshift,iJshift,iKshift); + if ((ny<0) || (ny>2)) + return GetPos(iIshift,iJshift,iKshift); + int pos[3] = {iIshift, iJshift, iKshift}; + pos[ny]+=step; + return GetPos(pos[0],pos[1],pos[2]); } void AdrOp::ResetShift() diff --git a/tools/AdrOp.h b/tools/AdrOp.h index 6568175..269866c 100644 --- a/tools/AdrOp.h +++ b/tools/AdrOp.h @@ -73,8 +73,8 @@ public: ///Set a checked shift in ny direction (e.g. 0 for i-direction) /*!Shift set by this methode will be ignored by methode GetPos*/ bool CheckShift(int ny, int step); - ///Returns the current 1-dim position including shift by methode "Shift" - unsigned int GetShiftedPos(); + ///Returns the current 1-dim position including shift by methode "Shift" + additional (transitory) shift + unsigned int GetShiftedPos(int ny=-1, int step=0); ///Reset shift set by "Shift"-methode void ResetShift(); ///Iterates through AdrOp; --- obsolete ---