lorentz material model now supports cylindrical coordinates

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; i<m_Op->GetOriginalNumLines(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;
 			}
 		}
 	}

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();

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]);
+						if (m_Op->GetYeeCoords(n,pos,coord,false)==false)
-						coord[1] = m_Op->GetDiscLine(1,pos[1]);
+							continue;
-						coord[2] = m_Op->GetDiscLine(2,pos[2]);
+						if (m_CC_R0_included && (n==2) && (pos[0]==0))
-						coord[n] = m_Op->GetDiscLine(n,pos[n],true); //pos of E_n
+							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);
+						if (m_Op->GetYeeCoords(n,pos,coord,true)==false)
-						coord[1] = m_Op->GetDiscLine(1,pos[1],true);
+							continue;
-						coord[2] = m_Op->GetDiscLine(2,pos[2],true);
-						coord[n] = m_Op->GetDiscLine(n,pos[n]); //pos of H_n
 
 						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
 							{

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");}
 

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;
 }

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;