From 138624d2d8efff02d53f42e5fb60f9af50250919 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Wed, 28 Apr 2010 21:45:05 +0200 Subject: [PATCH] mur_abc_extension: now allows material distribution --- FDTD/engine_ext_mur_abc.cpp | 14 +++---- FDTD/engine_ext_mur_abc.h | 3 +- FDTD/operator.h | 1 + FDTD/operator_ext_mur_abc.cpp | 69 +++++++++++++++++++++++++++++++++-- FDTD/operator_ext_mur_abc.h | 3 +- 5 files changed, 76 insertions(+), 14 deletions(-) diff --git a/FDTD/engine_ext_mur_abc.cpp b/FDTD/engine_ext_mur_abc.cpp index 17359b2..a5c3c19 100644 --- a/FDTD/engine_ext_mur_abc.cpp +++ b/FDTD/engine_ext_mur_abc.cpp @@ -31,7 +31,8 @@ Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext) : Engine_Ex m_LineNr = m_Op_mur->m_LineNr; m_LineNr_Shift = m_Op_mur->m_LineNr_Shift; - m_Mur_Coeff = m_Op_mur->m_Mur_Coeff; + m_Mur_Coeff_nyP = m_Op_mur->m_Mur_Coeff_nyP; + m_Mur_Coeff_nyPP = m_Op_mur->m_Mur_Coeff_nyPP; m_volt_nyP = Create2DArray(m_numLines); m_volt_nyPP = Create2DArray(m_numLines); @@ -66,7 +67,6 @@ void Engine_Ext_Mur_ABC::DoPreVoltageUpdates() { if (IsActive()==false) return; if (m_Eng==NULL) return; - if (m_Mur_Coeff==0) return; unsigned int pos[] = {0,0,0}; unsigned int pos_shift[] = {0,0,0}; pos[m_ny] = m_LineNr; @@ -78,8 +78,8 @@ void Engine_Ext_Mur_ABC::DoPreVoltageUpdates() for (pos[m_nyPP]=0;pos[m_nyPP]GetVolt(m_nyP,pos_shift) - m_Mur_Coeff * m_Eng->GetVolt(m_nyP,pos); - m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyPP,pos_shift) - m_Mur_Coeff * m_Eng->GetVolt(m_nyPP,pos); + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyPP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyPP,pos); } } } @@ -88,7 +88,6 @@ void Engine_Ext_Mur_ABC::DoPostVoltageUpdates() { if (IsActive()==false) return; if (m_Eng==NULL) return; - if (m_Mur_Coeff==0) return; unsigned int pos[] = {0,0,0}; unsigned int pos_shift[] = {0,0,0}; pos[m_ny] = m_LineNr; @@ -100,8 +99,8 @@ void Engine_Ext_Mur_ABC::DoPostVoltageUpdates() for (pos[m_nyPP]=0;pos[m_nyPP]GetVolt(m_nyP,pos_shift); - m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Mur_Coeff * m_Eng->GetVolt(m_nyPP,pos_shift); + m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos_shift); + m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyPP,pos_shift); } } } @@ -110,7 +109,6 @@ void Engine_Ext_Mur_ABC::Apply2Voltages() { if (IsActive()==false) return; if (m_Eng==NULL) return; - if (m_Mur_Coeff==0) return; unsigned int pos[] = {0,0,0}; pos[m_ny] = m_LineNr; diff --git a/FDTD/engine_ext_mur_abc.h b/FDTD/engine_ext_mur_abc.h index 9fb540f..930bc82 100644 --- a/FDTD/engine_ext_mur_abc.h +++ b/FDTD/engine_ext_mur_abc.h @@ -46,7 +46,8 @@ protected: int m_LineNr_Shift; unsigned int m_numLines[2]; - FDTD_FLOAT m_Mur_Coeff; + FDTD_FLOAT** m_Mur_Coeff_nyP; + FDTD_FLOAT** m_Mur_Coeff_nyPP; FDTD_FLOAT** m_volt_nyP; //n+1 direction FDTD_FLOAT** m_volt_nyPP; //n+2 direction }; diff --git a/FDTD/operator.h b/FDTD/operator.h index 851d08c..120439d 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -35,6 +35,7 @@ public: virtual ~Operator(); virtual bool SetGeometryCSX(ContinuousStructure* geo); + virtual ContinuousStructure* GetGeometryCSX() {return CSX;} virtual int CalcECOperator(); diff --git a/FDTD/operator_ext_mur_abc.cpp b/FDTD/operator_ext_mur_abc.cpp index e0ebd5a..29dde53 100644 --- a/FDTD/operator_ext_mur_abc.cpp +++ b/FDTD/operator_ext_mur_abc.cpp @@ -22,23 +22,35 @@ Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op) { - m_Mur_Coeff = NULL; m_ny = -1; m_nyP = -1; m_nyPP = -1; m_LineNr = 0; m_LineNr_Shift = 0; - m_Mur_Coeff = 0; + + m_Mur_Coeff_nyP = NULL; + m_Mur_Coeff_nyPP = NULL; + + m_numLines[0]=0; + m_numLines[1]=0; } Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC() { + Delete2DArray(m_Mur_Coeff_nyP,m_numLines); + m_Mur_Coeff_nyP = NULL; + Delete2DArray(m_Mur_Coeff_nyPP,m_numLines); + m_Mur_Coeff_nyPP = NULL; } void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny) { if ((ny<0) || (ny>2)) return; + + Delete2DArray(m_Mur_Coeff_nyP,m_numLines); + Delete2DArray(m_Mur_Coeff_nyPP,m_numLines); + m_ny = ny; m_nyP = (ny+1)%3; m_nyPP = (ny+2)%3; @@ -55,6 +67,10 @@ void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny) m_numLines[0] = m_Op->GetNumberOfLines(m_nyP); m_numLines[1] = m_Op->GetNumberOfLines(m_nyPP); + + m_Mur_Coeff_nyP = Create2DArray(m_numLines); + m_Mur_Coeff_nyPP = Create2DArray(m_numLines); + } bool Operator_Ext_Mur_ABC::BuildExtension() @@ -68,8 +84,53 @@ bool Operator_Ext_Mur_ABC::BuildExtension() unsigned int pos[] = {0,0,0}; pos[m_ny] = m_LineNr; double delta = fabs(m_Op->GetMeshDelta(m_ny,pos)); - m_Mur_Coeff = (__C0__ * dT - delta) / (__C0__ * dT + delta); -// cerr << "Operator_Ext_Mur_ABC::BuildExtension(): " << m_ny << " @ " << m_LineNr << " = "<< m_Mur_Coeff << endl; + double coord[] = {0,0,0}; + coord[0] = m_Op->GetDiscLine(0,pos[0]); + coord[1] = m_Op->GetDiscLine(1,pos[1]); + coord[2] = m_Op->GetDiscLine(2,pos[2]); + + double eps,mue; + double c0t; + + if (m_LineNr==0) + coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) + delta/2 / m_Op->GetGridDelta(); + else + coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) - delta/2 / m_Op->GetGridDelta(); + + for (pos[m_nyP]=0;pos[m_nyP]GetDiscLine(m_nyP,pos[m_nyP]); + for (pos[m_nyPP]=0;pos[m_nyPP]GetDiscLine(m_nyPP,pos[m_nyPP]); + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); + if (prop) + { + CSPropMaterial* mat = prop->ToMaterial(); + + //nP + eps = mat->GetEpsilonWeighted(m_nyP,coord); + mue = mat->GetMueWeighted(m_nyP,coord); + c0t = __C0__ * dT / sqrt(eps*mue); + m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta); + + //nPP + eps = mat->GetEpsilonWeighted(m_nyPP,coord); + mue = mat->GetMueWeighted(m_nyPP,coord); + c0t = __C0__ * dT / sqrt(eps*mue); + m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta); + + } + else + { + c0t = __C0__ * 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]]; + } +// cerr << m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] << " : " << m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] << endl; + } + } +// cerr << "Operator_Ext_Mur_ABC::BuildExtension(): " << m_ny << " @ " << m_LineNr << endl; return true; } diff --git a/FDTD/operator_ext_mur_abc.h b/FDTD/operator_ext_mur_abc.h index 52b15d1..bcef302 100644 --- a/FDTD/operator_ext_mur_abc.h +++ b/FDTD/operator_ext_mur_abc.h @@ -43,7 +43,8 @@ protected: unsigned int m_numLines[2]; - FDTD_FLOAT m_Mur_Coeff; + FDTD_FLOAT** m_Mur_Coeff_nyP; + FDTD_FLOAT** m_Mur_Coeff_nyPP; }; #endif // OPERATOR_EXT_MUR_ABC_H