From 8316b1c2bd2058bb6bdd9d076fc29b42ed734e75 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 29 Jul 2010 18:30:50 +0200 Subject: [PATCH] Operator: GetNodeArea & Update/fix in Calc_EffMatPos method - new method GetNodeArea will return the area of a node for a given direction - methods Calc_ECPos & Calc_EffMatPos now const - Calc_EffMatPos in Operator class updated to use new functions GetMeshDelta & GetNodeArea - Calc_EffMatPos introduced (fixed) in Operator_Cylinder (need some testing) - treatment of E_alpha=0 at r==0 moved from Calc_ECPos to Electric-BC --- FDTD/operator.cpp | 37 +++++++++------------- FDTD/operator.h | 9 ++++-- FDTD/operator_cylinder.cpp | 63 ++++++++++++++++++++++++++++++-------- FDTD/operator_cylinder.h | 9 ++++-- 4 files changed, 79 insertions(+), 39 deletions(-) diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 8e2d2a1..5217dbe 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -138,12 +138,14 @@ double Operator::GetDiscLine(int n, int pos, bool dualMesh) const if (dualMesh==false) return discLines[n][pos]; else - { - if (pos<(int)numLines[n]-1) - return 0.5*(discLines[n][pos+1]+discLines[n][pos]); - else - return 0.5*(discLines[n][pos]+discLines[n][pos-1]); - } + return (discLines[n][pos] + 0.5*fabs(MainOp->GetIndexDelta(n,pos))); +} + +double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const +{ + int nyP = (ny+1)%3; + int nyPP = (ny+2)%3; + return GetMeshDelta(nyP,pos,!dualMesh) * GetMeshDelta(nyPP,pos,!dualMesh); } bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, bool* inside) @@ -645,7 +647,7 @@ void Operator::ApplyMagneticBC(bool* dirs) } -bool Operator::Calc_ECPos(int n, unsigned int* pos, double* inEC) +bool Operator::Calc_ECPos(int n, const unsigned int* pos, double* inEC) const { double coord[3]; double shiftCoord[3]; @@ -778,26 +780,15 @@ bool Operator::Calc_ECPos(int n, unsigned int* pos, double* inEC) return true; } -bool Operator::Calc_EffMatPos(int n, unsigned int* pos, double* inMat) +bool Operator::Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const { - int nP = (n+1)%3; - int nPP = (n+2)%3; - - double delta=MainOp->GetIndexDelta(n,pos[n]); - double deltaP=MainOp->GetIndexDelta(nP,pos[nP]); - double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]); - - double delta_M=MainOp->GetIndexDelta(n,pos[n]-1); - double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1); - double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1); - this->Calc_ECPos(n,pos,inMat); - inMat[0] *= fabs(delta)/(0.25*(fabs(deltaP_M) + fabs(deltaP))*(fabs(deltaPP_M) + fabs(deltaPP)))/gridDelta; - inMat[1] *= fabs(delta)/(0.25*(fabs(deltaP_M) + fabs(deltaP))*(fabs(deltaPP_M) + fabs(deltaPP)))/gridDelta; + inMat[0] *= GetMeshDelta(n,pos)/GetNodeArea(n,pos); + inMat[1] *= GetMeshDelta(n,pos)/GetNodeArea(n,pos); - inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta; - inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta; + inMat[2] *= GetMeshDelta(n,pos,true)/GetNodeArea(n,pos,true); + inMat[3] *= GetMeshDelta(n,pos,true)/GetNodeArea(n,pos,true); return true; } diff --git a/FDTD/operator.h b/FDTD/operator.h index 463a168..8ee5a0e 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -83,6 +83,11 @@ public: //! 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 area for a given direction \a n and a given mesh posisition \a pos + virtual double GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {return GetNodeArea(ny,(const int*)pos,dualMesh);} + //! Get the node area for a given direction \a n and a given mesh posisition \a pos + virtual double GetNodeArea(int ny, const int pos[3], bool dualMesh = false) const; + virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL); virtual void AddExtension(Operator_Extension* op_ext); @@ -131,8 +136,8 @@ protected: //EC elements, internal only! virtual void Init_EC(); virtual bool Calc_EC(); - virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC); - virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat); + virtual bool Calc_ECPos(int n, const unsigned int* pos, double* inEC) const; + virtual bool Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const; double* EC_C[3]; double* EC_G[3]; double* EC_L[3]; diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index c8b0505..d956b3e 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -88,6 +88,28 @@ double Operator_Cylinder::GetMeshDelta(int n, const int* pos, bool dualMesh) con return delta; } +double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) const +{ + if (ny==2) + { + double da = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(1,pos,dualMesh)/gridDelta; + double r1,r2; + 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; + } + else + { + r1 = discLines[0][pos[0]]*gridDelta; + r2 = (discLines[0][pos[0]] + fabs(MainOp->GetIndexDelta(0,pos[0])))*gridDelta; + } + if (r1<0) + return da * pow(r2,2); + return da/2* (pow(r2,2) - pow(r1,2)); + } + return __OP_CYLINDER_BASE_CLASS__::GetNodeArea(ny,pos,dualMesh); +} bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) { @@ -138,8 +160,16 @@ void Operator_Cylinder::ApplyElectricBC(bool* dirs) } if (CC_R0_included) { - // no special treatment necessary - // operator for z-direction at r=0 will be calculated and set separately + // 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]