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

excitation moved to its own extension

Browse Source 
the benefit of this approach is, that the excitation
is subject to the extensions priorities
This commit is contained in:
Thorsten Liebig 2011-03-16 16:26:01 +01:00
parent a1a48b1c23
commit aba33c58bd
17 changed files with 517 additions and 253 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
FDTD/engine.cpp
View File

@ -86,6 +86,13 @@ void Engine::SortExtensionByPriority()
{
stable_sort(m_Eng_exts.begin(),m_Eng_exts.end(), CompareExtensions);
reverse(m_Eng_exts.begin(),m_Eng_exts.end());
if (g_settings.GetVerboseLevel()>0)
{
cout << "--- Engine::SortExtensionByPriority() ---" << endl;
for (size_t n=0; n<m_Eng_exts.size(); ++n)
cout << " #" << n << ": " << m_Eng_exts.at(n)->GetExtensionName() << " (" << m_Eng_exts.at(n)->GetPriority() << ")" << endl;
}
}
void Engine::Reset()
@ -132,25 +139,6 @@ void Engine::UpdateVoltages(unsigned int startX, unsigned int numX)
}
}
void Engine::ApplyVoltageExcite()
{
int exc_pos;
unsigned int ny;
unsigned int pos[3];
//soft voltage excitation here (E-field excite)
for (unsigned int n=0; n<Op->Exc->Volt_Count; ++n)
{
exc_pos = (int)numTS - (int)Op->Exc->Volt_delay[n];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
// if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl;
ny = Op->Exc->Volt_dir[n];
pos[0]=Op->Exc->Volt_index[0][n];
pos[1]=Op->Exc->Volt_index[1][n];
pos[2]=Op->Exc->Volt_index[2][n];
SetVolt(ny,pos, GetVolt(ny,pos) + Op->Exc->Volt_amp[n]*Op->Exc->Signal_volt[exc_pos]);
}
}
void Engine::UpdateCurrents(unsigned int startX, unsigned int numX)
{
unsigned int pos[3];
@ -179,25 +167,6 @@ void Engine::UpdateCurrents(unsigned int startX, unsigned int numX)
}
}
void Engine::ApplyCurrentExcite()
{
int exc_pos;
unsigned int ny;
unsigned int pos[3];
//soft current excitation here (H-field excite)
for (unsigned int n=0; n<Op->Exc->Curr_Count; ++n)
{
exc_pos = (int)numTS - (int)Op->Exc->Curr_delay[n];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
// if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl;
ny = Op->Exc->Curr_dir[n];
pos[0]=Op->Exc->Curr_index[0][n];
pos[1]=Op->Exc->Curr_index[1][n];
pos[2]=Op->Exc->Curr_index[2][n];
SetCurr(ny,pos, GetCurr(ny,pos) + Op->Exc->Curr_amp[n]*Op->Exc->Signal_curr[exc_pos]);
}
}
void Engine::DoPreVoltageUpdates()
{
//execute extensions in reverse order -> highest priority gets access to the voltages last
@ -250,14 +219,12 @@ bool Engine::IterateTS(unsigned int iterTS)
UpdateVoltages(0,numLines[0]);
DoPostVoltageUpdates();
Apply2Voltages();
ApplyVoltageExcite();
//current updates with extensions
DoPreCurrentUpdates();
UpdateCurrents(0,numLines[0]-1);
DoPostCurrentUpdates();
Apply2Current();
ApplyCurrentExcite();
++numTS;
}

4
FDTD/engine.h
View File

@ -66,8 +66,6 @@ public:
virtual void DoPostVoltageUpdates();
//! Apply extension voltage changes
virtual void Apply2Voltages();
//! Apply voltage excitations
virtual void ApplyVoltageExcite();
//! Execute Pre-Current extension updates
virtual void DoPreCurrentUpdates();
@ -77,8 +75,6 @@ public:
virtual void DoPostCurrentUpdates();
//! Apply extension current changes
virtual void Apply2Current();
//! Apply current excitations
virtual void ApplyCurrentExcite();
inline size_t GetExtensionCount() {return m_Eng_exts.size();}
inline Engine_Extension* GetExtension(size_t nr) {return m_Eng_exts.at(nr);}

2
FDTD/engine_mpi.cpp
View File

@ -201,7 +201,6 @@ bool Engine_MPI::IterateTS(unsigned int iterTS)
UpdateVoltages(0,numLines[0]);
DoPostVoltageUpdates();
Apply2Voltages();
ApplyVoltageExcite();
SendReceiveVoltages();
//current updates with extensions
@ -209,7 +208,6 @@ bool Engine_MPI::IterateTS(unsigned int iterTS)
UpdateCurrents(0,numLines[0]-1);
DoPostCurrentUpdates();
Apply2Current();
ApplyCurrentExcite();
SendReceiveCurrents();
++numTS;

13
FDTD/engine_multithread.cpp
View File

@ -278,12 +278,6 @@ void thread::operator()()
m_enginePtr->DoPostVoltageUpdates(m_threadID);
m_enginePtr->Apply2Voltages(m_threadID);
// voltage excitation (E-field excite) by the first thread
if (m_threadID==0)
m_enginePtr->ApplyVoltageExcite();
m_enginePtr->m_IterateBarrier->wait();
// voltage excitation finished
#ifdef MPI_SUPPORT
if (m_threadID==0)
{
@ -314,13 +308,6 @@ void thread::operator()()
m_enginePtr->DoPostCurrentUpdates(m_threadID);
m_enginePtr->Apply2Current(m_threadID);
// current excitation (H-field excite) by the first thread
if (m_threadID==0)
m_enginePtr->ApplyCurrentExcite();
m_enginePtr->m_IterateBarrier->wait();
// current excitation finished
#ifdef MPI_SUPPORT
if (m_threadID==0)
{

3
FDTD/excitation.h
View File

@ -49,6 +49,9 @@ public:
//! Dump current excitation signal to ASCII file
void DumpCurrentExcite(string filename);
//! Get the used timestep
double GetTimestep() const {return dT;}
//Excitation time-signal
unsigned int Length;
FDTD_FLOAT* Signal_volt;

2
FDTD/extensions/engine_ext_cylinder.cpp
View File

@ -30,7 +30,7 @@ Engine_Ext_Cylinder::Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext) : Engine
numLines[n] = op_ext->m_Op->GetOriginalNumLines(n);
//this cylindrical extension should be executed first?
m_Priority = 2000;
m_Priority = ENG_EXT_PRIO_CYLINDER;
}
void Engine_Ext_Cylinder::Apply2Voltages()

2
FDTD/extensions/engine_ext_cylindermultigrid.cpp
View File

@ -25,7 +25,7 @@ Engine_Ext_CylinderMultiGrid::Engine_Ext_CylinderMultiGrid(Operator_Extension* o
m_Eng_MG = NULL;
// the multi-grid should be applies last?
m_Priority = -1000;
m_Priority = ENG_EXT_PRIO_CYLINDERMULTIGRID;
}
Engine_Ext_CylinderMultiGrid::~Engine_Ext_CylinderMultiGrid()

146
FDTD/extensions/engine_ext_excitation.cpp Normal file
View File

@ -0,0 +1,146 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "engine_ext_excitation.h"
#include "operator_ext_excitation.h"
#include "FDTD/engine_sse.h"
Engine_Ext_Excitation::Engine_Ext_Excitation(Operator_Ext_Excitation* op_ext) : Engine_Extension(op_ext)
{
m_Op_Exc = op_ext;
m_Priority = ENG_EXT_PRIO_EXCITATION;
}
Engine_Ext_Excitation::~Engine_Ext_Excitation()
{
}
void Engine_Ext_Excitation::Apply2Voltages()
{
//soft voltage excitation here (E-field excite)
int exc_pos;
unsigned int ny;
unsigned int pos[3];
int numTS = m_Eng->GetNumberOfTimesteps();
//switch for different engine types to access faster inline engine functions
switch (m_Eng->GetType())
{
case Engine::BASIC:
{
for (unsigned int n=0; n<m_Op_Exc->m_Exc->Volt_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]);
}
break;
}
case Engine::SSE:
{
for (unsigned int n=0; n<m_Op_Exc->m_Exc->Volt_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]);
}
break;
}
default:
{
for (unsigned int n=0; n<m_Op_Exc->m_Exc->Volt_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]);
}
break;
}
}
}
void Engine_Ext_Excitation::Apply2Current()
{
//soft current excitation here (H-field excite)
int exc_pos;
unsigned int ny;
unsigned int pos[3];
int numTS = m_Eng->GetNumberOfTimesteps();
//switch for different engine types to access faster inline engine functions
switch (m_Eng->GetType())
{
case Engine::BASIC:
{
for (unsigned int n=0; n<m_Op_Exc->m_Exc->Volt_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]);
}
break;
}
case Engine::SSE:
{
for (unsigned int n=0; n<m_Op_Exc->m_Exc->Volt_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]);
}
break;
}
default:
{
for (unsigned int n=0; n<m_Op_Exc->m_Exc->Volt_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]);
}
break;
}
}
}

40
FDTD/extensions/engine_ext_excitation.h Normal file
View File

@ -0,0 +1,40 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef ENGINE_EXT_EXCITATION_H
#define ENGINE_EXT_EXCITATION_H
#include "engine_extension.h"
#include "FDTD/engine.h"
#include "FDTD/operator.h"
class Operator_Ext_Excitation;
class Engine_Ext_Excitation : public Engine_Extension
{
public:
Engine_Ext_Excitation(Operator_Ext_Excitation* op_ext);
virtual ~Engine_Ext_Excitation();
virtual void Apply2Voltages();
virtual void Apply2Current();
protected:
Operator_Ext_Excitation* m_Op_Exc;
};
#endif // ENGINE_EXT_EXCITATION_H

2
FDTD/extensions/engine_ext_upml.cpp
View File

@ -27,7 +27,7 @@ Engine_Ext_UPML::Engine_Ext_UPML(Operator_Ext_UPML* op_ext) : Engine_Extension(o
m_Op_UPML = op_ext;
//this ABC extension should be executed first!
m_Priority = 1e6;
m_Priority = ENG_EXT_PRIO_UPML;
volt_flux = Create_N_3DArray<FDTD_FLOAT>(m_Op_UPML->m_numLines);
curr_flux = Create_N_3DArray<FDTD_FLOAT>(m_Op_UPML->m_numLines);

11
FDTD/extensions/engine_extension.cpp
View File

@ -16,6 +16,7 @@
*/
#include "engine_extension.h"
#include "operator_extension.h"
#include "FDTD/engine.h"
@ -23,7 +24,7 @@ Engine_Extension::Engine_Extension(Operator_Extension* op_ext)
{
m_Op_ext = op_ext;
m_Eng = NULL;
m_Priority = 0;
m_Priority = ENG_EXT_PRIO_DEFAULT;
m_NrThreads = 1;
}
@ -38,6 +39,14 @@ void Engine_Extension::SetNumberOfThreads(int nrThread)
m_NrThreads=nrThread;
}
string Engine_Extension::GetExtensionName() const
{
if (m_Op_ext)
return m_Op_ext->GetExtensionName();
else
return "Unknown Extension";
}
void Engine_Extension::DoPreVoltageUpdates(int threadID)
{
//if this method gets called the derived extension obviously doesn't support multithrading, calling non-MT method...

12
FDTD/extensions/engine_extension.h
View File

@ -18,6 +18,16 @@
#ifndef ENGINE_EXTENSION_H
#define ENGINE_EXTENSION_H
#define ENG_EXT_PRIO_DEFAULT 0 //default engine extension priority
// priority definitions for some important extensions
#define ENG_EXT_PRIO_UPML 1e6; //unaxial pml extension priority
#define ENG_EXT_PRIO_EXCITATION 3000 //excitation priority
#define ENG_EXT_PRIO_CYLINDER 2000 //cylindrial extension priority
#define ENG_EXT_PRIO_CYLINDERMULTIGRID -1000; //cylindrial multi-grid extension priority
#include <string>
class Operator_Extension;
class Engine;
@ -60,6 +70,8 @@ public:
virtual bool operator< (const Engine_Extension& other);
virtual std::string GetExtensionName() const;
protected:
Engine_Extension(Operator_Extension* op_ext);

233
FDTD/extensions/operator_ext_excitation.cpp Normal file
View File

@ -0,0 +1,233 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "operator_ext_excitation.h"
#include "engine_ext_excitation.h"
#include "FDTD/excitation.h"
#include "ContinuousStructure.h"
Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op, Excitation* exc) : Operator_Extension(op)
{
m_Exc = exc;
}
Operator_Ext_Excitation::~Operator_Ext_Excitation()
{
}
Operator_Ext_Excitation::Operator_Ext_Excitation(Operator* op, Operator_Ext_Excitation* op_ext) : Operator_Extension(op, op_ext)
{
m_Exc = NULL;
}
bool Operator_Ext_Excitation::BuildExtension()
{
double dT = m_Exc->GetTimestep();
if (dT==0)
return false;
if (m_Exc==0)
return false;
ContinuousStructure* CSX = m_Op->GetGeometryCSX();
unsigned int pos[3];
double amp=0;
vector<unsigned int> volt_vIndex[3];
vector<FDTD_FLOAT> volt_vExcit;
vector<unsigned int> volt_vDelay;
vector<unsigned int> volt_vDir;
double volt_coord[3];
vector<unsigned int> curr_vIndex[3];
vector<FDTD_FLOAT> curr_vExcit;
vector<unsigned int> curr_vDelay;
vector<unsigned int> curr_vDir;
double curr_coord[3];
vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::ELECTRODE);
if (vec_prop.size()==0)
{
cerr << "Operator::CalcFieldExcitation: Warning, no excitation properties found" << endl;
return false;
}
CSPropElectrode* elec=NULL;
CSProperties* prop=NULL;
int priority=0;
unsigned int numLines[] = {m_Op->GetOriginalNumLines(0),m_Op->GetOriginalNumLines(1),m_Op->GetOriginalNumLines(2)};
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
//electric field excite
for (int n=0; n<3; ++n)
{
m_Op->GetYeeCoords(n,pos,volt_coord,false);
for (size_t p=0; p<vec_prop.size(); ++p)
{
prop = vec_prop.at(p);
elec = prop->ToElectrode();
if (elec==NULL)
continue;
if (prop->CheckCoordInPrimitive(volt_coord,priority,true))
{
if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))//&& (pos[n]<numLines[n]-1))
{
amp = elec->GetWeightedExcitation(n,volt_coord)*m_Op->GetEdgeLength(n,pos);// delta[n]*gridDelta;
if (amp!=0)
{
volt_vExcit.push_back(amp);
volt_vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
volt_vDir.push_back(n);
volt_vIndex[0].push_back(pos[0]);
volt_vIndex[1].push_back(pos[1]);
volt_vIndex[2].push_back(pos[2]);
}
if (elec->GetExcitType()==1) //hard excite
{
m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 );
}
}
}
}
}
//magnetic field excite
for (int n=0; n<3; ++n)
{
if ((pos[0]>=numLines[0]-1) || (pos[1]>=numLines[1]-1) || (pos[2]>=numLines[2]-1))
continue; //skip the last H-Line which is outside the FDTD-domain
m_Op->GetYeeCoords(n,pos,curr_coord,true);
for (size_t p=0; p<vec_prop.size(); ++p)
{
prop = vec_prop.at(p);
elec = prop->ToElectrode();
if (elec==NULL)
continue;
if (prop->CheckCoordInPrimitive(curr_coord,priority,true))
{
if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==2) || (elec->GetExcitType()==3) ))
{
amp = elec->GetWeightedExcitation(n,curr_coord)*m_Op->GetEdgeLength(n,pos,true);// delta[n]*gridDelta;
if (amp!=0)
{
curr_vExcit.push_back(amp);
curr_vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
curr_vDir.push_back(n);
curr_vIndex[0].push_back(pos[0]);
curr_vIndex[1].push_back(pos[1]);
curr_vIndex[2].push_back(pos[2]);
}
if (elec->GetExcitType()==3) //hard excite
{
m_Op->SetII(n,pos[0],pos[1],pos[2], 0 );
m_Op->SetIV(n,pos[0],pos[1],pos[2], 0 );
}
}
}
}
}
}
}
}
//special treatment for primitives of type curve (treated as wires) see also Calc_PEC
double p1[3];
double p2[3];
struct Operator::Grid_Path path;
for (size_t p=0; p<vec_prop.size(); ++p)
{
prop = vec_prop.at(p);
elec = prop->ToElectrode();
for (size_t n=0; n<prop->GetQtyPrimitives(); ++n)
{
CSPrimitives* prim = prop->GetPrimitive(n);
CSPrimCurve* curv = prim->ToCurve();
if (curv)
{
for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
{
curv->GetPoint(i-1,p1);
curv->GetPoint(i,p2);
path = m_Op->FindPath(p1,p2);
if (path.dir.size()>0)
prim->SetPrimitiveUsed(true);
for (size_t t=0; t<path.dir.size(); ++t)
{
n = path.dir.at(t);
pos[0] = path.posPath[0].at(t);
pos[1] = path.posPath[1].at(t);
pos[2] = path.posPath[2].at(t);
m_Op->GetYeeCoords(n,pos,volt_coord,false);
if (elec!=NULL)
{
if ((elec->GetActiveDir(n)) && (pos[n]<numLines[n]-1) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))
{
amp = elec->GetWeightedExcitation(n,volt_coord)*m_Op->GetEdgeLength(n,pos);
if (amp!=0)
{
volt_vExcit.push_back(amp);
volt_vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
volt_vDir.push_back(n);
volt_vIndex[0].push_back(pos[0]);
volt_vIndex[1].push_back(pos[1]);
volt_vIndex[2].push_back(pos[2]);
}
if (elec->GetExcitType()==1) //hard excite
{
m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
m_Op->SetVI(n,pos[0],pos[1],pos[2], 0 );
}
}
}
}
}
}
}
}
// set voltage excitations
m_Exc->setupVoltageExcitation( volt_vIndex, volt_vExcit, volt_vDelay, volt_vDir );
// set current excitations
m_Exc->setupCurrentExcitation( curr_vIndex, curr_vExcit, curr_vDelay, curr_vDir );
return true;
}
Engine_Extension* Operator_Ext_Excitation::CreateEngineExtention()
{
return new Engine_Ext_Excitation(this);
}
void Operator_Ext_Excitation::ShowStat(ostream &ostr) const
{
Operator_Extension::ShowStat(ostr);
cout << "Voltage excitations\t: " << m_Exc->Volt_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;
}

52
FDTD/extensions/operator_ext_excitation.h Normal file
View File

@ -0,0 +1,52 @@
/*
* Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPERATOR_EXT_EXCITATION_H
#define OPERATOR_EXT_EXCITATION_H
#include "operator_extension.h"
#include "FDTD/operator.h"
class Excitation;
class Operator_Ext_Excitation : public Operator_Extension
{
friend class Engine_Ext_Excitation;
public:
Operator_Ext_Excitation(Operator* op, Excitation* exc);
~Operator_Ext_Excitation();
virtual Operator_Extension* Clone(Operator* op) {return NULL;}
virtual bool BuildExtension();
virtual Engine_Extension* CreateEngineExtention();
virtual bool IsCylinderCoordsSave() const {return true;}
virtual bool IsCylindricalMultiGridSave(bool child) const {return true;}
virtual string GetExtensionName() const {return string("Excitation Extension");}
virtual void ShowStat(ostream &ostr) const;
protected:
Operator_Ext_Excitation(Operator* op, Operator_Ext_Excitation* op_ext);
Excitation* m_Exc;
};
#endif // OPERATOR_EXT_EXCITATION_H

187
FDTD/operator.cpp
View File

@ -19,6 +19,7 @@
#include "operator.h"
#include "engine.h"
#include "extensions/operator_extension.h"
#include "extensions/operator_ext_excitation.h"
#include "Common/processfields.h"
#include "tools/array_ops.h"
#include "fparser.hh"
@ -351,9 +352,6 @@ void Operator::ShowStat() const
cout << "Size of Operator\t: " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MiB) " << endl;
cout << "Size of Field-Data\t: " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MiB) " << endl;
cout << "-----------------------------------" << endl;
cout << "Voltage excitations\t: " << Exc->Volt_Count << "\t (" << Exc->Volt_Count_Dir[0] << ", " << Exc->Volt_Count_Dir[1] << ", " << Exc->Volt_Count_Dir[2] << ")" << endl;
cout << "Current excitations\t: " << Exc->Curr_Count << "\t (" << Exc->Curr_Count_Dir[0] << ", " << Exc->Curr_Count_Dir[1] << ", " << Exc->Curr_Count_Dir[2] << ")" << endl;
cout << "-----------------------------------" << endl;
cout << "Number of PEC edges\t: " << m_Nr_PEC[0]+m_Nr_PEC[1]+m_Nr_PEC[2] << endl;
cout << "in " << GetDirName(0) << " direction\t\t: " << m_Nr_PEC[0] << endl;
cout << "in " << GetDirName(1) << " direction\t\t: " << m_Nr_PEC[1] << endl;
@ -366,8 +364,6 @@ void Operator::ShowStat() const
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 << "Excitation Length (TS)\t: " << Exc->Length << endl;
cout << "Excitation Length (s)\t: " << Exc->Length*dT << endl;
cout << "-----------------------------------" << endl;
}
@ -717,6 +713,8 @@ void Operator::InitExcitation()
{
delete Exc;
Exc = new Excitation( dT );
this->AddExtension(new Operator_Ext_Excitation(this,Exc));
}
void Operator::Calc_ECOperatorPos(int n, unsigned int* pos)
@ -804,7 +802,6 @@ int Operator::CalcECOperator( DebugFlags debugFlags )
ApplyMagneticBC(PMC);
InitExcitation();
CalcFieldExcitation();
//all information available for extension... create now...
for (size_t n=0; n<m_Op_exts.size(); ++n)
@ -1352,184 +1349,6 @@ double Operator::CalcTimestep_Var3()
return 0;
}
bool Operator::CalcFieldExcitation()
{
if (dT==0)
return false;
if (Exc==0)
return false;
unsigned int pos[3];
double amp=0;
vector<unsigned int> volt_vIndex[3];
vector<FDTD_FLOAT> volt_vExcit;
vector<unsigned int> volt_vDelay;
vector<unsigned int> volt_vDir;
double volt_coord[3];
vector<unsigned int> curr_vIndex[3];
vector<FDTD_FLOAT> curr_vExcit;
vector<unsigned int> curr_vDelay;
vector<unsigned int> curr_vDir;
double curr_coord[3];
vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::ELECTRODE);
if (vec_prop.size()==0)
{
cerr << "Operator::CalcFieldExcitation: Warning, no excitation properties found" << endl;
return false;
}
CSPropElectrode* elec=NULL;
CSProperties* prop=NULL;
int priority=0;
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
//electric field excite
for (int n=0; n<3; ++n)
{
GetYeeCoords(n,pos,volt_coord,false);
for (size_t p=0; p<vec_prop.size(); ++p)
{
prop = vec_prop.at(p);
elec = prop->ToElectrode();
if (elec==NULL)
continue;
if (prop->CheckCoordInPrimitive(volt_coord,priority,true))
{
if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))//&& (pos[n]<numLines[n]-1))
{
amp = elec->GetWeightedExcitation(n,volt_coord)*GetEdgeLength(n,pos);// delta[n]*gridDelta;
if (amp!=0)
{
volt_vExcit.push_back(amp);
volt_vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
volt_vDir.push_back(n);
volt_vIndex[0].push_back(pos[0]);
volt_vIndex[1].push_back(pos[1]);
volt_vIndex[2].push_back(pos[2]);
}
if (elec->GetExcitType()==1) //hard excite
{
SetVV(n,pos[0],pos[1],pos[2], 0 );
SetVI(n,pos[0],pos[1],pos[2], 0 );
}
}
}
}
}
//magnetic field excite
for (int n=0; n<3; ++n)
{
if ((pos[0]>=numLines[0]-1) || (pos[1]>=numLines[1]-1) || (pos[2]>=numLines[2]-1))
continue; //skip the last H-Line which is outside the FDTD-domain
GetYeeCoords(n,pos,curr_coord,true);
for (size_t p=0; p<vec_prop.size(); ++p)
{
prop = vec_prop.at(p);
elec = prop->ToElectrode();
if (elec==NULL)
continue;
if (prop->CheckCoordInPrimitive(curr_coord,priority,true))
{
if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==2) || (elec->GetExcitType()==3) ))
{
amp = elec->GetWeightedExcitation(n,curr_coord)*GetEdgeLength(n,pos,true);// delta[n]*gridDelta;
if (amp!=0)
{
curr_vExcit.push_back(amp);
curr_vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
curr_vDir.push_back(n);
curr_vIndex[0].push_back(pos[0]);
curr_vIndex[1].push_back(pos[1]);
curr_vIndex[2].push_back(pos[2]);
}
if (elec->GetExcitType()==3) //hard excite
{
SetII(n,pos[0],pos[1],pos[2], 0 );
SetIV(n,pos[0],pos[1],pos[2], 0 );
}
}
}
}
}
}
}
}
//special treatment for primitives of type curve (treated as wires) see also Calc_PEC
double p1[3];
double p2[3];
struct Grid_Path path;
for (size_t p=0; p<vec_prop.size(); ++p)
{
prop = vec_prop.at(p);
elec = prop->ToElectrode();
for (size_t n=0; n<prop->GetQtyPrimitives(); ++n)
{
CSPrimitives* prim = prop->GetPrimitive(n);
CSPrimCurve* curv = prim->ToCurve();
if (curv)
{
for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
{
curv->GetPoint(i-1,p1);
curv->GetPoint(i,p2);
path = FindPath(p1,p2);
if (path.dir.size()>0)
prim->SetPrimitiveUsed(true);
for (size_t t=0; t<path.dir.size(); ++t)
{
n = path.dir.at(t);
pos[0] = path.posPath[0].at(t);
pos[1] = path.posPath[1].at(t);
pos[2] = path.posPath[2].at(t);
GetYeeCoords(n,pos,volt_coord,false);
if (elec!=NULL)
{
if ((elec->GetActiveDir(n)) && (pos[n]<numLines[n]-1) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))
{
amp = elec->GetWeightedExcitation(n,volt_coord)*GetEdgeLength(n,pos);
if (amp!=0)
{
volt_vExcit.push_back(amp);
volt_vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
volt_vDir.push_back(n);
volt_vIndex[0].push_back(pos[0]);
volt_vIndex[1].push_back(pos[1]);
volt_vIndex[2].push_back(pos[2]);
}
if (elec->GetExcitType()==1) //hard excite
{
SetVV(n,pos[0],pos[1],pos[2], 0 );
SetVI(n,pos[0],pos[1],pos[2], 0 );
}
}
}
}
}
}
}
}
// set voltage excitations
Exc->setupVoltageExcitation( volt_vIndex, volt_vExcit, volt_vDelay, volt_vDir );
// set current excitations
Exc->setupCurrentExcitation( curr_vIndex, curr_vExcit, curr_vDelay, curr_vDir );
return true;
}
bool Operator::CalcPEC()
{
m_Nr_PEC[0]=0;

4
FDTD/operator.h
View File

@ -34,6 +34,7 @@ class Operator : public Operator_Base
friend class Engine_Interface_FDTD;
friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice
friend class Operator_Ext_PML_SF_Plane;
friend class Operator_Ext_Excitation;
public:
enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4};
@ -153,9 +154,6 @@ protected:
};
struct Grid_Path FindPath(double start[], double stop[]);
//! Calculate the field excitations.
virtual bool CalcFieldExcitation();
// debug
virtual void DumpOperator2File(string filename);
virtual void DumpMaterial2File(string filename);

10
openEMS.pro
View File

@ -84,7 +84,9 @@ SOURCES += FDTD/extensions/engine_extension.cpp \
FDTD/extensions/engine_ext_mur_abc.cpp \
FDTD/extensions/operator_ext_mur_abc.cpp \
FDTD/extensions/operator_ext_cylinder.cpp \
FDTD/extensions/engine_ext_cylinder.cpp
FDTD/extensions/engine_ext_cylinder.cpp \
FDTD/extensions/operator_ext_excitation.cpp \
FDTD/extensions/engine_ext_excitation.cpp
# Common source files
SOURCES += Common/operator_base.cpp \
@ -123,7 +125,7 @@ HEADERS += tools/ErrorMsg.h \
FDTD/engine_interface_fdtd.h
# FDTD/extensions header files
HEADERS += FDTD/operator_extension.h \
HEADERS += FDTD/extensions/operator_extension.h \
FDTD/extensions/engine_extension.h \
FDTD/extensions/engine_ext_mur_abc.h \
FDTD/extensions/operator_ext_mur_abc.h \
@ -137,7 +139,9 @@ HEADERS += FDTD/operator_extension.h \
FDTD/extensions/engine_ext_pml_sf.h \
FDTD/extensions/engine_ext_cylindermultigrid.h \
FDTD/extensions/operator_ext_upml.h \
FDTD/extensions/engine_ext_upml.h
FDTD/extensions/engine_ext_upml.h \
FDTD/extensions/operator_ext_excitation.h \
FDTD/extensions/engine_ext_excitation.h
# Common header files
HEADERS += Common/operator_base.h \