Added support for background material in CSXCAD

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

Common/operator_base.cpp

@@ -57,6 +57,8 @@ void Operator_Base::Init()
 		discLines[n]=NULL;
 	for (int n=0; n<6; ++n)
 		m_BC[n]=0;
+
+	SetBackgroundMaterial(1,1,0,0);
 }
 
 void Operator_Base::Delete()
@@ -82,3 +84,52 @@ void Operator_Base::SetMaterialStoreFlags(int type, bool val)
 		return;
 	m_StoreMaterial[type]=val;
 }
+
+
+void Operator_Base::SetBackgroundMaterial(double epsR, double mueR, double kappa, double sigma)
+{
+	SetBackgroundEpsR(epsR);
+	SetBackgroundMueR(mueR);
+	SetBackgroundKappa(kappa);
+	SetBackgroundSigma(sigma);
+}
+
+void Operator_Base::SetBackgroundEpsR(double val)
+{
+	if (val<1)
+	{
+		cerr << __func__ << ": Warning, a relative electric permittivity <1 it not supported, skipping" << endl;
+		return;
+	}
+	m_BG_epsR=val;
+}
+
+void Operator_Base::SetBackgroundMueR(double val)
+{
+	if (val<1)
+	{
+		cerr << __func__ << ": Warning, a relative magnetic permeability <1 it not supported, skipping" << endl;
+		return;
+	}
+	m_BG_mueR=val;
+}
+
+void Operator_Base::SetBackgroundKappa(double val)
+{
+	if (val<0)
+	{
+		cerr << __func__ << ": Warning, an electric conductivity <0 it not supported, skipping" << endl;
+		return;
+	}
+	m_BG_kappa=val;
+}
+
+void Operator_Base::SetBackgroundSigma(double val)
+{
+	if (val<0)
+	{
+		cerr << __func__ << ": Warning, an artifival magnetic conductivity <0 it not supported, skipping" << endl;
+		return;
+	}
+	m_BG_sigma=val;
+}

Common/operator_base.h

@@ -108,6 +108,29 @@ public:
 	//! Get stored discrete material (if storage is enabled).
 	virtual double GetDiscMaterial(int type, int ny, const unsigned int pos[3]) const = 0;
 
+	//! Set the background material (default is vacuum)
+	virtual void SetBackgroundMaterial(double epsR=0, double mueR=0, double kappa=0, double sigma=0);
+
+	//! Get background rel. electric permittivity
+	double GetBackgroundEpsR() const {return m_BG_epsR;}
+	//! Set background rel. electric permittivity
+	void SetBackgroundEpsR(double val);
+
+	//! Get background rel. magnetic permeability
+	double GetBackgroundMueR() const {return m_BG_mueR;}
+	//! Set background rel. magnetic permeability
+	void SetBackgroundMueR(double val);
+
+	//! Get background electric conductivity
+	double GetBackgroundKappa() const {return m_BG_kappa;}
+	//! Set background electric conductivity
+	void SetBackgroundKappa(double val);
+
+	//! Get background magnetic conductivity (artificial)
+	double GetBackgroundSigma() const {return m_BG_sigma;}
+	//! Set background magnetic conductivity (artificial)
+	void SetBackgroundSigma(double val);
+
 protected:
 	Operator_Base();
 
@@ -127,6 +150,12 @@ protected:
 	//! bool flag array to store material data for post-processing
 	bool m_StoreMaterial[4];
 
+	//! background materials
+	double m_BG_epsR;
+	double m_BG_mueR;
+	double m_BG_kappa;
+	double m_BG_sigma;
+
 	CoordinateSystem m_MeshType;
 	unsigned int numLines[3];
 	double* discLines[3];

FDTD/extensions/operator_ext_mur_abc.cpp

@@ -175,7 +175,7 @@ bool Operator_Ext_Mur_ABC::BuildExtension()
 				if (m_v_phase>0.0)
 					c0t = m_v_phase * dT;
 				else
-					c0t = __C0__ * dT;
+					c0t = __C0__ / sqrt(m_Op->GetBackgroundEpsR()*m_Op->GetBackgroundMueR()) * dT;
 				m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta);
 				m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]];
 			}

FDTD/extensions/operator_ext_tfsf.cpp

@@ -100,7 +100,8 @@ bool Operator_Ext_TFSF::BuildExtension()
 		return false;
 	}
 
-	m_PhVel = __C0__;
+	double ref_index = sqrt(m_Op->GetBackgroundEpsR()*m_Op->GetBackgroundMueR());
+	m_PhVel = __C0__/ref_index;
 	CSPropExcitation* elec=NULL;
 	CSProperties* prop=NULL;
 	CSPrimitives* prim=NULL;
@@ -154,14 +155,14 @@ bool Operator_Ext_TFSF::BuildExtension()
 //		if (m_Frequency<=0)
 //			m_Frequency = m_Op->GetExcitationSignal()->GetFrequencyOfInterest();
 		if (m_Frequency<=0)
-			m_PhVel=__C0__;
+			m_PhVel=__C0__/ref_index;
 		else
 			m_PhVel=m_Op->CalcNumericPhaseVelocity(m_Start,m_Stop,m_PropDir,m_Frequency);
 
-		if ((m_PhVel<0) || (m_PhVel>__C0__) || isnan(m_PhVel))
+		if ((m_PhVel<0) || (m_PhVel>__C0__/ref_index) || isnan(m_PhVel))
 		{
 			cerr << "Operator_Ext_TFSF::BuildExtension: Warning, invalid phase velocity found, resetting to c0! " << endl;
-			m_PhVel = __C0__;
+			m_PhVel = __C0__/ref_index;
 		}
 
 		double origin[3];
@@ -215,7 +216,7 @@ bool Operator_Ext_TFSF::BuildExtension()
 			nP = (n+1)%3;
 			nPP = (n+2)%3;
 			m_H_Amp[n] = m_PropDir[nP]*m_E_Amp[nPP] - m_PropDir[nPP]*m_E_Amp[nP];
-			m_H_Amp[n] /= __Z0__;
+			m_H_Amp[n] /= __Z0__*sqrt(m_Op->GetBackgroundMueR()/m_Op->GetBackgroundEpsR());
 		}
 
 		double coord[3];
@@ -419,8 +420,8 @@ void Operator_Ext_TFSF::ShowStat(ostream &ostr) const
 	cout << "Active directions\t: " << "(" << m_ActiveDir[0][0] << "/" << m_ActiveDir[0][1] << ", " << m_ActiveDir[1][0] << "/" << m_ActiveDir[1][1]  << ", " << m_ActiveDir[2][0] << "/" << m_ActiveDir[2][1]  << ")" << endl;
 	cout << "Propagation direction\t: " << "(" << m_PropDir[0] << ", " << m_PropDir[1] << ", " << m_PropDir[2] << ")" << endl;
 	cout << "Rel. propagation speed\t: " << m_PhVel/__C0__ << "*c0  @ " << m_Frequency << " Hz" << endl;
-	cout << "E-field amplitude\t: " << "(" << m_E_Amp[0] << ", " << m_E_Amp[1] << ", " << m_E_Amp[2] << ")" << endl;
+	cout << "E-field amplitude (V/m)\t: " << "(" << m_E_Amp[0] << ", " << m_E_Amp[1] << ", " << m_E_Amp[2] << ")" << endl;
-	cout << "H-field amplitude\t: " << "(" << m_H_Amp[0]*__Z0__ << ", " << m_H_Amp[1]*__Z0__ << ", " << m_H_Amp[2]*__Z0__ << ")/Z0" << endl;
+	cout << "H-field amplitude (A/m)\t: " << "(" << m_H_Amp[0] << ", " << m_H_Amp[1] << ", " << m_H_Amp[2] << ")" << endl;
 	cout << "Box Dimensions\t\t: " << m_numLines[0] << " x " << m_numLines[1] << " x " << m_numLines[2] << endl;
 	cout << "Max. Delay (TS)\t\t: " << m_maxDelay << endl;
 	int dirs = m_ActiveDir[0][0] + m_ActiveDir[0][1] + m_ActiveDir[1][0] + m_ActiveDir[1][1] + m_ActiveDir[2][0] + m_ActiveDir[2][1] ;

FDTD/operator.cpp

@@ -505,6 +505,8 @@ void Operator::ShowStat() const
 	cout << "Size of Operator\t: " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MiB) " << endl;
 	cout << "Size of Field-Data\t: " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MiB) " << endl;
 	cout << "-----------------------------------" << endl;
+	cout << "Background materials (epsR/mueR/kappa/sigma): " << GetBackgroundEpsR() << "/" << GetBackgroundMueR() << "/" << GetBackgroundKappa() << "/" << GetBackgroundSigma() << endl;
+	cout << "-----------------------------------" << endl;
 	cout << "Number of PEC edges\t: " << m_Nr_PEC[0]+m_Nr_PEC[1]+m_Nr_PEC[2] << endl;
 	cout << "in " << GetDirName(0) << " direction\t\t: " << m_Nr_PEC[0] << endl;
 	cout << "in " << GetDirName(1) << " direction\t\t: " << m_Nr_PEC[1] << endl;
@@ -805,8 +807,13 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo)
 
 	CSX = geo;
 
-	CSRectGrid* grid=CSX->GetGrid();
+	CSBackgroundMaterial* bg_mat=CSX->GetBackgroundMaterial();
+	SetBackgroundEpsR(bg_mat->GetEpsilon());
+	SetBackgroundMueR(bg_mat->GetMue());
+	SetBackgroundKappa(bg_mat->GetKappa());
+	SetBackgroundSigma(bg_mat->GetSigma());
 
+	CSRectGrid* grid=CSX->GetGrid();
 	return SetupCSXGrid(CSRectGrid::Clone(grid));
 }
 
@@ -1227,13 +1234,13 @@ double Operator::GetMaterial(int ny, const double* coords, int MatType, bool mar
 	switch (MatType)
 	{
 	case 0:
-		return 1;
+		return GetBackgroundEpsR();
 	case 1:
-		return 0;
+		return GetBackgroundKappa();
 	case 2:
-		return 1;
+		return GetBackgroundMueR();
 	case 3:
-		return 0;
+		return GetBackgroundSigma();
 	default:
 		cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
 		return 0;
@@ -1850,6 +1857,7 @@ void Operator::DeleteExtension(Operator_Extension* op_ext)
 double Operator::CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const
 {
 	double average_mesh_disc[3];
+	double c0 = __C0__/sqrt(GetBackgroundEpsR()*GetBackgroundMueR());
 
 	//calculate average mesh deltas
 	for (int n=0;n<3;++n)
@@ -1864,18 +1872,18 @@ double Operator::CalcNumericPhaseVelocity(unsigned int start[3], unsigned int st
 		int nPP = (n+2)%3;
 		if ((fabs(propDir[n])==1) && (propDir[nP]==0) && (propDir[nPP]==0))
 		{
-			double kx = asin(average_mesh_disc[0]/__C0__/dT*sin(2*PI*freq*dT/2))*2/average_mesh_disc[0];
+			double kx = asin(average_mesh_disc[0]/c0/dT*sin(2*PI*freq*dT/2))*2/average_mesh_disc[0];
 			return 2*PI*freq/kx;
 		}
 	}
 
 	// else, do an newton iterative estimation
-	double k0=2*PI*freq/__C0__;
+	double k0=2*PI*freq/c0;
 	double k=k0;
-	double RHS = pow(sin(2*PI*freq*dT/2)/__C0__/dT,2);
+	double RHS = pow(sin(2*PI*freq*dT/2)/c0/dT,2);
 	double fk=1,fdk=0;
 	double old_phv=0;
-	double phv=__C0__;
+	double phv=c0;
 	double err_est = 1e-6;
 	int it_count=0;
 	while (fabs(old_phv-phv)>err_est)

matlab/CalcNF2FF.m

@@ -44,6 +44,13 @@ nf2ff_xml.Planes = {};
 
 nf2ff_xml.ATTRIBUTE.Outfile = [filename '.h5'];
 
+if (isfield(nf2ff,'Eps_r'))
+    nf2ff_xml.ATTRIBUTE.Eps_r = nf2ff.Eps_r;
+end
+if (isfield(nf2ff,'Mue_r'))
+    nf2ff_xml.ATTRIBUTE.Mue_r = nf2ff.Mue_r;
+end
+
 for n=1:2:numel(varargin)-1
     if (strcmp(varargin{n},'Mode'))
         mode = varargin{n+1};

matlab/CreateNF2FFBox.m

@@ -57,6 +57,14 @@ if (isfield(CSX,'ATTRIBUTE'))
     if (isfield(CSX.ATTRIBUTE,'CoordSystem'))
         nf2ff.CoordSystem = CSX.ATTRIBUTE.CoordSystem;
     end
+    if (isfield(CSX,'BackgroundMaterial'))
+        if (isfield(CSX.ATTRIBUTE,'Epsilon'))
+            nf2ff.Eps_r = CSX.ATTRIBUTE.BG_epsR;
+        end
+        if (isfield(CSX.ATTRIBUTE,'Mue'))
+            nf2ff.Mue_r = CSX.ATTRIBUTE.BG_mueR;
+        end
+    end
 end
 
 for nd = 1:3

nf2ff/nf2ff.cpp

@@ -84,12 +84,19 @@ void nf2ff::SetPermittivity(vector<float> permittivity)
 	if (permittivity.size()==0)
 		return;
 
+	m_permittivity = permittivity;
+	if (permittivity.size()==1)
+	{
+		for (size_t fn=0;fn<m_nf2ff.size();++fn)
+			m_nf2ff.at(fn)->SetPermittivity(permittivity.at(0));
+		return;
+	}
+
 	if (permittivity.size()!=m_freq.size())
 	{
 		cerr << __func__ << ": Error, permittivity vector size must match number of set frequencies! skipping!" << endl;
 		return;
 	}
-	m_permittivity = permittivity;
 	for (size_t fn=0;fn<m_nf2ff.size();++fn)
 		m_nf2ff.at(fn)->SetPermittivity(permittivity.at(fn));
 }
@@ -99,12 +106,19 @@ void nf2ff::SetPermeability(vector<float> permeability)
 	if (permeability.size()==0)
 		return;
 
+	m_permeability = permeability;
+	if (permeability.size()==1)
+	{
+		for (size_t fn=0;fn<m_nf2ff.size();++fn)
+			m_nf2ff.at(fn)->SetPermeability(permeability.at(0));
+		return;
+	}
+
 	if (permeability.size()!=m_freq.size())
 	{
 		cerr << __func__ << ": Error, permeability vector size must match number of set frequencies! skipping!" << endl;
 		return;
 	}
-	m_permeability = permeability;
 	for (size_t fn=0;fn<m_nf2ff.size();++fn)
 		m_nf2ff.at(fn)->SetPermeability(permeability.at(fn));