diff --git a/FDTD/excitation.cpp b/FDTD/excitation.cpp index 893d832..e11cda6 100644 --- a/FDTD/excitation.cpp +++ b/FDTD/excitation.cpp @@ -27,25 +27,9 @@ Excitation::Excitation( double timestep ) { Signal_volt = 0; Signal_curr = 0; - Volt_delay = 0; - Volt_amp = 0; - Volt_dir = 0; - Volt_Count = 0; - Curr_delay = 0; - Curr_amp = 0; - Curr_dir = 0; - Curr_Count = 0; m_Excit_Type = -1; - for (int n=0; n<3; ++n) - { - Volt_index[n] = 0; - Curr_index[n] = 0; - Volt_Count_Dir[n] = 0; - Curr_Count_Dir[n] = 0; - } - dT = timestep; m_nyquistTS = 0; } @@ -54,17 +38,6 @@ Excitation::~Excitation() { delete[] Signal_volt; delete[] Signal_curr; - delete[] Volt_delay; - delete[] Volt_dir; - delete[] Volt_amp; - delete[] Curr_delay; - delete[] Curr_dir; - delete[] Curr_amp; - for (int n=0; n<3; ++n) - { - delete[] Volt_index[n]; - delete[] Curr_index[n]; - } } void Excitation::Reset( double timestep ) @@ -73,35 +46,9 @@ void Excitation::Reset( double timestep ) Signal_volt = 0; delete[] Signal_curr; Signal_curr = 0; - delete[] Volt_delay; - Volt_delay = 0; - delete[] Volt_dir; - Volt_dir = 0; - delete[] Volt_amp; - Volt_amp = 0; - delete[] Curr_delay; - Curr_delay = 0; - delete[] Curr_dir; - Curr_dir = 0; - delete[] Curr_amp; - Curr_amp = 0; - - Volt_Count = 0; - Curr_Count = 0; m_Excit_Type = -1; - for (int n=0; n<3; ++n) - { - delete[] Volt_index[n]; - Volt_index[n] = 0; - delete[] Curr_index[n]; - Curr_index[n] = 0; - - Volt_Count_Dir[n] = 0; - Curr_Count_Dir[n] = 0; - } - dT = timestep; m_nyquistTS = 0; } @@ -311,67 +258,3 @@ void Excitation::DumpCurrentExcite(string filename) file.close(); } -void Excitation::setupVoltageExcitation( vector const volt_vIndex[3], vector const& volt_vExcit, - vector const& volt_vDelay, vector const& volt_vDir ) -{ - Volt_Count = volt_vIndex[0].size(); - for (int n=0; n<3; n++) - { - Volt_Count_Dir[n]=0; - delete[] Volt_index[n]; - Volt_index[n] = new unsigned int[Volt_Count]; - } - delete[] Volt_delay; - delete[] Volt_amp; - delete[] Volt_dir; - Volt_delay = new unsigned int[Volt_Count]; - Volt_amp = new FDTD_FLOAT[Volt_Count]; - Volt_dir = new unsigned short[Volt_Count]; - -// cerr << "Excitation::setupVoltageExcitation(): Number of voltage excitation points: " << Volt_Count << endl; -// if (Volt_Count==0) -// cerr << "No E-Field/voltage excitation found!" << endl; - for (int n=0; n<3; n++) - for (unsigned int i=0; i const curr_vIndex[3], vector const& curr_vExcit, - vector const& curr_vDelay, vector const& curr_vDir ) -{ - Curr_Count = curr_vIndex[0].size(); - for (int n=0; n<3; n++) - { - Curr_Count_Dir[n]=0; - delete[] Curr_index[n]; - Curr_index[n] = new unsigned int[Curr_Count]; - } - delete[] Curr_delay; - delete[] Curr_amp; - delete[] Curr_dir; - Curr_delay = new unsigned int[Curr_Count]; - Curr_amp = new FDTD_FLOAT[Curr_Count]; - Curr_dir = new unsigned short[Curr_Count]; - -// cerr << "Excitation::setupCurrentExcitation(): Number of current excitation points: " << Curr_Count << endl; -// if (Curr_Count==0) -// cerr << "No H-Field/current excitation found!" << endl; - for (int n=0; n<3; ++n) - for (unsigned int i=0; i const volt_vIndex[3], vector const& volt_vExcit, - vector const& volt_vDelay, vector const& volt_vDir ); - void setupCurrentExcitation( vector const curr_vIndex[3], vector const& curr_vExcit, - vector const& curr_vDelay, vector const& curr_vDir ); - void SetNyquistNum(unsigned int nyquist) {m_nyquistTS=nyquist;} unsigned int GetNyquistNum() const {return m_nyquistTS;} @@ -57,32 +52,22 @@ public: //! Get the type of excitation int GetExciteType() const {return m_Excit_Type;} - //Excitation time-signal - unsigned int Length; - FDTD_FLOAT* Signal_volt; - FDTD_FLOAT* Signal_curr; + //! Get the length of the excitation signal + unsigned int GetLength() const {return Length;} - //E-Field/voltage Excitation - unsigned int Volt_Count; - unsigned int Volt_Count_Dir[3]; - unsigned int* Volt_index[3]; - unsigned short* Volt_dir; - FDTD_FLOAT* Volt_amp; //represented as edge-voltages!! - unsigned int* Volt_delay; - - //H-Field/current Excitation - unsigned int Curr_Count; - unsigned int Curr_Count_Dir[3]; - unsigned int* Curr_index[3]; - unsigned short* Curr_dir; - FDTD_FLOAT* Curr_amp; //represented as edge-currents!! - unsigned int* Curr_delay; + FDTD_FLOAT* GetVoltageSignal() const {return Signal_volt;} + FDTD_FLOAT* GetCurrentSignal() const {return Signal_curr;} protected: double dT; unsigned int m_nyquistTS; int m_Excit_Type; + //Excitation time-signal + unsigned int Length; + FDTD_FLOAT* Signal_volt; + FDTD_FLOAT* Signal_curr; + //! Calculate a custom signal virtual void CalcCustomExcitation(double f0, int nTS, string signal); //! Calculate an excitation with center of \a f0 and the half bandwidth \a fc diff --git a/FDTD/extensions/engine_ext_excitation.cpp b/FDTD/extensions/engine_ext_excitation.cpp index ec328d2..bcb7b6e 100644 --- a/FDTD/extensions/engine_ext_excitation.cpp +++ b/FDTD/extensions/engine_ext_excitation.cpp @@ -37,50 +37,52 @@ void Engine_Ext_Excitation::Apply2Voltages() unsigned int ny; unsigned int pos[3]; int numTS = m_Eng->GetNumberOfTimesteps(); + unsigned int length = m_Op_Exc->m_Exc->GetLength(); + FDTD_FLOAT* exc_volt = m_Op_Exc->m_Exc->GetVoltageSignal(); //switch for different engine types to access faster inline engine functions switch (m_Eng->GetType()) { case Engine::BASIC: { - for (unsigned int n=0; nm_Exc->Volt_Count; ++n) + for (unsigned int n=0; nVolt_Count; ++n) { - exc_pos = numTS - (int)m_Op_Exc->m_Exc->Volt_delay[n]; - exc_pos *= (exc_pos>0 && exc_pos<=(int)m_Op_Exc->m_Exc->Length); - ny = m_Op_Exc->m_Exc->Volt_dir[n]; - pos[0]=m_Op_Exc->m_Exc->Volt_index[0][n]; - pos[1]=m_Op_Exc->m_Exc->Volt_index[1][n]; - pos[2]=m_Op_Exc->m_Exc->Volt_index[2][n]; - m_Eng->Engine::SetVolt(ny,pos, m_Eng->Engine::GetVolt(ny,pos) + m_Op_Exc->m_Exc->Volt_amp[n]*m_Op_Exc->m_Exc->Signal_volt[exc_pos]); + exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; + exc_pos *= (exc_pos>0 && exc_pos<=(int)length); + ny = m_Op_Exc->Volt_dir[n]; + pos[0]=m_Op_Exc->Volt_index[0][n]; + pos[1]=m_Op_Exc->Volt_index[1][n]; + pos[2]=m_Op_Exc->Volt_index[2][n]; + m_Eng->Engine::SetVolt(ny,pos, m_Eng->Engine::GetVolt(ny,pos) + m_Op_Exc->Volt_amp[n]*exc_volt[exc_pos]); } break; } case Engine::SSE: { - for (unsigned int n=0; nm_Exc->Volt_Count; ++n) + for (unsigned int n=0; nVolt_Count; ++n) { Engine_sse* eng_sse = (Engine_sse*) m_Eng; - exc_pos = numTS - (int)m_Op_Exc->m_Exc->Volt_delay[n]; - exc_pos *= (exc_pos>0 && exc_pos<=(int)m_Op_Exc->m_Exc->Length); - ny = m_Op_Exc->m_Exc->Volt_dir[n]; - pos[0]=m_Op_Exc->m_Exc->Volt_index[0][n]; - pos[1]=m_Op_Exc->m_Exc->Volt_index[1][n]; - pos[2]=m_Op_Exc->m_Exc->Volt_index[2][n]; - eng_sse->Engine_sse::SetVolt(ny,pos, eng_sse->Engine_sse::GetVolt(ny,pos) + m_Op_Exc->m_Exc->Volt_amp[n]*m_Op_Exc->m_Exc->Signal_volt[exc_pos]); + exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; + exc_pos *= (exc_pos>0 && exc_pos<=(int)length); + ny = m_Op_Exc->Volt_dir[n]; + pos[0]=m_Op_Exc->Volt_index[0][n]; + pos[1]=m_Op_Exc->Volt_index[1][n]; + pos[2]=m_Op_Exc->Volt_index[2][n]; + eng_sse->Engine_sse::SetVolt(ny,pos, eng_sse->Engine_sse::GetVolt(ny,pos) + m_Op_Exc->Volt_amp[n]*exc_volt[exc_pos]); } break; } default: { - for (unsigned int n=0; nm_Exc->Volt_Count; ++n) + for (unsigned int n=0; nVolt_Count; ++n) { - exc_pos = numTS - (int)m_Op_Exc->m_Exc->Volt_delay[n]; - exc_pos *= (exc_pos>0 && exc_pos<=(int)m_Op_Exc->m_Exc->Length); - ny = m_Op_Exc->m_Exc->Volt_dir[n]; - pos[0]=m_Op_Exc->m_Exc->Volt_index[0][n]; - pos[1]=m_Op_Exc->m_Exc->Volt_index[1][n]; - pos[2]=m_Op_Exc->m_Exc->Volt_index[2][n]; - m_Eng->SetVolt(ny,pos, m_Eng->GetVolt(ny,pos) + m_Op_Exc->m_Exc->Volt_amp[n]*m_Op_Exc->m_Exc->Signal_volt[exc_pos]); + exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n]; + exc_pos *= (exc_pos>0 && exc_pos<=(int)length); + ny = m_Op_Exc->Volt_dir[n]; + pos[0]=m_Op_Exc->Volt_index[0][n]; + pos[1]=m_Op_Exc->Volt_index[1][n]; + pos[2]=m_Op_Exc->Volt_index[2][n]; + m_Eng->SetVolt(ny,pos, m_Eng->GetVolt(ny,pos) + m_Op_Exc->Volt_amp[n]*exc_volt[exc_pos]); } break; } @@ -95,50 +97,52 @@ void Engine_Ext_Excitation::Apply2Current() unsigned int ny; unsigned int pos[3]; int numTS = m_Eng->GetNumberOfTimesteps(); + unsigned int length = m_Op_Exc->m_Exc->GetLength(); + FDTD_FLOAT* exc_curr = m_Op_Exc->m_Exc->GetCurrentSignal(); //switch for different engine types to access faster inline engine functions switch (m_Eng->GetType()) { case Engine::BASIC: { - for (unsigned int n=0; nm_Exc->Curr_Count; ++n) + for (unsigned int n=0; nCurr_Count; ++n) { - exc_pos = numTS - (int)m_Op_Exc->m_Exc->Curr_delay[n]; - exc_pos *= (exc_pos>0 && exc_pos<=(int)m_Op_Exc->m_Exc->Length); - ny = m_Op_Exc->m_Exc->Curr_dir[n]; - pos[0]=m_Op_Exc->m_Exc->Curr_index[0][n]; - pos[1]=m_Op_Exc->m_Exc->Curr_index[1][n]; - pos[2]=m_Op_Exc->m_Exc->Curr_index[2][n]; - m_Eng->Engine::SetCurr(ny,pos, m_Eng->Engine::GetCurr(ny,pos) + m_Op_Exc->m_Exc->Curr_amp[n]*m_Op_Exc->m_Exc->Signal_curr[exc_pos]); + exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; + exc_pos *= (exc_pos>0 && exc_pos<=(int)length); + ny = m_Op_Exc->Curr_dir[n]; + pos[0]=m_Op_Exc->Curr_index[0][n]; + pos[1]=m_Op_Exc->Curr_index[1][n]; + pos[2]=m_Op_Exc->Curr_index[2][n]; + m_Eng->Engine::SetCurr(ny,pos, m_Eng->Engine::GetCurr(ny,pos) + m_Op_Exc->Curr_amp[n]*exc_curr[exc_pos]); } break; } case Engine::SSE: { - for (unsigned int n=0; nm_Exc->Curr_Count; ++n) + for (unsigned int n=0; nCurr_Count; ++n) { Engine_sse* eng_sse = (Engine_sse*) m_Eng; - exc_pos = numTS - (int)m_Op_Exc->m_Exc->Curr_delay[n]; - exc_pos *= (exc_pos>0 && exc_pos<=(int)m_Op_Exc->m_Exc->Length); - ny = m_Op_Exc->m_Exc->Curr_dir[n]; - pos[0]=m_Op_Exc->m_Exc->Curr_index[0][n]; - pos[1]=m_Op_Exc->m_Exc->Curr_index[1][n]; - pos[2]=m_Op_Exc->m_Exc->Curr_index[2][n]; - eng_sse->Engine_sse::SetCurr(ny,pos, eng_sse->Engine_sse::GetCurr(ny,pos) + m_Op_Exc->m_Exc->Curr_amp[n]*m_Op_Exc->m_Exc->Signal_curr[exc_pos]); + exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; + exc_pos *= (exc_pos>0 && exc_pos<=(int)length); + ny = m_Op_Exc->Curr_dir[n]; + pos[0]=m_Op_Exc->Curr_index[0][n]; + pos[1]=m_Op_Exc->Curr_index[1][n]; + pos[2]=m_Op_Exc->Curr_index[2][n]; + eng_sse->Engine_sse::SetCurr(ny,pos, eng_sse->Engine_sse::GetCurr(ny,pos) + m_Op_Exc->Curr_amp[n]*exc_curr[exc_pos]); } break; } default: { - for (unsigned int n=0; nm_Exc->Curr_Count; ++n) + for (unsigned int n=0; nCurr_Count; ++n) { - exc_pos = numTS - (int)m_Op_Exc->m_Exc->Curr_delay[n]; - exc_pos *= (exc_pos>0 && exc_pos<=(int)m_Op_Exc->m_Exc->Length); - ny = m_Op_Exc->m_Exc->Curr_dir[n]; - pos[0]=m_Op_Exc->m_Exc->Curr_index[0][n]; - pos[1]=m_Op_Exc->m_Exc->Curr_index[1][n]; - pos[2]=m_Op_Exc->m_Exc->Curr_index[2][n]; - m_Eng->SetCurr(ny,pos, m_Eng->GetCurr(ny,pos) + m_Op_Exc->m_Exc->Curr_amp[n]*m_Op_Exc->m_Exc->Signal_curr[exc_pos]); + exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n]; + exc_pos *= (exc_pos>0 && exc_pos<=(int)length); + ny = m_Op_Exc->Curr_dir[n]; + pos[0]=m_Op_Exc->Curr_index[0][n]; + pos[1]=m_Op_Exc->Curr_index[1][n]; + pos[2]=m_Op_Exc->Curr_index[2][n]; + m_Eng->SetCurr(ny,pos, m_Eng->GetCurr(ny,pos) + m_Op_Exc->Curr_amp[n]*exc_curr[exc_pos]); } break; } diff --git a/FDTD/extensions/engine_ext_mur_abc.cpp b/FDTD/extensions/engine_ext_mur_abc.cpp index d877497..9e826be 100644 --- a/FDTD/extensions/engine_ext_mur_abc.cpp +++ b/FDTD/extensions/engine_ext_mur_abc.cpp @@ -21,6 +21,7 @@ #include "FDTD/engine_sse.h" #include "tools/array_ops.h" #include "tools/useful.h" +#include "operator_ext_excitation.h" Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext) : Engine_Extension(op_ext) { @@ -41,18 +42,20 @@ Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC(Operator_Ext_Mur_ABC* op_ext) : Engine_Ex //find if some excitation is on this mur-abc and find the max length of this excite, so that the abc can start after the excitation is done... int maxDelay=-1; - for (unsigned int n=0; nm_Op->Exc->Volt_Count; ++n) + Excitation* Exc = m_Op_mur->m_Op->GetExcitationSignal(); + Operator_Ext_Excitation* Exc_ext = m_Op_mur->m_Op->GetExcitationExtension(); + for (unsigned int n=0; nGetVoltCount(); ++n) { - if ( ((m_Op_mur->m_Op->Exc->Volt_dir[n]==m_nyP) || (m_Op_mur->m_Op->Exc->Volt_dir[n]==m_nyPP)) && (m_Op_mur->m_Op->Exc->Volt_index[m_ny][n]==m_LineNr) ) + if ( ((Exc_ext->Volt_dir[n]==m_nyP) || (Exc_ext->Volt_dir[n]==m_nyPP)) && (Exc_ext->Volt_index[m_ny][n]==m_LineNr) ) { - if ((int)m_Op_mur->m_Op->Exc->Volt_delay[n]>maxDelay) - maxDelay = (int)m_Op_mur->m_Op->Exc->Volt_delay[n]; + if ((int)Exc_ext->Volt_delay[n]>maxDelay) + maxDelay = (int)Exc_ext->Volt_delay[n]; } } m_start_TS = 0; if (maxDelay>=0) { - m_start_TS = maxDelay + m_Op_mur->m_Op->Exc->Length + 10; //give it some extra timesteps, for the excitation to travel at least one cell away + m_start_TS = maxDelay + m_Op_mur->m_Op->GetExcitationSignal()->GetLength() + 10; //give it some extra timesteps, for the excitation to travel at least one cell away cerr << "Engine_Ext_Mur_ABC::Engine_Ext_Mur_ABC: Warning: Excitation inside the Mur-ABC #" << m_ny << "-" << (int)(m_LineNr>0) << " found!!!! Mur-ABC will be switched on after excitation is done at " << m_start_TS << " timesteps!!! " << endl; } diff --git a/FDTD/extensions/operator_ext_excitation.cpp b/FDTD/extensions/operator_ext_excitation.cpp index 3e9e641..be1a0ff 100644 --- a/FDTD/extensions/operator_ext_excitation.cpp +++ b/FDTD/extensions/operator_ext_excitation.cpp @@ -20,15 +20,64 @@ #include "FDTD/excitation.h" #include "ContinuousStructure.h" -Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op, Excitation* exc) : Operator_Extension(op) +Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op) : Operator_Extension(op) { - m_Exc = exc; + Volt_delay = 0; + Volt_amp = 0; + Volt_dir = 0; + Volt_Count = 0; + Curr_delay = 0; + Curr_amp = 0; + Curr_dir = 0; + Curr_Count = 0; + + for (int n=0; n<3; ++n) + { + Volt_index[n] = 0; + Curr_index[n] = 0; + Volt_Count_Dir[n] = 0; + Curr_Count_Dir[n] = 0; + } + + m_Exc = m_Op->m_Exc; } Operator_Ext_Excitation::~Operator_Ext_Excitation() { + Reset(); } +void Operator_Ext_Excitation::Reset( ) +{ + delete[] Volt_delay; + Volt_delay = 0; + delete[] Volt_dir; + Volt_dir = 0; + delete[] Volt_amp; + Volt_amp = 0; + delete[] Curr_delay; + Curr_delay = 0; + delete[] Curr_dir; + Curr_dir = 0; + delete[] Curr_amp; + Curr_amp = 0; + + Volt_Count = 0; + Curr_Count = 0; + + for (int n=0; n<3; ++n) + { + delete[] Volt_index[n]; + Volt_index[n] = 0; + delete[] Curr_index[n]; + Curr_index[n] = 0; + + Volt_Count_Dir[n] = 0; + Curr_Count_Dir[n] = 0; + } +} + + Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op, Operator_Ext_Excitation* op_ext) : Operator_Extension(op, op_ext) { m_Exc = NULL; @@ -36,7 +85,7 @@ Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op, Operator_Ext_Exci bool Operator_Ext_Excitation::BuildExtension() { - double dT = m_Exc->GetTimestep(); + double dT = m_Op->GetTimestep(); if (dT==0) return false; if (m_Exc==0) @@ -217,14 +266,79 @@ bool Operator_Ext_Excitation::BuildExtension() } // set voltage excitations - m_Exc->setupVoltageExcitation( volt_vIndex, volt_vExcit, volt_vDelay, volt_vDir ); + setupVoltageExcitation( volt_vIndex, volt_vExcit, volt_vDelay, volt_vDir ); // set current excitations - m_Exc->setupCurrentExcitation( curr_vIndex, curr_vExcit, curr_vDelay, curr_vDir ); + setupCurrentExcitation( curr_vIndex, curr_vExcit, curr_vDelay, curr_vDir ); return true; } +void Operator_Ext_Excitation::setupVoltageExcitation( vector const volt_vIndex[3], vector const& volt_vExcit, + vector const& volt_vDelay, vector const& volt_vDir ) +{ + Volt_Count = volt_vIndex[0].size(); + for (int n=0; n<3; n++) + { + Volt_Count_Dir[n]=0; + delete[] Volt_index[n]; + Volt_index[n] = new unsigned int[Volt_Count]; + } + delete[] Volt_delay; + delete[] Volt_amp; + delete[] Volt_dir; + Volt_delay = new unsigned int[Volt_Count]; + Volt_amp = new FDTD_FLOAT[Volt_Count]; + Volt_dir = new unsigned short[Volt_Count]; + +// cerr << "Excitation::setupVoltageExcitation(): Number of voltage excitation points: " << Volt_Count << endl; +// if (Volt_Count==0) +// cerr << "No E-Field/voltage excitation found!" << endl; + for (int n=0; n<3; n++) + for (unsigned int i=0; i const curr_vIndex[3], vector const& curr_vExcit, + vector const& curr_vDelay, vector const& curr_vDir ) +{ + Curr_Count = curr_vIndex[0].size(); + for (int n=0; n<3; n++) + { + Curr_Count_Dir[n]=0; + delete[] Curr_index[n]; + Curr_index[n] = new unsigned int[Curr_Count]; + } + delete[] Curr_delay; + delete[] Curr_amp; + delete[] Curr_dir; + Curr_delay = new unsigned int[Curr_Count]; + Curr_amp = new FDTD_FLOAT[Curr_Count]; + Curr_dir = new unsigned short[Curr_Count]; + +// cerr << "Excitation::setupCurrentExcitation(): Number of current excitation points: " << Curr_Count << endl; +// if (Curr_Count==0) +// cerr << "No H-Field/current excitation found!" << endl; + for (int n=0; n<3; ++n) + for (unsigned int i=0; iVolt_Count << "\t (" << m_Exc->Volt_Count_Dir[0] << ", " << m_Exc->Volt_Count_Dir[1] << ", " << m_Exc->Volt_Count_Dir[2] << ")" << endl; - cout << "Current excitations\t: " << m_Exc->Curr_Count << "\t (" << m_Exc->Curr_Count_Dir[0] << ", " << m_Exc->Curr_Count_Dir[1] << ", " << m_Exc->Curr_Count_Dir[2] << ")" << endl; - cout << "Excitation Length (TS)\t: " << m_Exc->Length << endl; - cout << "Excitation Length (s)\t: " << m_Exc->Length*m_Op->GetTimestep() << endl; + cout << "Voltage excitations\t: " << Volt_Count << "\t (" << Volt_Count_Dir[0] << ", " << Volt_Count_Dir[1] << ", " << Volt_Count_Dir[2] << ")" << endl; + cout << "Current excitations\t: " << Curr_Count << "\t (" << Curr_Count_Dir[0] << ", " << Curr_Count_Dir[1] << ", " << Curr_Count_Dir[2] << ")" << endl; + cout << "Excitation Length (TS)\t: " << m_Exc->GetLength() << endl; + cout << "Excitation Length (s)\t: " << m_Exc->GetLength()*m_Op->GetTimestep() << endl; } diff --git a/FDTD/extensions/operator_ext_excitation.h b/FDTD/extensions/operator_ext_excitation.h index 721def1..32b5e87 100644 --- a/FDTD/extensions/operator_ext_excitation.h +++ b/FDTD/extensions/operator_ext_excitation.h @@ -26,8 +26,10 @@ class Excitation; class Operator_Ext_Excitation : public Operator_Extension { friend class Engine_Ext_Excitation; + friend class Engine_Ext_Mur_ABC; + friend class Operator; public: - Operator_Ext_Excitation(Operator* op, Excitation* exc); + Operator_Ext_Excitation(Operator* op); ~Operator_Ext_Excitation(); virtual Operator_Extension* Clone(Operator* op) {UNUSED(op);return NULL;} @@ -43,10 +45,38 @@ public: virtual void ShowStat(ostream &ostr) const; + virtual void Reset(); + + unsigned int GetVoltCount() const {return Volt_Count;} + unsigned int GetVoltCount(int ny) const {return Volt_Count_Dir[ny];} + + unsigned int GetCurrCount() const {return Curr_Count;} + unsigned int GetCurrCount(int ny) const {return Curr_Count_Dir[ny];} + protected: Operator_Ext_Excitation(Operator* op, Operator_Ext_Excitation* op_ext); Excitation* m_Exc; + + void setupVoltageExcitation( vector const volt_vIndex[3], vector const& volt_vExcit, + vector const& volt_vDelay, vector const& volt_vDir ); + void setupCurrentExcitation( vector const curr_vIndex[3], vector const& curr_vExcit, + vector const& curr_vDelay, vector const& curr_vDir ); + //E-Field/voltage Excitation + unsigned int Volt_Count; + unsigned int Volt_Count_Dir[3]; + unsigned int* Volt_index[3]; + unsigned short* Volt_dir; + FDTD_FLOAT* Volt_amp; //represented as edge-voltages!! + unsigned int* Volt_delay; + + //H-Field/current Excitation + unsigned int Curr_Count; + unsigned int Curr_Count_Dir[3]; + unsigned int* Curr_index[3]; + unsigned short* Curr_dir; + FDTD_FLOAT* Curr_amp; //represented as edge-currents!! + unsigned int* Curr_delay; }; #endif // OPERATOR_EXT_EXCITATION_H diff --git a/FDTD/extensions/operator_extension.h b/FDTD/extensions/operator_extension.h index 27ade51..ac5956a 100644 --- a/FDTD/extensions/operator_extension.h +++ b/FDTD/extensions/operator_extension.h @@ -60,6 +60,8 @@ public: virtual void ShowStat(ostream &ostr) const; + virtual void Reset() {} + protected: Operator_Extension(Operator* op); //! Copy constructor diff --git a/FDTD/openems_fdtd_mpi.cpp b/FDTD/openems_fdtd_mpi.cpp index 015dc40..c8ae4cc 100644 --- a/FDTD/openems_fdtd_mpi.cpp +++ b/FDTD/openems_fdtd_mpi.cpp @@ -443,9 +443,10 @@ void openEMS_FDTD_MPI::RunFDTD() double currE=0; //add all timesteps to end-crit field processing with max excite amplitude - unsigned int maxExcite = FDTD_Op->Exc->GetMaxExcitationTimestep(); - for (unsigned int n=0; nExc->Volt_Count; ++n) - m_ProcField->AddStep(FDTD_Op->Exc->Volt_delay[n]+maxExcite); + unsigned int maxExcite = FDTD_Op->GetExcitationSignal()->GetMaxExcitationTimestep(); +// for (unsigned int n=0; nExc->Volt_Count; ++n) +// m_ProcField->AddStep(FDTD_Op->Exc->Volt_delay[n]+maxExcite); + m_ProcField->AddStep(maxExcite); int prevTS=0,currTS=0; double speed = m_NumberCells/1e6; @@ -457,7 +458,6 @@ void openEMS_FDTD_MPI::RunFDTD() timeval prevTime= currTime; //*************** simulate ************// - PA->PreProcess(); int step = GetNextStep(); @@ -501,7 +501,7 @@ void openEMS_FDTD_MPI::RunFDTD() PA->FlushNext(); } } - if ((m_MyID==0) && (m_EnergyDecrement>endCrit) && (FDTD_Op->Exc->GetExciteType()==0)) + if ((m_MyID==0) && (m_EnergyDecrement>endCrit) && (FDTD_Op->GetExcitationSignal()->GetExciteType()==0)) cerr << "RunFDTD: max. number of timesteps was reached before the end-criteria of -" << fabs(10.0*log10(endCrit)) << "dB was reached... " << endl << \ "\tYou may want to choose a higher number of max. timesteps... " << endl; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 0a6e572..2565765 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -25,6 +25,7 @@ #include "tools/array_ops.h" #include "tools/vtk_file_writer.h" #include "fparser.hh" +#include "extensions/operator_ext_excitation.h" Operator* Operator::New() { @@ -36,7 +37,7 @@ Operator* Operator::New() Operator::Operator() : Operator_Base() { - Exc = 0; + m_Exc = 0; m_InvaildTimestep = false; m_TimeStepVar = 3; } @@ -82,7 +83,7 @@ void Operator::Init() EC_R[n]=NULL; } - Exc = 0; + m_Exc = 0; } void Operator::Delete() @@ -103,7 +104,7 @@ void Operator::Delete() delete[] EC_R[n];EC_R[n]=0; } - delete Exc;Exc=0; + delete m_Exc;m_Exc=0; Delete_N_3DArray(m_epsR,numLines); m_epsR=0; @@ -458,8 +459,8 @@ void Operator::ShowStat() const cout <<"\t(" << opt_dT << ")"; cout << endl; cout << "Timestep method name\t: " << m_Used_TS_Name << endl; - cout << "Nyquist criteria (TS)\t: " << Exc->GetNyquistNum() << endl; - cout << "Nyquist criteria (s)\t: " << Exc->GetNyquistNum()*dT << endl; + cout << "Nyquist criteria (TS)\t: " << m_Exc->GetNyquistNum() << endl; + cout << "Nyquist criteria (s)\t: " << m_Exc->GetNyquistNum()*dT << endl; cout << "-----------------------------------" << endl; } @@ -482,11 +483,32 @@ void Operator::DumpOperator2File(string filename) cout << "Operator: Dumping FDTD operator information to vtk file: " << filename << " ..." << flush; - FDTD_FLOAT**** exc = Create_N_3DArray(numLines); - if (Exc) + VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType); + vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling); + vtk_Writer->SetHeader("openEMS - Operator dump"); + + vtk_Writer->SetNativeDump(true); + + if (m_Op_Ext_Exc) { - for (unsigned int n=0; nVolt_Count; ++n) - exc[Exc->Volt_dir[n]][Exc->Volt_index[0][n]][Exc->Volt_index[1][n]][Exc->Volt_index[2][n]] = Exc->Volt_amp[n]; + FDTD_FLOAT**** exc = NULL; + if (m_Op_Ext_Exc->Volt_Count>0) + { + exc = Create_N_3DArray(numLines); + for (unsigned int n=0; nVolt_Count; ++n) + exc[m_Op_Ext_Exc->Volt_dir[n]][m_Op_Ext_Exc->Volt_index[0][n]][m_Op_Ext_Exc->Volt_index[1][n]][m_Op_Ext_Exc->Volt_index[2][n]] = m_Op_Ext_Exc->Volt_amp[n]; + vtk_Writer->AddVectorField("exc_volt",exc); + Delete_N_3DArray(exc,numLines); + } + + if (m_Op_Ext_Exc->Curr_Count>0) + { + exc = Create_N_3DArray(numLines); + for (unsigned int n=0; nCurr_Count; ++n) + exc[m_Op_Ext_Exc->Curr_dir[n]][m_Op_Ext_Exc->Curr_index[0][n]][m_Op_Ext_Exc->Curr_index[1][n]][m_Op_Ext_Exc->Curr_index[2][n]] = m_Op_Ext_Exc->Curr_amp[n]; + vtk_Writer->AddVectorField("exc_curr",exc); + Delete_N_3DArray(exc,numLines); + } } FDTD_FLOAT**** vv_temp = Create_N_3DArray(numLines); @@ -506,11 +528,6 @@ void Operator::DumpOperator2File(string filename) ii_temp[n][pos[0]][pos[1]][pos[2]] = GetII(n,pos); } - VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType); - vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling); - vtk_Writer->SetHeader("openEMS - Operator dump"); - - vtk_Writer->SetNativeDump(true); vtk_Writer->AddVectorField("vv",vv_temp); Delete_N_3DArray(vv_temp,numLines); @@ -520,13 +537,11 @@ void Operator::DumpOperator2File(string filename) Delete_N_3DArray(iv_temp,numLines); vtk_Writer->AddVectorField("ii",ii_temp); Delete_N_3DArray(ii_temp,numLines); - vtk_Writer->AddVectorField("exc",exc); - Delete_N_3DArray(exc,numLines); if (vtk_Writer->Write()==false) cerr << "Operator::DumpOperator2File: Error: Can't write file... skipping!" << endl; - cout << " done!" << endl; + delete vtk_Writer; } //! \brief dump PEC (perfect electric conductor) information (into VTK-file) @@ -615,7 +630,7 @@ void Operator::DumpPEC2File( string filename ) if (vtk_Writer->Write()==false) cerr << "Operator::DumpPEC2File: Error: Can't write file... skipping!" << endl; - cout << " done!" << endl; + delete vtk_Writer; } void Operator::DumpMaterial2File(string filename) @@ -671,7 +686,7 @@ void Operator::DumpMaterial2File(string filename) if (vtk_Writer->Write()==false) cerr << "Operator::DumpMaterial2File: Error: Can't write file... skipping!" << endl; - cout << " done!" << endl; + delete vtk_Writer; } bool Operator::SetupCSXGrid(CSRectGrid* grid) @@ -816,12 +831,13 @@ double Operator::GetDiscMaterial(int type, int n, const unsigned int pos[3]) con void Operator::InitExcitation() { - if (Exc!=NULL) - Exc->Reset(dT); + if (m_Exc!=NULL) + m_Exc->Reset(dT); else { - Exc = new Excitation( dT ); - this->AddExtension(new Operator_Ext_Excitation(this,Exc)); + m_Exc = new Excitation( dT ); + m_Op_Ext_Exc = new Operator_Ext_Excitation(this); + this->AddExtension(m_Op_Ext_Exc); } } @@ -1408,8 +1424,8 @@ bool Operator::Calc_LumpedElements() void Operator::DumpExciationSignals() { - Exc->DumpVoltageExcite("et"); - Exc->DumpCurrentExcite("ht"); + m_Exc->DumpVoltageExcite("et"); + m_Exc->DumpCurrentExcite("ht"); } void Operator::Init_EC() diff --git a/FDTD/operator.h b/FDTD/operator.h index ba5e0be..9390a2a 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -24,6 +24,7 @@ #include "Common/operator_base.h" class Operator_Extension; +class Operator_Ext_Excitation; class Engine; class TiXmlElement; @@ -51,7 +52,7 @@ public: virtual int CalcECOperator( DebugFlags debugFlags = None ); - virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) {return Exc->setupExcitation(Excite,maxTS);} + virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) {return m_Exc->setupExcitation(Excite,maxTS);} virtual void DumpExciationSignals(); @@ -84,7 +85,7 @@ public: bool GetTimestepValid() const {return !m_InvaildTimestep;} virtual double GetNumberCells() const; - virtual unsigned int GetNumberOfNyquistTimesteps() const {return Exc->GetNyquistNum();} + virtual unsigned int GetNumberOfNyquistTimesteps() const {return m_Exc->GetNyquistNum();} virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];} @@ -145,6 +146,10 @@ public: virtual double GetDiscMaterial(int type, int ny, const unsigned int pos[3]) const; + Excitation* GetExcitationSignal() const {return m_Exc;} + + Operator_Ext_Excitation* GetExcitationExtension() const {return m_Op_Ext_Exc;} + protected: //! use New() for creating a new Operator Operator(); @@ -229,6 +234,10 @@ protected: vector m_Op_exts; + // excitation classes + Excitation* m_Exc; // excitation time signal class + Operator_Ext_Excitation* m_Op_Ext_Exc; // excitation extension + // engine/post-proc needs access public: //EC operator @@ -236,8 +245,6 @@ public: FDTD_FLOAT**** vi; //calc new voltage from old current FDTD_FLOAT**** ii; //calc new current from old current FDTD_FLOAT**** iv; //calc new current from old voltage - - Excitation* Exc; }; inline Operator::DebugFlags operator|( Operator::DebugFlags a, Operator::DebugFlags b ) { return static_cast(static_cast(a) | static_cast(b)); } diff --git a/FDTD/operator_cylindermultigrid.cpp b/FDTD/operator_cylindermultigrid.cpp index ee18de1..d1b50b3 100644 --- a/FDTD/operator_cylindermultigrid.cpp +++ b/FDTD/operator_cylindermultigrid.cpp @@ -449,7 +449,7 @@ bool Operator_CylinderMultiGrid::SetupExcitation(TiXmlElement* Excite, unsigned { if (!m_InnerOp->SetupExcitation(Excite,maxTS)) return false; - return Exc->setupExcitation(Excite,maxTS); + return m_Exc->setupExcitation(Excite,maxTS); } void Operator_CylinderMultiGrid::Delete() diff --git a/openems.cpp b/openems.cpp index b76f1f9..1298b57 100644 --- a/openems.cpp +++ b/openems.cpp @@ -306,7 +306,7 @@ bool openEMS::SetupProcessing() if (g_settings.GetVerboseLevel()>0) cout << "Setting up processing..." << endl; - unsigned int Nyquist = FDTD_Op->Exc->GetNyquistNum(); + unsigned int Nyquist = FDTD_Op->GetExcitationSignal()->GetNyquistNum(); PA = new ProcessingArray(Nyquist); double start[3]; @@ -681,14 +681,15 @@ int openEMS::SetupFDTD(const char* file) FDTD_Op->ShowExtStat(); cout << "Creation time for operator: " << CalcDiffTime(OpDoneTime,startTime) << " s" << endl; } + Excitation* Exc=FDTD_Op->GetExcitationSignal(); cout << "FDTD simulation size: " << FDTD_Op->GetNumberOfLines(0) << "x" << FDTD_Op->GetNumberOfLines(1) << "x" << FDTD_Op->GetNumberOfLines(2) << " --> " << FDTD_Op->GetNumberCells() << " FDTD cells " << endl; - cout << "FDTD timestep is: " <GetTimestep() << " s; Nyquist rate: " << FDTD_Op->Exc->GetNyquistNum() << " timesteps @" << CalcNyquistFrequency(FDTD_Op->Exc->GetNyquistNum(),FDTD_Op->GetTimestep()) << " Hz" << endl; - if (FDTD_Op->Exc->GetNyquistNum()>1000) + cout << "FDTD timestep is: " <GetTimestep() << " s; Nyquist rate: " << Exc->GetNyquistNum() << " timesteps @" << CalcNyquistFrequency(Exc->GetNyquistNum(),FDTD_Op->GetTimestep()) << " Hz" << endl; + if (Exc->GetNyquistNum()>1000) cerr << "openEMS::SetupFDTD: Warning, the timestep seems to be very small --> long simulation. Check your mesh!?" << endl; - cout << "Excitation signal length is: " << FDTD_Op->Exc->Length << " timesteps (" << FDTD_Op->Exc->Length*FDTD_Op->GetTimestep() << "s)" << endl; - cout << "Max. number of timesteps: " << NrTS << " ( --> " << (double)NrTS/(double)(FDTD_Op->Exc->Length) << " * Excitation signal length)" << endl; - if ( ((double)NrTS/(double)FDTD_Op->Exc->Length < 3) && (FDTD_Op->Exc->GetExciteType()==0)) + cout << "Excitation signal length is: " << Exc->GetLength() << " timesteps (" << Exc->GetLength()*FDTD_Op->GetTimestep() << "s)" << endl; + cout << "Max. number of timesteps: " << NrTS << " ( --> " << (double)NrTS/(double)(Exc->GetLength()) << " * Excitation signal length)" << endl; + if ( ((double)NrTS/(double)Exc->GetLength() < 3) && (Exc->GetExciteType()==0)) cerr << "openEMS::SetupFDTD: Warning, max. number of timesteps is smaller than three times the excitation. " << endl << \ "\tYou may want to choose a higher number of max. timesteps... " << endl; @@ -767,9 +768,10 @@ void openEMS::RunFDTD() PA->InitAll(); //add all timesteps to end-crit field processing with max excite amplitude - unsigned int maxExcite = FDTD_Op->Exc->GetMaxExcitationTimestep(); - for (unsigned int n=0; nExc->Volt_Count; ++n) - ProcField->AddStep(FDTD_Op->Exc->Volt_delay[n]+maxExcite); + unsigned int maxExcite = FDTD_Op->GetExcitationSignal()->GetMaxExcitationTimestep(); +// for (unsigned int n=0; nExc->Volt_Count; ++n) +// ProcField->AddStep(FDTD_Op->Exc->Volt_delay[n]+maxExcite); + ProcField->AddStep(maxExcite); double change=1; int prevTS=0,currTS=0; @@ -821,7 +823,7 @@ void openEMS::RunFDTD() PA->FlushNext(); } } - if ((change>endCrit) && (FDTD_Op->Exc->GetExciteType()==0)) + if ((change>endCrit) && (FDTD_Op->GetExcitationSignal()->GetExciteType()==0)) cerr << "RunFDTD: Warning: Max. number of timesteps was reached before the end-criteria of -" << fabs(10.0*log10(endCrit)) << "dB was reached... " << endl << \ "\tYou may want to choose a higher number of max. timesteps... " << endl;