From e113afb6566823caad8470b71ac7c5ae1aa5d8a9 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Tue, 19 Mar 2013 14:02:06 +0100 Subject: [PATCH] new: support debye dispersive material Signed-off-by: Thorsten Liebig --- .../operator_ext_lorentzmaterial.cpp | 45 ++++++++++++++++--- openems.cpp | 2 +- 2 files changed, 41 insertions(+), 6 deletions(-) diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/FDTD/extensions/operator_ext_lorentzmaterial.cpp index 3b303b8..7ef355c 100644 --- a/FDTD/extensions/operator_ext_lorentzmaterial.cpp +++ b/FDTD/extensions/operator_ext_lorentzmaterial.cpp @@ -21,6 +21,7 @@ #include "../operator_cylinder.h" #include "CSPropLorentzMaterial.h" +#include "CSPropDebyeMaterial.h" Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Operator_Ext_Dispersive(op) { @@ -120,6 +121,9 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() double coord[3]; unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)}; CSPropLorentzMaterial* mat = NULL; + CSPropDebyeMaterial* debye_mat = NULL; + + bool warn_once = true; bool b_pos_on; vector v_pos[3]; @@ -150,6 +154,15 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() if (LorMat->GetDispersionOrder()>m_Order) m_Order=LorMat->GetDispersionOrder(); } + LD_props = m_Op->CSX->GetPropertyByType(CSProperties::DEBYEMATERIAL); + for (size_t n=0;n(LD_props.at(n)); + if (DebyeMat==NULL) + return false; //sanity check, this should not happen + if (DebyeMat->GetDispersionOrder()>m_Order) + m_Order=DebyeMat->GetDispersionOrder(); + } m_LM_pos = new unsigned int**[m_Order]; @@ -231,6 +244,23 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() C_L[n] = 1/(w_Lor_Pol*w_Lor_Pol*L_D[n]); } } + if ((debye_mat = prop->ToDebyeMaterial())) + { + C_L[n] = 8.85418781762e-12*debye_mat->GetEpsDeltaWeighted(order,n,coord) * m_Op->GetEdgeArea(n, pos) / m_Op->GetEdgeLength(n,pos); + t_relax = debye_mat->GetEpsRelaxTimeWeighted(order,n,coord); + if ((t_relax<2.0*dT) && warn_once) + { + warn_once = false; + cerr << "Operator_Ext_LorentzMaterial::BuildExtension(): Warning, debye relaxation time is to small, skipping..." << endl; + } + if ((C_L[n]>0) && (t_relax>0) && (t_relax>2.0*dT)) + { + R_D[n] = t_relax/C_L[n]; + b_pos_on = true; + m_volt_ADE_On[order] = true; + m_volt_Lor_ADE_On[order] = true; + } + } } for (int n=0; n<3; ++n) @@ -274,12 +304,17 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() v_pos[n].push_back(pos[n]); 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_int[n].push_back((2.0*L_D[n]-dT*R_D[n])/(2.0*L_D[n]+dT*R_D[n])); // 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]]); + v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2.0)*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])); + v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2.0)*m_Op->GetVI(n,pos[0],pos[1],pos[2])); + } + else if ((R_D[n]>0) && (C_L[n]>0)) + { + v_int[n].push_back((2.0*dT-R_D[n]*C_L[n])/(C_L[n]*R_D[n])); + v_ext[n].push_back(2.0/R_D[n]*m_Op->GetVI(n,pos[0],pos[1],pos[2])); } else { @@ -288,8 +323,8 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() } if (C_D[n]>0) { - i_int[n].push_back((2*C_D[n]-dT*G_D[n])/(2*C_D[n]+dT*G_D[n])); - i_ext[n].push_back(dT/(C_D[n]+dT*G_D[n]/2)*m_Op->GetIV(n,pos[0],pos[1],pos[2])); + i_int[n].push_back((2.0*C_D[n]-dT*G_D[n])/(2.0*C_D[n]+dT*G_D[n])); + i_ext[n].push_back(dT/(C_D[n]+dT*G_D[n]/2.0)*m_Op->GetIV(n,pos[0],pos[1],pos[2])); } else { diff --git a/openems.cpp b/openems.cpp index a34e00a..75ec38d 100644 --- a/openems.cpp +++ b/openems.cpp @@ -674,7 +674,7 @@ int openEMS::SetupFDTD(const char* file) if (FDTD_Opts->QueryDoubleAttribute("TimeStepFactor",×tepfactor)==TIXML_SUCCESS) FDTD_Op->SetTimestepFactor(timestepfactor); - if (m_CSX->GetQtyPropertyType(CSProperties::LORENTZMATERIAL)>0) + if ((m_CSX->GetQtyPropertyType(CSProperties::LORENTZMATERIAL)>0) || (m_CSX->GetQtyPropertyType(CSProperties::DEBYEMATERIAL)>0)) FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op)); if (m_CSX->GetQtyPropertyType(CSProperties::CONDUCTINGSHEET)>0) FDTD_Op->AddExtension(new Operator_Ext_ConductingSheet(FDTD_Op,m_Exc->GetMaxFrequency()));