From 57dfd64c9b24ff80f2baedbf1593bf27848b0136 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 14 Mar 2013 16:52:44 +0100 Subject: [PATCH] added support for Lorentz dispersive material Signed-off-by: Thorsten Liebig --- .../extensions/engine_ext_lorentzmaterial.cpp | 312 ++++++++++++++---- FDTD/extensions/engine_ext_lorentzmaterial.h | 9 + .../operator_ext_lorentzmaterial.cpp | 125 ++++++- .../extensions/operator_ext_lorentzmaterial.h | 6 + 4 files changed, 375 insertions(+), 77 deletions(-) diff --git a/FDTD/extensions/engine_ext_lorentzmaterial.cpp b/FDTD/extensions/engine_ext_lorentzmaterial.cpp index 783829d..9aaa1ff 100644 --- a/FDTD/extensions/engine_ext_lorentzmaterial.cpp +++ b/FDTD/extensions/engine_ext_lorentzmaterial.cpp @@ -23,11 +23,57 @@ Engine_Ext_LorentzMaterial::Engine_Ext_LorentzMaterial(Operator_Ext_LorentzMater { m_Op_Ext_Lor = op_ext_lorentz; m_Order = m_Op_Ext_Lor->GetDispersionOrder(); + int order = m_Op_Ext_Lor->m_Order; + + curr_Lor_ADE = new FDTD_FLOAT**[order]; + volt_Lor_ADE = new FDTD_FLOAT**[order]; + for (int o=0;om_curr_Lor_ADE_On[o]==true) + { + curr_Lor_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Lor->m_LM_Count[o]]; + for (unsigned int i=0; im_LM_Count[o]; ++i) + curr_Lor_ADE[o][n][i]=0.0; + } + else + curr_Lor_ADE[o][n] = NULL; + + if (m_Op_Ext_Lor->m_volt_Lor_ADE_On[o]==true) + { + volt_Lor_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Lor->m_LM_Count[o]]; + for (unsigned int i=0; im_LM_Count[o]; ++i) + volt_Lor_ADE[o][n][i]=0.0; + } + else + volt_Lor_ADE[o][n] = NULL; + } + } } Engine_Ext_LorentzMaterial::~Engine_Ext_LorentzMaterial() { + if (curr_Lor_ADE==NULL && volt_Lor_ADE==NULL) + return; + for (int o=0;om_Order;++o) + { + for (int n=0; n<3; ++n) + { + delete[] curr_Lor_ADE[o][n]; + delete[] volt_Lor_ADE[o][n]; + } + delete[] curr_Lor_ADE[o]; + delete[] volt_Lor_ADE[o]; + } + delete[] curr_Lor_ADE; + curr_Lor_ADE=NULL; + + delete[] volt_Lor_ADE; + volt_Lor_ADE=NULL; } void Engine_Ext_LorentzMaterial::DoPreVoltageUpdates() @@ -38,53 +84,116 @@ void Engine_Ext_LorentzMaterial::DoPreVoltageUpdates() unsigned int **pos = m_Op_Ext_Lor->m_LM_pos[o]; - //switch for different engine types to access faster inline engine functions - switch (m_Eng->GetType()) + if (m_Op_Ext_Lor->m_volt_Lor_ADE_On[o]) { - case Engine::BASIC: - { - for (unsigned int i=0; im_LM_Count.at(o); ++i) + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) { - volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; - volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); + case Engine::BASIC: + { + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + volt_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][0][i]*volt_ADE[o][0][i]; + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * (m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][0][i]); - volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; - volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); + volt_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][1][i]*volt_ADE[o][1][i]; + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * (m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); - volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; - volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + volt_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][2][i]*volt_ADE[o][2][i]; + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * (m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + volt_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][0][i]*volt_ADE[o][0][i]; + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * (eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][0][i]); + + volt_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][1][i]*volt_ADE[o][1][i]; + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * (eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][1][i]); + + volt_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][2][i]*volt_ADE[o][2][i]; + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * (eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); + } + break; + } + default: + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + volt_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][0][i]*volt_ADE[o][0][i]; + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * (m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][0][i]); + + volt_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][1][i]*volt_ADE[o][1][i]; + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * (m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][1][i]); + + volt_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->v_Lor_ADE[o][2][i]*volt_ADE[o][2][i]; + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * (m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i])-volt_Lor_ADE[o][2][i]); + } + break; } - break; } - case Engine::SSE: + else { - Engine_sse* eng_sse = (Engine_sse*)m_Eng; - for (unsigned int i=0; im_LM_Count.at(o); ++i) + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) { - volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; - volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); - - volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; - volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); - - volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; - volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); - } - break; - } - default: - for (unsigned int i=0; im_LM_Count.at(o); ++i) + case Engine::BASIC: { - volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; - volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); - volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; - volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); - volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; - volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + default: + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i]; + volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i]; + volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]); + + volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i]; + volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; } - break; } } } @@ -97,53 +206,116 @@ void Engine_Ext_LorentzMaterial::DoPreCurrentUpdates() unsigned int **pos = m_Op_Ext_Lor->m_LM_pos[o]; - //switch for different engine types to access faster inline engine functions - switch (m_Eng->GetType()) + if (m_Op_Ext_Lor->m_curr_Lor_ADE_On[o]) { - case Engine::BASIC: - { - for (unsigned int i=0; im_LM_Count.at(o); ++i) + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) { - curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; - curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); + case Engine::BASIC: + { + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + curr_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][0][i]*curr_ADE[o][0][i]; + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * (m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][0][i]); - curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; - curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); + curr_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][1][i]*curr_ADE[o][1][i]; + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * (m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][1][i]); - curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; - curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + curr_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][2][i]*curr_ADE[o][2][i]; + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * (m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + curr_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][0][i]*curr_ADE[o][0][i]; + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * (eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][0][i]); + + curr_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][1][i]*curr_ADE[o][1][i]; + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * (eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][1][i]); + + curr_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][2][i]*curr_ADE[o][2][i]; + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * (eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][2][i]); + } + break; + } + default: + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + curr_Lor_ADE[o][0][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][0][i]*curr_ADE[o][0][i]; + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * (m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][0][i]); + + curr_Lor_ADE[o][1][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][1][i]*curr_ADE[o][1][i]; + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * (m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][1][i]); + + curr_Lor_ADE[o][2][i]+=m_Op_Ext_Lor->i_Lor_ADE[o][2][i]*curr_ADE[o][2][i]; + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * (m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i])-curr_Lor_ADE[o][2][i]); + } + break; } - break; } - case Engine::SSE: + else { - Engine_sse* eng_sse = (Engine_sse*)m_Eng; - for (unsigned int i=0; im_LM_Count.at(o); ++i) + //switch for different engine types to access faster inline engine functions + switch (m_Eng->GetType()) { - curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; - curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); - - curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; - curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); - - curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; - curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); - } - break; - } - default: - for (unsigned int i=0; im_LM_Count.at(o); ++i) + case Engine::BASIC: { - curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; - curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); - curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; - curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); - curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; - curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + case Engine::SSE: + { + Engine_sse* eng_sse = (Engine_sse*)m_Eng; + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; + } + default: + for (unsigned int i=0; im_LM_Count.at(o); ++i) + { + curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i]; + curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i]; + curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]); + + curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i]; + curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]); + } + break; } - break; } } } diff --git a/FDTD/extensions/engine_ext_lorentzmaterial.h b/FDTD/extensions/engine_ext_lorentzmaterial.h index 3ac6850..804f49a 100644 --- a/FDTD/extensions/engine_ext_lorentzmaterial.h +++ b/FDTD/extensions/engine_ext_lorentzmaterial.h @@ -34,6 +34,15 @@ public: protected: Operator_Ext_LorentzMaterial* m_Op_Ext_Lor; + + //! ADE Lorentz voltages + // Array setup: volt_Lor_ADE[N_order][direction][mesh_pos] + FDTD_FLOAT ***volt_Lor_ADE; + + //! ADE Lorentz currents + // Array setup: curr_Lor_ADE[N_order][direction][mesh_pos] + FDTD_FLOAT ***curr_Lor_ADE; + }; #endif // ENGINE_EXT_LORENTZMATERIAL_H diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/FDTD/extensions/operator_ext_lorentzmaterial.cpp index 723bc36..3b303b8 100644 --- a/FDTD/extensions/operator_ext_lorentzmaterial.cpp +++ b/FDTD/extensions/operator_ext_lorentzmaterial.cpp @@ -28,6 +28,12 @@ Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Opera v_ext_ADE = NULL; i_int_ADE = NULL; i_ext_ADE = NULL; + + v_Lor_ADE = NULL; + i_Lor_ADE = NULL; + + m_curr_Lor_ADE_On = NULL; + m_curr_Lor_ADE_On = NULL; } Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op, Operator_Ext_LorentzMaterial* op_ext) : Operator_Ext_Dispersive(op,op_ext) @@ -36,6 +42,12 @@ Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op, Operato v_ext_ADE = NULL; i_int_ADE = NULL; i_ext_ADE = NULL; + + v_Lor_ADE = NULL; + i_Lor_ADE = NULL; + + m_curr_Lor_ADE_On = NULL; + m_curr_Lor_ADE_On = NULL; } Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial() @@ -54,6 +66,10 @@ Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial() delete[] i_int_ADE[i][n]; delete[] i_ext_ADE[i][n]; } + if (m_volt_Lor_ADE_On[i]) + delete[] v_Lor_ADE[i][n]; + if (m_curr_Lor_ADE_On[i]) + delete[] i_Lor_ADE[i][n]; } if (m_volt_ADE_On[i]) { @@ -65,6 +81,10 @@ Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial() delete[] i_int_ADE[i]; delete[] i_ext_ADE[i]; } + if (m_volt_Lor_ADE_On[i]) + delete[] v_Lor_ADE[i]; + if (m_curr_Lor_ADE_On[i]) + delete[] i_Lor_ADE[i]; } delete[] v_int_ADE; delete[] v_ext_ADE; @@ -74,6 +94,16 @@ Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial() v_ext_ADE = NULL; i_int_ADE = NULL; i_ext_ADE = NULL; + + delete[] v_Lor_ADE; + delete[] i_Lor_ADE; + v_Lor_ADE = NULL; + i_Lor_ADE = NULL; + + delete[] m_curr_Lor_ADE_On; + delete[] m_volt_Lor_ADE_On; + m_curr_Lor_ADE_On = NULL; + m_curr_Lor_ADE_On = NULL; } Operator_Extension* Operator_Ext_LorentzMaterial::Clone(Operator* op) @@ -91,15 +121,24 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() unsigned int numLines[3] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)}; CSPropLorentzMaterial* mat = NULL; - double w_plasma,t_relax; bool b_pos_on; + vector v_pos[3]; + + // drude material parameter + double w_plasma,t_relax; double L_D[3], C_D[3]; double R_D[3], G_D[3]; vector v_int[3]; vector v_ext[3]; vector i_int[3]; vector i_ext[3]; - vector v_pos[3]; + + //additional Dorentz material parameter + double w_Lor_Pol; + double C_L[3]; + double L_L[3]; + vector v_Lor[3]; + vector i_Lor[3]; m_Order = 0; vector LD_props = m_Op->CSX->GetPropertyByType(CSProperties::LORENTZMATERIAL); @@ -112,27 +151,40 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() m_Order=LorMat->GetDispersionOrder(); } + m_LM_pos = new unsigned int**[m_Order]; + m_volt_ADE_On = new bool[m_Order]; m_curr_ADE_On = new bool[m_Order]; - m_LM_pos = new unsigned int**[m_Order]; + m_volt_Lor_ADE_On = new bool[m_Order]; + m_curr_Lor_ADE_On = new bool[m_Order]; v_int_ADE = new FDTD_FLOAT**[m_Order]; v_ext_ADE = new FDTD_FLOAT**[m_Order]; i_int_ADE = new FDTD_FLOAT**[m_Order]; i_ext_ADE = new FDTD_FLOAT**[m_Order]; + v_Lor_ADE = new FDTD_FLOAT**[m_Order]; + i_Lor_ADE = new FDTD_FLOAT**[m_Order]; + for (int order=0;orderGetYeeCoords(n,pos,coord,false)==false) continue; if (m_CC_R0_included && (n==2) && (pos[0]==0)) coord[1] = m_Op->GetDiscLine(1,0); + if (m_Op->GetVI(n,pos[0],pos[1],pos[2])==0) + continue; + CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); if ((mat = prop->ToLorentzMaterial())) { @@ -168,6 +224,12 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() { R_D[n] = L_D[n]/t_relax; } + w_Lor_Pol = mat->GetEpsLorPoleFreqWeighted(order,n,coord) * 2 * PI; + if ((w_Lor_Pol>0) && (L_D[n]>0)) + { + m_volt_Lor_ADE_On[order] = true; + C_L[n] = 1/(w_Lor_Pol*w_Lor_Pol*L_D[n]); + } } } @@ -175,8 +237,11 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() { C_D[n]=0; G_D[n]=0; + L_L[n]=0; if (m_Op->GetYeeCoords(n,pos,coord,true)==false) continue; + if (m_Op->GetIV(n,pos[0],pos[1],pos[2])==0) + continue; CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true); if ((mat = prop->ToLorentzMaterial())) @@ -193,10 +258,16 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() { G_D[n] = C_D[n]/t_relax; } + w_Lor_Pol = mat->GetMueLorPoleFreqWeighted(order,n,coord) * 2 * PI; + if ((w_Lor_Pol>0) && (C_D[n]>0)) + { + m_curr_Lor_ADE_On[order] = true; + L_L[n] = 1/(w_Lor_Pol*w_Lor_Pol*C_D[n]); + } } } - if (b_pos_on) //this position has active lorentz material + if (b_pos_on) //this position has active drude material { for (unsigned int n=0; n<3; ++n) { @@ -225,7 +296,14 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() i_int[n].push_back(1); i_ext[n].push_back(0); } - // cerr << v_int[n].back() << " " << v_ext[n].back() << " " << i_int[n].back() << " " << i_ext[n].back() << endl; + if (C_L[n]>0) + v_Lor[n].push_back(dT/C_L[n]/m_Op->GetVI(n,pos[0],pos[1],pos[2])); + else + v_Lor[n].push_back(0); + if (L_L[n]>0) + i_Lor[n].push_back(dT/L_L[n]/m_Op->GetIV(n,pos[0],pos[1],pos[2])); + else + i_Lor[n].push_back(0); } } } @@ -259,6 +337,16 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() i_ext_ADE[order] = NULL; } + if (m_volt_Lor_ADE_On[order]) + v_Lor_ADE[order] = new FDTD_FLOAT*[3]; + else + v_Lor_ADE[order] = NULL; + + if (m_curr_Lor_ADE_On[order]) + i_Lor_ADE[order] = new FDTD_FLOAT*[3]; + else + i_Lor_ADE[order] = NULL; + for (int n=0; n<3; ++n) { m_LM_pos[order][n] = new unsigned int[m_LM_Count.at(order)]; @@ -286,6 +374,19 @@ bool Operator_Ext_LorentzMaterial::BuildExtension() i_ext_ADE[order][n][i] = i_ext[n].at(i); } } + + if (m_volt_Lor_ADE_On[order]) + { + v_Lor_ADE[order][n] = new FDTD_FLOAT[m_LM_Count.at(order)]; + for (unsigned int i=0; i