diff --git a/Common/operator_base.cpp b/Common/operator_base.cpp index 3271008..bea1581 100644 --- a/Common/operator_base.cpp +++ b/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; +} diff --git a/Common/operator_base.h b/Common/operator_base.h index 9da358e..9a9fdfa 100644 --- a/Common/operator_base.h +++ b/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]; diff --git a/FDTD/extensions/operator_ext_mur_abc.cpp b/FDTD/extensions/operator_ext_mur_abc.cpp index a9697d2..19a6389 100644 --- a/FDTD/extensions/operator_ext_mur_abc.cpp +++ b/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]]; } diff --git a/FDTD/extensions/operator_ext_tfsf.cpp b/FDTD/extensions/operator_ext_tfsf.cpp index e57e16d..663dc1f 100644 --- a/FDTD/extensions/operator_ext_tfsf.cpp +++ b/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 << "H-field amplitude\t: " << "(" << m_H_Amp[0]*__Z0__ << ", " << m_H_Amp[1]*__Z0__ << ", " << m_H_Amp[2]*__Z0__ << ")/Z0" << 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 (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] ; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 71d11d4..2ab5dcd 100644 --- a/FDTD/operator.cpp +++ b/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) diff --git a/matlab/CalcNF2FF.m b/matlab/CalcNF2FF.m index a6fb555..2e1fc37 100644 --- a/matlab/CalcNF2FF.m +++ b/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}; diff --git a/matlab/CreateNF2FFBox.m b/matlab/CreateNF2FFBox.m index 1a87b37..a97f0d2 100644 --- a/matlab/CreateNF2FFBox.m +++ b/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 diff --git a/nf2ff/nf2ff.cpp b/nf2ff/nf2ff.cpp index 3c48f5e..6a6fae9 100644 --- a/nf2ff/nf2ff.cpp +++ b/nf2ff/nf2ff.cpp @@ -84,12 +84,19 @@ void nf2ff::SetPermittivity(vector permittivity) if (permittivity.size()==0) return; + m_permittivity = permittivity; + if (permittivity.size()==1) + { + for (size_t fn=0;fnSetPermittivity(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;fnSetPermittivity(permittivity.at(fn)); } @@ -99,12 +106,19 @@ void nf2ff::SetPermeability(vector permeability) if (permeability.size()==0) return; + m_permeability = permeability; + if (permeability.size()==1) + { + for (size_t fn=0;fnSetPermeability(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;fnSetPermeability(permeability.at(fn));