zcy
/
openEMS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

new: support debye dispersive material 

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
pull/1/head
Thorsten Liebig 2013-03-19 14:02:06 +01:00
parent 57dfd64c9b
commit e113afb656
2 changed files with 41 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
FDTD/extensions/operator_ext_lorentzmaterial.cpp
View File

@ -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<unsigned int> 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.size();++n)
{
CSPropDebyeMaterial* DebyeMat = dynamic_cast<CSPropDebyeMaterial*>(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
{

2
openems.cpp
View File

@ -674,7 +674,7 @@ int openEMS::SetupFDTD(const char* file)
if (FDTD_Opts->QueryDoubleAttribute("TimeStepFactor",&timestepfactor)==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()));