diff --git a/FDTD/extensions/operator_ext_cylinder.cpp b/FDTD/extensions/operator_ext_cylinder.cpp index 9e336a4..dead36b 100644 --- a/FDTD/extensions/operator_ext_cylinder.cpp +++ b/FDTD/extensions/operator_ext_cylinder.cpp @@ -72,6 +72,12 @@ bool Operator_Ext_Cylinder::BuildExtension() vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C); vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C); + for (unsigned int i=0; iGetOriginalNumLines(1); ++i) + { + m_Op->EC_C[2][m_Op->MainOp->SetPos(0,i,pos[2])] = C; + m_Op->EC_G[2][m_Op->MainOp->SetPos(0,i,pos[2])] = G; + } + //search for metal on z-axis m_Op_Cyl->GetYeeCoords(2,pos,coord,false); CSProperties* prop = m_Op->CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL), true); @@ -82,6 +88,8 @@ bool Operator_Ext_Cylinder::BuildExtension() m_Op->SetVV(2,0,0,pos[2], 0); vv_R0[pos[2]] = 0; vi_R0[pos[2]] = 0; + m_Op->EC_C[2][m_Op->MainOp->SetPos(0,0,pos[2])] = 0; + m_Op->EC_G[2][m_Op->MainOp->SetPos(0,0,pos[2])] = 0; } } } diff --git a/FDTD/extensions/operator_ext_cylinder.h b/FDTD/extensions/operator_ext_cylinder.h index cc9919d..c47a6bf 100644 --- a/FDTD/extensions/operator_ext_cylinder.h +++ b/FDTD/extensions/operator_ext_cylinder.h @@ -26,6 +26,7 @@ class Operator_Cylinder; class Operator_Ext_Cylinder : public Operator_Extension { friend class Engine_Ext_Cylinder; + friend class Operator_Ext_LorentzMaterial; public: Operator_Ext_Cylinder(Operator_Cylinder* op); ~Operator_Ext_Cylinder(); diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/FDTD/extensions/operator_ext_lorentzmaterial.cpp index d2b93ad..e536f7b 100644 --- a/FDTD/extensions/operator_ext_lorentzmaterial.cpp +++ b/FDTD/extensions/operator_ext_lorentzmaterial.cpp @@ -17,6 +17,8 @@ #include "operator_ext_lorentzmaterial.h" #include "engine_ext_lorentzmaterial.h" +#include "operator_ext_cylinder.h" +#include "../operator_cylinder.h" Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Operator_Ext_Dispersive(op) { @@ -69,7 +71,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() double dT = m_Op->GetTimestep(); unsigned int pos[] = {0,0,0}; double coord[3]; - unsigned int numLines[3] = {m_Op->GetNumberOfLines(0),m_Op->GetNumberOfLines(1),m_Op->GetNumberOfLines(2)}; + unsigned int numLines[3] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)}; CSPropLorentzMaterial* mat = NULL; double w_plasma,t_relax; @@ -129,16 +131,16 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() { L_D[n]=0; R_D[n]=0; - coord[0] = m_Op->GetDiscLine(0,pos[0]); - coord[1] = m_Op->GetDiscLine(1,pos[1]); - coord[2] = m_Op->GetDiscLine(2,pos[2]); - coord[n] = m_Op->GetDiscLine(n,pos[n],true); //pos of E_n + if (m_Op->GetYeeCoords(n,pos,coord,false)==false) + continue; + if (m_CC_R0_included && (n==2) && (pos[0]==0)) + coord[1] = m_Op->GetDiscLine(1,0); CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,CSProperties::LORENTZMATERIAL, true); if ((mat = prop->ToLorentzMaterial())) { w_plasma = mat->GetEpsPlasmaFreqWeighted(order,n,coord) * 2 * PI; - if (w_plasma>0) + if ((w_plasma>0) && (m_Op->EC_C[n][index]>0)) { b_pos_on = true; m_volt_ADE_On[order] = true; @@ -156,16 +158,14 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() { C_D[n]=0; G_D[n]=0; - coord[0] = m_Op->GetDiscLine(0,pos[0],true); - coord[1] = m_Op->GetDiscLine(1,pos[1],true); - coord[2] = m_Op->GetDiscLine(2,pos[2],true); - coord[n] = m_Op->GetDiscLine(n,pos[n]); //pos of H_n + if (m_Op->GetYeeCoords(n,pos,coord,true)==false) + continue; CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,CSProperties::LORENTZMATERIAL, true); if ((mat = prop->ToLorentzMaterial())) { w_plasma = mat->GetMuePlasmaFreqWeighted(order,n,coord) * 2 * PI; - if (w_plasma>0) + if ((w_plasma>0) && (m_Op->EC_L[n][index]>0)) { b_pos_on = true; m_curr_ADE_On[order] = true; @@ -187,7 +187,11 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() if (L_D[n]>0) { v_int[n].push_back((2*L_D[n]-dT*R_D[n])/(2*L_D[n]+dT*R_D[n])); - v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2)*m_Op->GetVI(n,pos[0],pos[1],pos[2])); + // check for r==0 in clyindrical coords and get special VI cooefficient + if (m_CC_R0_included && n==2 && pos[0]==0) + v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2)*m_Op_Cyl->m_Cyl_Ext->vi_R0[pos[2]]); + else + v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2)*m_Op->GetVI(n,pos[0],pos[1],pos[2])); } else { diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.h b/FDTD/extensions/operator_ext_lorentzmaterial.h index 9c34002..ee4386c 100644 --- a/FDTD/extensions/operator_ext_lorentzmaterial.h +++ b/FDTD/extensions/operator_ext_lorentzmaterial.h @@ -32,7 +32,7 @@ public: virtual Engine_Extension* CreateEngineExtention(); - virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {return false;} + virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {return true;} virtual string GetExtensionName() const {return string("Drude/Lorentz Dispersive Material Extension");} diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index af89ece..0b5eefd 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -34,6 +34,7 @@ Operator_Cylinder* Operator_Cylinder::New(unsigned int numThreads) Operator_Cylinder::Operator_Cylinder() : Operator_Multithread() { m_MeshType = CYLINDRICAL; + m_Cyl_Ext = NULL; } Operator_Cylinder::~Operator_Cylinder() @@ -280,7 +281,10 @@ bool Operator_Cylinder::SetupCSXGrid(CSRectGrid* grid) return false; if (CC_closedAlpha || CC_R0_included) - this->AddExtension(new Operator_Ext_Cylinder(this)); + { + m_Cyl_Ext = new Operator_Ext_Cylinder(this); + this->AddExtension(m_Cyl_Ext); + } return true; } diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h index f3c2420..fd946e7 100644 --- a/FDTD/operator_cylinder.h +++ b/FDTD/operator_cylinder.h @@ -22,6 +22,8 @@ #include "operator_multithread.h" +class Operator_Ext_Cylinder; + //! This class creates an operator for a cylindrical FDTD. /*! This class creates an operator for a cylindrical FDTD. No special engine is necessary, @@ -31,6 +33,7 @@ class Operator_Cylinder : public Operator_Multithread { friend class Operator_CylinderMultiGrid; friend class Operator_Ext_Cylinder; + friend class Operator_Ext_LorentzMaterial; public: static Operator_Cylinder* New(unsigned int numThreads = 0); virtual ~Operator_Cylinder(); @@ -89,6 +92,7 @@ protected: bool CC_closedAlpha; bool CC_R0_included; + Operator_Ext_Cylinder* m_Cyl_Ext; #ifdef MPI_SUPPORT bool CC_MPI_Alpha;