diff --git a/FDTD/engine_cylinder.cpp b/FDTD/engine_cylinder.cpp index 8796937..579f079 100644 --- a/FDTD/engine_cylinder.cpp +++ b/FDTD/engine_cylinder.cpp @@ -46,7 +46,28 @@ inline void Engine_Cylinder::CloseAlphaVoltages() volt[1][pos[0]][0][pos[2]] = volt[1][pos[0]][last_A_Line][pos[2]]; volt[2][pos[0]][0][pos[2]] = volt[2][pos[0]][last_A_Line][pos[2]]; } + } +} +void Engine_Cylinder::R0IncludeVoltages() +{ + unsigned int pos[3]; + pos[0] = 0; + for (pos[2]=0;pos[2]vv_R0[pos[2]]; + for (pos[1]=0;pos[1]GetClosedAlpha();++pos[1]) + { + volt[2][0][0][pos[2]] += cyl_Op->vi_R0[pos[2]] * curr[1][0][pos[1]][pos[2]]; + } + } + for (pos[1]=0;pos[1]GetR0Included()) + R0IncludeVoltages(); + ApplyVoltageExcite(); CloseAlphaVoltages(); diff --git a/FDTD/engine_cylinder.h b/FDTD/engine_cylinder.h index 88910ef..df93b30 100644 --- a/FDTD/engine_cylinder.h +++ b/FDTD/engine_cylinder.h @@ -35,6 +35,8 @@ protected: virtual inline void CloseAlphaVoltages(); virtual inline void CloseAlphaCurrents(); + virtual inline void R0IncludeVoltages(); + const Operator_Cylinder* cyl_Op; }; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 559f7f6..a9ecbd7 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -102,6 +102,42 @@ unsigned int Operator::CalcNyquistNum(double fmax) return floor(T0/2/dT); } +double Operator::GetMeshDelta(int n, unsigned int* pos, bool dualMesh) const +{ + if ((n<0) || (n>2)) return 0.0; + int i_pos[] = {pos[0],pos[1],pos[2]}; + return GetMeshDelta(n,i_pos,dualMesh); +} + +double Operator::GetMeshDelta(int n, int* pos, bool dualMesh) const +{ + if ((n<0) || (n>2)) 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 +{ + return GetDiscLine(n,(int)pos,dualMesh); +} + +double Operator::GetDiscLine(int n, int pos, bool dualMesh) const +{ + if ((n<0) || (n>2)) return 0.0; + if ((pos<0) || (pos>=(int)numLines[n])) return 0.0; + 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]); + } +} + bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, bool* inside) { bool ok=true; @@ -451,6 +487,7 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo) for (int n=0;n<3;++n) { discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true); + if (n==1) if (numLines[n]<3) {cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl; Reset(); return false;} } MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]); @@ -821,11 +858,16 @@ double Operator::CalcTimestep() ipos_PPM= MainOp->Shift(nPP,-1); MainOp->ResetShift(); newT = 2/sqrt( ( 4/EC_L[nP][ipos] + 4/EC_L[nP][ipos_PPM] + 4/EC_L[nPP][ipos] + 4/EC_L[nPP][ipos_PM]) / EC_C[n][ipos] ); - 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; } @@ -852,13 +894,11 @@ bool Operator::CalcEFieldExcitation() for (pos[0]=0;pos[0]GetIndexDelta(0,pos[0])); - coord[0] = discLines[0][pos[0]]; - coord[1] = discLines[1][pos[1]]; - coord[2] = discLines[2][pos[2]]; -// CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::ELECTRODE | CSProperties::METAL)); - CSProperties* prop = NULL; for (int n=0;n<3;++n) { + coord[0] = discLines[0][pos[0]]; + coord[1] = discLines[1][pos[1]]; + coord[2] = discLines[2][pos[2]]; coord[n]+=delta[n]*0.5; CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::ELECTRODE)); if (prop) @@ -866,9 +906,9 @@ bool Operator::CalcEFieldExcitation() CSPropElectrode* elec = prop->ToElectrode(); if (elec!=NULL) { - if ((elec->GetActiveDir(n)) && (pos[n]GetActiveDir(n)) )//&& (pos[n]GetWeightedExcitation(n,coord)*delta[n]*gridDelta; + amp = elec->GetWeightedExcitation(n,coord)*GetMeshDelta(n,pos);// delta[n]*gridDelta; if (amp!=0) { vExcit.push_back(amp); @@ -886,7 +926,6 @@ bool Operator::CalcEFieldExcitation() } } } - coord[n]-=delta[n]*0.5; } } } diff --git a/FDTD/operator.h b/FDTD/operator.h index 16151cd..c4fd927 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -59,7 +59,7 @@ public: virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];} void SetNyquistNum(unsigned int nyquist) {m_nyquistTS=nyquist;} - unsigned int GetNyquistNum() const {return m_nyquistTS;}; + unsigned int GetNyquistNum() const {return m_nyquistTS;} unsigned int CalcNyquistNum(double fmax); void ShowStat() const; @@ -67,15 +67,22 @@ public: void DumpOperator2File(string filename); void DumpMaterial2File(string filename); - virtual void Reset(); + virtual double GetGridDelta() const {return gridDelta;} + //! Get the mesh delta times the grid delta for a 3D position + virtual double GetMeshDelta(int n, int* pos, bool dualMesh=false) const; + virtual double GetMeshDelta(int n, unsigned int* pos, bool dualMesh=false) const; - bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL); + //! Get the disc line in n direction + virtual double GetDiscLine(int n, int pos, bool dualMesh=false) const; + virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const; + virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL); protected: //! use New() for creating a new Operator Operator(); virtual void Init(); + virtual void Reset(); virtual void InitOperator(); struct Grid_Path @@ -111,14 +118,13 @@ protected: double* EC_R[3]; unsigned int numLines[3]; - - // engine/post-proc needs access -public: double* discLines[3]; double gridDelta; AdrOp* MainOp; AdrOp* DualOp; + // engine/post-proc needs access +public: //EC operator FDTD_FLOAT**** vv; //calc new voltage from old voltage FDTD_FLOAT**** vi; //calc new voltage from old current diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index c643162..1b583cf 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -37,12 +37,26 @@ void Operator_Cylinder::Init() { CC_closedAlpha = false; CC_R0_included = false; + vv_R0 = NULL; + vi_R0 = NULL; Operator::Init(); } void Operator_Cylinder::Reset() { Operator::Reset(); + delete[] vv_R0;vv_R0=NULL; + delete[] vi_R0;vi_R0=NULL; +} + +void Operator_Cylinder::InitOperator() +{ + if (CC_R0_included) + { + vv_R0 = new FDTD_FLOAT[numLines[2]]; + vi_R0 = new FDTD_FLOAT[numLines[2]]; + } + Operator::InitOperator(); } inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const @@ -54,13 +68,23 @@ inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const return numLines[ny]; } +double Operator_Cylinder::GetMeshDelta(int n, int* pos, bool dualMesh) const +{ + double delta = Operator::GetMeshDelta(n,pos,dualMesh); + if (delta==0) return delta; + if (n==1) + { + return delta * GetDiscLine(n,pos[0],dualMesh); + } + return delta; +} + bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) { if (Operator::SetGeometryCSX(geo)==false) return false; double minmaxA = fabs(discLines[1][numLines[1]-1]-discLines[1][0]); -// cerr << minmaxA -2*PI << " < " << (2*PI)/10/numLines[1] << endl; if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width { cout << "Operator_Cylinder::SetGeometryCSX: Alpha is a full 2*PI => closed Cylinder..." << endl; @@ -87,13 +111,41 @@ bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) else if (discLines[0][0]==0.0) { cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl; - cerr << "Operator_Cylinder::SetGeometryCSX: r=0 included not yet implemented... exit... " << endl; exit(1); CC_R0_included=true; } return true; } +int Operator_Cylinder::CalcECOperator() +{ + int val = Operator::CalcECOperator(); + if (val) + return val; + + if (CC_R0_included==false) + return val; + + unsigned int pos[3]; + double inEC[4]; + pos[0]=0; + for (pos[2]=0;pos[2]SetPos(pos[0],pos[1],pos[2]); @@ -118,6 +170,12 @@ inline void Operator_Cylinder::Calc_ECOperatorPos(int n, unsigned int* pos) ii[n][pos[0]][pos[1]][pos[2]] = 0; iv[n][pos[0]][pos[1]][pos[2]] = 0; } + + if (CC_R0_included && (n==2) && (pos[0]==0)) + { + vv[n][pos[0]][pos[1]][pos[2]] = 1; + vi[n][pos[0]][pos[1]][pos[2]] = 0; + } } void Operator_Cylinder::ApplyElectricBC(bool* dirs) @@ -165,7 +223,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) 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); - double geom_factor,A_n; + double geom_factor=0,A_n=0; //******************************* epsilon,kappa averaging *****************************// //shift up-right @@ -176,13 +234,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) switch (n) { case 0: - geom_factor = fabs(deltaPP*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25; + geom_factor = fabs((deltaPP*deltaP/delta)*(coord[0]+fabs(delta)/2))*0.25; break; case 1: geom_factor = fabs(deltaP*deltaPP/(delta*coord[0]))*0.25; break; case 2: - geom_factor = fabs(deltaPP/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0))*0.25; + geom_factor = fabs((deltaPP/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0)))*0.25; break; } geom_factor*=gridDelta; @@ -206,13 +264,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) switch (n) { case 0: - geom_factor = fabs(deltaPP*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25; + geom_factor = fabs((deltaPP*deltaP_M/delta)*(coord[0]+fabs(delta)/2))*0.25; break; case 1: geom_factor = fabs(deltaP_M*deltaPP/(delta*coord[0]))*0.25; break; case 2: - geom_factor = fabs(deltaPP/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0))*0.25; + geom_factor = fabs((deltaPP/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0)))*0.25; break; } geom_factor*=gridDelta; @@ -236,13 +294,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) switch (n) { case 0: - geom_factor = fabs(deltaPP_M*deltaP/delta)*(coord[0]+fabs(delta)/2)*0.25; + geom_factor = fabs((deltaPP_M*deltaP/delta)*(coord[0]+fabs(delta)/2))*0.25; break; case 1: geom_factor = fabs(deltaP*deltaPP_M/(delta*coord[0]))*0.25; break; case 2: - geom_factor = fabs(deltaPP_M/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0))*0.25; + geom_factor = fabs((deltaPP_M/delta) * (pow(coord[0]+fabs(deltaP)/2.0,2.0) - pow(coord[0],2.0)))*0.25; break; } geom_factor*=gridDelta; @@ -266,13 +324,13 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) switch (n) { case 0: - geom_factor = fabs(deltaPP_M*deltaP_M/delta)*(coord[0]+fabs(delta)/2)*0.25; + geom_factor = fabs((deltaPP_M*deltaP_M/delta)*(coord[0]+fabs(delta)/2))*0.25; break; case 1: geom_factor = fabs(deltaP_M*deltaPP_M/(delta*coord[0]))*0.25; break; case 2: - geom_factor = fabs(deltaPP_M/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0))*0.25; + geom_factor = fabs((deltaPP_M/delta) * (pow(coord[0],2.0) - pow(coord[0]-fabs(deltaP_M)/2.0,2.0)))*0.25; break; } geom_factor*=gridDelta; @@ -288,15 +346,12 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) inEC[1] += 0; } - if (CC_R0_included && (n>0) && (coord[0]==0)) + if (CC_R0_included && (n==1) && (coord[0]==0)) { inEC[0]=0; inEC[1]=0; } -// if ((n==2) && (pos[1]==0) && (pos[2]==0)) -// cerr << n << " -> " << pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[0] << endl; - //******************************* mu,sigma averaging *****************************// //shift down shiftCoord[n] = coord[n]-delta_M*0.25; @@ -306,7 +361,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) double delta_n = fabs(delta_M); if (n==1) { - delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius + delta_n = delta_n * fabs(coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius } if (prop) { @@ -330,7 +385,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) delta_n = fabs(delta); if (n==1) { - delta_n = delta_n * (coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius + delta_n = delta_n * fabs(coord[0]+0.5*fabs(deltaPP)); //multiply delta-angle by radius } if (prop) { @@ -360,7 +415,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) inEC[2] = gridDelta * A_n * 2 * __MUE0__ / inEC[2]; if (inEC[3]) inEC[3]=gridDelta * A_n * 2 / inEC[3]; -// if ((n==0) && (pos[1]==0) && (pos[2]==0)) +// if ((n==1) && (pos[1]==0) && (pos[2]==0)) // cerr << inEC[2]/(coord[0]) << endl; // cerr << n << " -> " << pos[0] << " " << pos[1] << " " << pos[2] << " " << inEC[2] << endl; @@ -369,28 +424,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, unsigned int* pos, double* inEC) bool Operator_Cylinder::Calc_EffMatPos(int n, unsigned int* pos, double* inMat) { - return false; //fixme - -// int nP = (n+1)%3; -// int nPP = (n+2)%3; -// -// unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); -// 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] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta; -// inMat[1] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta; -// -// 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; -// -// return true; + cerr << "Operator_Cylinder::Calc_EffMatPos: Warning! method not implemented yet..." << endl; + return false; } diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h index 8a98a20..62fa279 100644 --- a/FDTD/operator_cylinder.h +++ b/FDTD/operator_cylinder.h @@ -28,19 +28,24 @@ public: virtual bool SetGeometryCSX(ContinuousStructure* geo); + virtual int CalcECOperator(); + virtual void ApplyElectricBC(bool* dirs); virtual void ApplyMagneticBC(bool* dirs); - virtual void Reset(); - virtual unsigned int GetNumberOfLines(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, int* pos, bool dualMesh=false) const; + bool GetClosedAlpha() const {return CC_closedAlpha;} bool GetR0Included() const {return CC_R0_included;} protected: Operator_Cylinder(); virtual void Init(); + virtual void InitOperator(); + virtual void Reset(); virtual inline void Calc_ECOperatorPos(int n, unsigned int* pos); @@ -49,6 +54,12 @@ protected: virtual bool Calc_ECPos(int n, unsigned int* pos, double* inEC); virtual bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat); + +public: + //special EC operator for R0 + FDTD_FLOAT* vv_R0; //calc new voltage from old voltage + FDTD_FLOAT* vi_R0; //calc new voltage from old current + }; #endif // OPERATOR_CYLINDER_H diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp index 46d9c7d..5bb3382 100644 --- a/FDTD/processfields.cpp +++ b/FDTD/processfields.cpp @@ -128,7 +128,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) lines.clear(); for (unsigned int i=start[n];i<=stop[n];i+=subSample[n]) { - lines.push_back(Op->discLines[n][i]); + lines.push_back(Op->GetDiscLine(n,i));//Op->discLines[n][i]); } numLines[n] = lines.size(); delete[] discLines[n]; @@ -156,7 +156,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) lines.clear(); for (unsigned int i=start[n];idiscLines[n][i+1] + Op->discLines[n][i])); + lines.push_back(Op->GetDiscLine(n,i,true));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i])); } numDLines[n] = lines.size(); delete[] discDLines[n]; @@ -280,7 +280,7 @@ void ProcessFields::WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT* { for (pos[0]=0;pos[0]discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; - E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]+1]; - E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); + delta = Op->GetMeshDelta(0,OpPos); //Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; + if (delta) + { + E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]+1]; + E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); + } //in y - delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; - E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]+1]; - E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); + delta = Op->GetMeshDelta(1,OpPos); //Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; + if (delta) + { + E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]+1]; + E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); + } //in z - delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; - E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]+1][OpPos[2]]; - E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); + delta = Op->GetMeshDelta(2,OpPos); //Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; + if (delta) + { + E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]+1][OpPos[2]]; + E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); + } } } } @@ -104,18 +113,30 @@ void ProcessFieldsTD::DumpCellInterpol(string filename) { OpPos[2]=start[2]+pos[2]*subSample[2]; //in x - if (OpPos[0]==0) delta = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; - else delta = 0.5* (Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]-1]); - H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + curr[0][OpPos[0]+1][OpPos[1]][OpPos[2]]; - H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); +// if (OpPos[0]==0) delta = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; +// else delta = 0.5* (Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]-1]); + delta = Op->GetDiscLine(0,OpPos[0],true); + if (delta) + { + H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + curr[0][OpPos[0]+1][OpPos[1]][OpPos[2]]; + H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta); + } //in y - delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; - H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]] + curr[1][OpPos[0]][OpPos[1]+1][OpPos[2]]; - H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); +// delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; + delta = Op->GetDiscLine(1,OpPos[1],true); + if (delta) + { + H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]] + curr[1][OpPos[0]][OpPos[1]+1][OpPos[2]]; + H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta); + } //in z - delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; - H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]] + curr[2][OpPos[0]][OpPos[1]][OpPos[2]+1]; - H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); +// delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; + delta = Op->GetDiscLine(2,OpPos[2],true); + if (delta) + { + H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]] + curr[2][OpPos[0]][OpPos[1]][OpPos[2]+1]; + H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta); + } } } } @@ -154,18 +175,21 @@ void ProcessFieldsTD::DumpNoInterpol(string filename) for (pos[0]=0;pos[0]MainOp->GetIndexDelta(0,OpPos[0])); + delta[0]=Op->GetMeshDelta(0,OpPos);//fabs(Op->MainOp->GetIndexDelta(0,OpPos[0])); for (pos[1]=0;pos[1]MainOp->GetIndexDelta(1,OpPos[1])); + delta[1]=Op->GetMeshDelta(1,OpPos);//fabs(Op->MainOp->GetIndexDelta(1,OpPos[1])); for (pos[2]=0;pos[2]MainOp->GetIndexDelta(2,OpPos[2])); - E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0]/Op->gridDelta; - E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1]/Op->gridDelta; - E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2]/Op->gridDelta; + delta[2]=Op->GetMeshDelta(2,OpPos);//fabs(Op->MainOp->GetIndexDelta(2,OpPos[2])); + if (delta[0]) + E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0];// /Op->gridDelta; + if (delta[1]) + E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1];// /Op->gridDelta; + if (delta[2]) + E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2];// /Op->gridDelta; } } } @@ -196,19 +220,22 @@ void ProcessFieldsTD::DumpNoInterpol(string filename) for (pos[0]=0;pos[0]MainOp->GetIndexWidth(0,OpPos[0])); + delta[0]=Op->GetMeshDelta(0,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(0,OpPos[0])); for (pos[1]=0;pos[1]MainOp->GetIndexWidth(1,OpPos[1])); + delta[1]=Op->GetMeshDelta(1,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(1,OpPos[1])); for (pos[2]=0;pos[2]MainOp->GetIndexWidth(2,OpPos[2])); + delta[2]=Op->GetMeshDelta(2,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(2,OpPos[2])); //in x - H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0]/Op->gridDelta; - H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1]/Op->gridDelta; - H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2]/Op->gridDelta; + if (delta[0]) + H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0];// /Op->gridDelta; + if (delta[1]) + H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1];// /Op->gridDelta; + if (delta[2]) + H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2];// /Op->gridDelta; } } } diff --git a/openems.cpp b/openems.cpp index a156960..3143f52 100644 --- a/openems.cpp +++ b/openems.cpp @@ -237,6 +237,7 @@ int openEMS::SetupFDTD(const char* file) if (CylinderCoords) { FDTD_Op = Operator_Cylinder::New(); + CSX.SetCoordInputType(1); //tell CSX to use cylinder-coords } else {