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

split excitation from operator 

The excitation variables and functions are separated into a class Excitation.
This allows completely different operator implementations (e.g. TLM) to use the excitations functions without deriving from class Operator.
pull/1/head
Sebastian Held 2010-05-03 18:33:14 +02:00
parent 38c03d5b0c
commit f762214eae
11 changed files with 498 additions and 378 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
FDTD/engine.cpp
View File

@ -113,19 +113,19 @@ void Engine::ApplyVoltageExcite()
{
int exc_pos;
//soft voltage excitation here (E-field excite)
for (unsigned int n=0;n<Op->E_Exc_Count;++n)
for (unsigned int n=0;n<Op->Exc->E_Count;++n)
{
exc_pos = (int)numTS - (int)Op->E_Exc_delay[n];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->ExciteLength);
exc_pos = (int)numTS - (int)Op->Exc->E_delay[n];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
// if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl;
GetVolt(Op->E_Exc_dir[n],Op->E_Exc_index[0][n],Op->E_Exc_index[1][n],Op->E_Exc_index[2][n]) += Op->E_Exc_amp[n]*Op->ExciteSignal_volt[exc_pos];
GetVolt(Op->Exc->E_dir[n],Op->Exc->E_index[0][n],Op->Exc->E_index[1][n],Op->Exc->E_index[2][n]) += Op->Exc->E_amp[n]*Op->Exc->Signal_volt[exc_pos];
}
// write the first excitation into the file "et1"
if (Op->E_Exc_Count >= 1) {
exc_pos = (int)numTS - (int)Op->E_Exc_delay[0];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->ExciteLength);
file_et1 << numTS * Op->GetTimestep() << "\t" << Op->E_Exc_amp[0]*Op->ExciteSignal_volt[exc_pos] << "\n"; // do not use std::endl here, because it will do an implicit flush
if (Op->Exc->E_Count >= 1) {
exc_pos = (int)numTS - (int)Op->Exc->E_delay[0];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
file_et1 << numTS * Op->GetTimestep() << "\t" << Op->Exc->E_amp[0]*Op->Exc->Signal_volt[exc_pos] << "\n"; // do not use std::endl here, because it will do an implicit flush
}
}
@ -159,12 +159,12 @@ void Engine::ApplyCurrentExcite()
{
int exc_pos;
//soft current excitation here (H-field excite)
for (unsigned int n=0;n<Op->Curr_Exc_Count;++n)
for (unsigned int n=0;n<Op->Exc->Curr_Count;++n)
{
exc_pos = (int)numTS - (int)Op->Curr_Exc_delay[n];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->ExciteLength);
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;
GetCurr(Op->Curr_Exc_dir[n],Op->Curr_Exc_index[0][n],Op->Curr_Exc_index[1][n],Op->Curr_Exc_index[2][n]) += Op->Curr_Exc_amp[n]*Op->ExciteSignal_curr[exc_pos];
GetCurr(Op->Exc->Curr_dir[n],Op->Exc->Curr_index[0][n],Op->Exc->Curr_index[1][n],Op->Exc->Curr_index[2][n]) += Op->Exc->Curr_amp[n]*Op->Exc->Signal_curr[exc_pos];
}
}

10
FDTD/engine_ext_mur_abc.cpp
View File

@ -39,18 +39,18 @@ 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;n<m_Op_mur->m_Op->E_Exc_Count;++n)
for (unsigned int n=0;n<m_Op_mur->m_Op->Exc->E_Count;++n)
{
if ( ((m_Op_mur->m_Op->E_Exc_dir[n]==m_nyP) || (m_Op_mur->m_Op->E_Exc_dir[n]==m_nyPP)) && (m_Op_mur->m_Op->E_Exc_index[m_ny][n]==m_LineNr) )
if ( ((m_Op_mur->m_Op->Exc->E_dir[n]==m_nyP) || (m_Op_mur->m_Op->Exc->E_dir[n]==m_nyPP)) && (m_Op_mur->m_Op->Exc->E_index[m_ny][n]==m_LineNr) )
{
if ((int)m_Op_mur->m_Op->E_Exc_delay[n]>maxDelay)
maxDelay = (int)m_Op_mur->m_Op->E_Exc_delay[n];
if ((int)m_Op_mur->m_Op->Exc->E_delay[n]>maxDelay)
maxDelay = (int)m_Op_mur->m_Op->Exc->E_delay[n];
}
}
m_start_TS = 0;
if (maxDelay>=0)
{
m_start_TS = maxDelay + m_Op_mur->m_Op->ExciteLength + 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->Exc->Length + 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;
}
}

296
FDTD/excitation.cpp Normal file
View File

@ -0,0 +1,296 @@
/*
* Copyright (C) 2010 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 "tools/array_ops.h"
#include "fparser.hh"
#include "tinyxml.h"
#include "excitation.h"
Excitation::Excitation( double timestep )
{
Signal_volt = 0;
Signal_curr = 0;
E_delay = 0;
E_amp = 0;
E_dir = 0;
Curr_delay = 0;
Curr_amp = 0;
Curr_dir = 0;
for (int n=0;n<3;++n) {
E_index[n] = 0;
Curr_index[n] = 0;
}
dT = timestep;
m_nyquistTS = nan("");
}
Excitation::~Excitation()
{
delete[] Signal_volt;
delete[] Signal_curr;
delete[] E_delay;
delete[] E_dir;
delete[] E_amp;
delete[] Curr_delay;
delete[] Curr_dir;
delete[] Curr_amp;
for (int n=0;n<3;++n) {
delete[] E_index[n];
delete[] Curr_index[n];
}
}
bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS )
{
if (!Excite) {
cerr << "Can't read openEMS excitation settings... " << endl;
return false;
}
int Excit_Type=0;
double f0=0;
double fc=0;
Excite->QueryIntAttribute("Type",&Excit_Type);
switch (Excit_Type)
{
case 0:
Excite->QueryDoubleAttribute("f0",&f0);
Excite->QueryDoubleAttribute("fc",&fc);
CalcGaussianPulsExcitation(f0,fc);
break;
case 1:
Excite->QueryDoubleAttribute("f0",&f0);
CalcSinusExcitation(f0,maxTS);
break;
case 2:
CalcDiracPulsExcitation();
break;
case 3:
CalcStepExcitation();
break;
case 10:
Excite->QueryDoubleAttribute("f0",&f0);
CalcCustomExcitation(f0,maxTS,Excite->Attribute("Function"));
break;
}
if (GetNyquistNum() == 0) {
cerr << "openEMS: excitation setup failed!!" << endl;
return false;
}
return true;
}
unsigned int Excitation::CalcNyquistNum(double fmax) const
{
if (dT==0) return 1;
double T0 = 1/fmax;
return floor(T0/2/dT);
}
unsigned int Excitation::GetMaxExcitationTimestep() const
{
FDTD_FLOAT maxAmp=0;
unsigned int maxStep=0;
for (unsigned int n=1;n<Length+1;++n)
{
if (fabs(Signal_volt[n])>maxAmp)
{
maxAmp = fabs(Signal_volt[n]);
maxStep = n;
}
}
return maxStep;
}
void Excitation::CalcGaussianPulsExcitation(double f0, double fc)
{
if (dT==0) return;
Length = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
// cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0;
Signal_curr[0]=0.0;
for (unsigned int n=1;n<Length+1;++n)
{
double t = (n-1)*dT;
Signal_volt[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
t += 0.5*dT;
Signal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
}
SetNyquistNum( CalcNyquistNum(f0+fc) );
}
void Excitation::CalcDiracPulsExcitation()
{
if (dT==0) return;
Length = 1;
// cerr << "Operator::CalcDiracPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0;
Signal_volt[1]=1.0;
Signal_curr[0]=0.0;
Signal_curr[1]=1.0;
SetNyquistNum( 1 );
}
void Excitation::CalcStepExcitation()
{
if (dT==0) return;
Length = 1;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=1.0;
Signal_volt[1]=1.0;
Signal_curr[0]=1.0;
Signal_curr[1]=1.0;
SetNyquistNum( 1 );
}
void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
{
if (dT==0) return;
if (nTS<=0) return;
Length = (unsigned int)(nTS);
// cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0;
Signal_curr[0]=0.0;
FunctionParser fParse;
fParse.AddConstant("pi", 3.14159265358979323846);
fParse.AddConstant("e", 2.71828182845904523536);
fParse.Parse(signal,"t");
if (fParse.GetParseErrorType()!=FunctionParser::FP_NO_ERROR)
{
cerr << "Operator::CalcCustomExcitation: Function Parser error: " << fParse.ErrorMsg() << endl;
exit(1);
}
double vars[1];
for (unsigned int n=1;n<Length+1;++n)
{
vars[0] = (n-1)*dT;
Signal_volt[n] = fParse.Eval(vars);
vars[0] += 0.5*dT;
Signal_curr[n] = fParse.Eval(vars);
}
SetNyquistNum( CalcNyquistNum(f0) );
}
void Excitation::CalcSinusExcitation(double f0, int nTS)
{
if (dT==0) return;
if (nTS<=0) return;
Length = (unsigned int)(nTS);
// cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0;
Signal_curr[0]=0.0;
for (unsigned int n=1;n<Length+1;++n)
{
double t = (n-1)*dT;
Signal_volt[n] = sin(2.0*PI*f0*t);
t += 0.5*dT;
Signal_curr[n] = sin(2.0*PI*f0*t);
}
SetNyquistNum( CalcNyquistNum(f0) );
}
void Excitation::setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit,
vector<unsigned int> const& volt_vDelay, vector<unsigned int> const& volt_vDir )
{
E_Count = volt_vIndex[0].size();
for (int n=0; n<3; n++) {
delete[] E_index[n];
E_index[n] = new unsigned int[E_Count];
}
delete[] E_delay;
delete[] E_amp;
delete[] E_dir;
E_delay = new unsigned int[E_Count];
E_amp = new FDTD_FLOAT[E_Count];
E_dir = new unsigned short[E_Count];
cerr << "Excitation::setupVoltageExcitation(): Number of voltage excitation points: " << E_Count << endl;
// if (E_Count==0)
// cerr << "No E-Field/voltage excitation found!" << endl;
for (int n=0; n<3; n++)
for (unsigned int i=0; i<E_Count; i++)
E_index[n][i] = volt_vIndex[n].at(i);
for (unsigned int i=0; i<E_Count; i++) {
E_delay[i] = volt_vDelay.at(i);
E_amp[i] = volt_vExcit.at(i);
E_dir[i] = volt_vDir.at(i);
}
}
void Excitation::setupCurrentExcitation( vector<unsigned int> const curr_vIndex[3], vector<FDTD_FLOAT> const& curr_vExcit,
vector<unsigned int> const& curr_vDelay, vector<unsigned int> const& curr_vDir )
{
Curr_Count = curr_vIndex[0].size();
for (int n=0; n<3; n++) {
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<Curr_Count; i++)
Curr_index[n][i] = curr_vIndex[n].at(i);
for (unsigned int i=0; i<Curr_Count; i++) {
Curr_delay[i] = curr_vDelay.at(i);
Curr_amp[i] = curr_vExcit.at(i);
Curr_dir[i] = curr_vDir.at(i);
}
}

82
FDTD/excitation.h Normal file
View File

@ -0,0 +1,82 @@
/*
* Copyright (C) 2010 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 EXCITATION_H
#define EXCITATION_H
#include <vector>
#include <string>
#include "tools/constants.h"
class TiXmlElement;
class Excitation
{
public:
Excitation( double timestep );
virtual ~Excitation();
bool setupExcitation( TiXmlElement* Excite, unsigned int maxTS );
//! Get the excitation timestep with the (first) max amplitude
virtual unsigned int GetMaxExcitationTimestep() const;
unsigned int CalcNyquistNum(double fmax) const;
void setupVoltageExcitation( vector<unsigned int> const volt_vIndex[3], vector<FDTD_FLOAT> const& volt_vExcit,
vector<unsigned int> const& volt_vDelay, vector<unsigned int> const& volt_vDir );
void setupCurrentExcitation( vector<unsigned int> const curr_vIndex[3], vector<FDTD_FLOAT> const& curr_vExcit,
vector<unsigned int> const& curr_vDelay, vector<unsigned int> const& curr_vDir );
void SetNyquistNum(unsigned int nyquist) {m_nyquistTS=nyquist;}
unsigned int GetNyquistNum() const {return m_nyquistTS;}
//Excitation time-signal
unsigned int Length;
FDTD_FLOAT* Signal_volt;
FDTD_FLOAT* Signal_curr;
//E-Field/voltage Excitation
unsigned int E_Count;
unsigned int* E_index[3];
unsigned short* E_dir;
FDTD_FLOAT* E_amp; //represented as edge-voltages!!
unsigned int* E_delay;
//H-Field/current Excitation
unsigned int Curr_Count;
unsigned int* Curr_index[3];
unsigned short* Curr_dir;
FDTD_FLOAT* Curr_amp; //represented as edge-currents!!
unsigned int* Curr_delay;
protected:
double dT;
unsigned int m_nyquistTS;
//! Calculate a custom signal
virtual void CalcCustomExcitation(double f0, int nTS, string signal);
//! Calculate an excitation with center of f0 and the half bandwidth fc
virtual void CalcGaussianPulsExcitation(double f0, double fc);
//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps
virtual void CalcSinusExcitation(double f0, int nTS);
//! Calculate a dirac impuls excitation
virtual void CalcDiracPulsExcitation();
//! Calculate a step excitation
virtual void CalcStepExcitation();
};
#endif // EXCITATION_H

243
FDTD/operator.cpp
View File

@ -22,6 +22,7 @@
#include "tools/array_ops.h"
#include "fparser.hh"
Operator* Operator::New()
{
Operator* op = new Operator();
@ -31,6 +32,7 @@ Operator* Operator::New()
Operator::Operator()
{
Exc = 0;
}
Operator::~Operator()
@ -45,24 +47,12 @@ void Operator::Init()
{
CSX = NULL;
ExciteSignal_volt = NULL;
ExciteSignal_curr = NULL;
E_Exc_delay = NULL;
E_Exc_amp=NULL;
E_Exc_dir=NULL;
vv=NULL;
vi=NULL;
Curr_Exc_delay = NULL;
Curr_Exc_amp=NULL;
Curr_Exc_dir=NULL;
iv=NULL;
ii=NULL;
for (int n=0;n<3;++n)
{
discLines[n]=NULL;
E_Exc_index[n]=NULL;
Curr_Exc_index[n]=NULL;
}
MainOp=NULL;
DualOp=NULL;
@ -77,28 +67,18 @@ void Operator::Init()
for (int n=0;n<6;++n)
m_BC[n]=0;
Exc = 0;
}
void Operator::Reset()
{
delete[] ExciteSignal_volt;
delete[] ExciteSignal_curr;
delete[] E_Exc_delay;
delete[] E_Exc_dir;
delete[] E_Exc_amp;
delete[] Curr_Exc_delay;
delete[] Curr_Exc_dir;
delete[] Curr_Exc_amp;
Delete_N_3DArray(vv,numLines);
Delete_N_3DArray(vi,numLines);
Delete_N_3DArray(iv,numLines);
Delete_N_3DArray(ii,numLines);
for (int n=0;n<3;++n)
{
delete[] discLines[n];
delete[] E_Exc_index[n];
delete[] Curr_Exc_index[n];
}
delete MainOp;
delete DualOp;
for (int n=0;n<3;++n)
@ -109,14 +89,9 @@ void Operator::Reset()
delete[] EC_R[n];
}
Init();
}
delete Exc;
unsigned int Operator::CalcNyquistNum(double fmax)
{
if (dT==0) return 1;
double T0 = 1/fmax;
return floor(T0/2/dT);
Init();
}
string Operator::GetDirName(int ny) const
@ -310,140 +285,13 @@ void Operator::ShowStat() const
cout << "Size of Field-Data : " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MB) " << endl;
cout << "-----------------------------------" << endl;
cout << "Timestep (s) : " << dT << endl;
cout << "Nyquist criteria (TS) : " << m_nyquistTS << endl;
cout << "Nyquist criteria (s) : " << m_nyquistTS*dT << endl;
cout << "Excitation Length (TS) : " << ExciteLength << endl;
cout << "Excitation Length (s) : " << ExciteLength*dT << endl;
cout << "Nyquist criteria (TS) : " << Exc->GetNyquistNum() << endl;
cout << "Nyquist criteria (s) : " << Exc->GetNyquistNum()*dT << endl;
cout << "Excitation Length (TS) : " << Exc->Length << endl;
cout << "Excitation Length (s) : " << Exc->Length*dT << endl;
cout << "-----------------------------------" << endl;
}
unsigned int Operator::GetMaxExcitationTimestep() const
{
FDTD_FLOAT maxAmp=0;
unsigned int maxStep=0;
for (unsigned int n=1;n<ExciteLength+1;++n)
{
if (fabs(ExciteSignal_volt[n])>maxAmp)
{
maxAmp = fabs(ExciteSignal_volt[n]);
maxStep = n;
}
}
return maxStep;
}
unsigned int Operator::CalcGaussianPulsExcitation(double f0, double fc)
{
if (dT==0) return 0;
ExciteLength = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
// cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] ExciteSignal_volt;
delete[] ExciteSignal_curr;
ExciteSignal_volt = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_curr = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_volt[0]=0.0;
ExciteSignal_curr[0]=0.0;
for (unsigned int n=1;n<ExciteLength+1;++n)
{
double t = (n-1)*dT;
ExciteSignal_volt[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
t += 0.5*dT;
ExciteSignal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
}
return CalcNyquistNum(f0+fc);
}
unsigned int Operator::CalcDiracPulsExcitation()
{
if (dT==0) return 0;
ExciteLength = 1;
// cerr << "Operator::CalcDiracPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] ExciteSignal_volt;
delete[] ExciteSignal_curr;
ExciteSignal_volt = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_curr = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_volt[0]=0.0;
ExciteSignal_volt[1]=1.0;
ExciteSignal_curr[0]=0.0;
ExciteSignal_curr[1]=1.0;
return 1;
}
unsigned int Operator::CalcStepExcitation()
{
if (dT==0) return 0;
ExciteLength = 1;
delete[] ExciteSignal_volt;
delete[] ExciteSignal_curr;
ExciteSignal_volt = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_curr = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_volt[0]=1.0;
ExciteSignal_volt[1]=1.0;
ExciteSignal_curr[0]=0.0;
ExciteSignal_curr[1]=1.0;
return 1;
}
unsigned int Operator::CalcCustomExcitation(double f0, int nTS, string signal)
{
if (dT==0) return 0;
if (nTS<=0) return 0;
ExciteLength = (unsigned int)(nTS);
// cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] ExciteSignal_volt;
delete[] ExciteSignal_curr;
ExciteSignal_volt = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_curr = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_volt[0]=0.0;
ExciteSignal_curr[0]=0.0;
FunctionParser fParse;
fParse.AddConstant("pi", 3.14159265358979323846);
fParse.AddConstant("e", 2.71828182845904523536);
fParse.Parse(signal,"t");
if (fParse.GetParseErrorType()!=FunctionParser::FP_NO_ERROR)
{
cerr << "Operator::CalcCustomExcitation: Function Parser error: " << fParse.ErrorMsg() << endl;
exit(1);
}
double vars[1];
for (unsigned int n=1;n<ExciteLength+1;++n)
{
vars[0] = (n-1)*GetTimestep();
ExciteSignal_volt[n] = fParse.Eval(vars);
vars[0] += 0.5*GetTimestep();
ExciteSignal_curr[n] = fParse.Eval(vars);
}
return CalcNyquistNum(f0);
}
unsigned int Operator::CalcSinusExcitation(double f0, int nTS)
{
if (dT==0) return 0;
if (nTS<=0) return 0;
ExciteLength = (unsigned int)(nTS);
// cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] ExciteSignal_volt;
delete[] ExciteSignal_curr;
ExciteSignal_volt = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_curr = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal_volt[0]=0.0;
ExciteSignal_curr[0]=0.0;
for (unsigned int n=1;n<ExciteLength+1;++n)
{
double t = (n-1)*dT;
ExciteSignal_volt[n] = sin(2.0*PI*f0*t);
t += 0.5*dT;
ExciteSignal_curr[n] = sin(2.0*PI*f0*t);
}
return CalcNyquistNum(f0);
}
void Operator::DumpOperator2File(string filename)
{
@ -455,9 +303,9 @@ void Operator::DumpOperator2File(string filename)
}
FDTD_FLOAT**** exc = Create_N_3DArray(numLines);
for (unsigned int n=0;n<E_Exc_Count;++n)
{
exc[E_Exc_dir[n]][E_Exc_index[0][n]][E_Exc_index[1][n]][E_Exc_index[2][n]] = E_Exc_amp[n];
if (Exc) {
for (unsigned int n=0;n<Exc->E_Count;++n)
exc[Exc->E_dir[n]][Exc->E_index[0][n]][Exc->E_index[1][n]][Exc->E_index[2][n]] = Exc->E_amp[n];
}
string names[] = {"vv", "vi", "iv" , "ii", "exc"};
@ -557,6 +405,12 @@ void Operator::InitOperator()
ii = Create_N_3DArray(numLines);
}
void Operator::InitExcitation()
{
delete Exc;
Exc = new Excitation( dT );
}
void Operator::Calc_ECOperatorPos(int n, unsigned int* pos)
{
unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]);
@ -576,6 +430,8 @@ int Operator::CalcECOperator()
InitOperator();
InitExcitation();
unsigned int pos[3];
for (int n=0;n<3;++n)
@ -925,7 +781,11 @@ double Operator::CalcTimestep()
bool Operator::CalcFieldExcitation()
{
if (dT==0) return false;
if (dT==0)
return false;
if (Exc==0)
return false;
unsigned int pos[3];
double delta[3];
double amp=0;
@ -1094,59 +954,16 @@ bool Operator::CalcFieldExcitation()
// set voltage excitations
E_Exc_Count = volt_vExcit.size();
cerr << "Operator::CalcFieldExcitation: Number of voltage excitation points: " << E_Exc_Count << endl;
if (E_Exc_Count==0)
cerr << "No E-Field/voltage excitation found!" << endl;
for (int n=0;n<3;++n)
{
delete[] E_Exc_index[n];
E_Exc_index[n] = new unsigned int[E_Exc_Count];
for (unsigned int i=0;i<E_Exc_Count;++i)
E_Exc_index[n][i]=volt_vIndex[n].at(i);
}
delete[] E_Exc_delay;
E_Exc_delay = new unsigned int[E_Exc_Count];
delete[] E_Exc_amp;
E_Exc_amp = new FDTD_FLOAT[E_Exc_Count];
delete[] E_Exc_dir;
E_Exc_dir = new unsigned short[E_Exc_Count];
for (unsigned int i=0;i<E_Exc_Count;++i)
{
E_Exc_delay[i]=volt_vDelay.at(i);
E_Exc_amp[i]=volt_vExcit.at(i);
E_Exc_dir[i]=volt_vDir.at(i);
}
Exc->setupVoltageExcitation( volt_vIndex, volt_vExcit, volt_vDelay, volt_vDir );
// set current excitations
Curr_Exc_Count = curr_vExcit.size();
cerr << "Operator::CalcFieldExcitation: Number of current excitation points: " << Curr_Exc_Count << endl;
if (Curr_Exc_Count==0)
cerr << "No H-Field/current excitation found!" << endl;
for (int n=0;n<3;++n)
{
delete[] Curr_Exc_index[n];
Curr_Exc_index[n] = new unsigned int[Curr_Exc_Count];
for (unsigned int i=0;i<Curr_Exc_Count;++i)
Curr_Exc_index[n][i]=curr_vIndex[n].at(i);
}
delete[] Curr_Exc_delay;
Curr_Exc_delay = new unsigned int[Curr_Exc_Count];
delete[] Curr_Exc_amp;
Curr_Exc_amp = new FDTD_FLOAT[Curr_Exc_Count];
delete[] Curr_Exc_dir;
Curr_Exc_dir = new unsigned short[Curr_Exc_Count];
for (unsigned int i=0;i<Curr_Exc_Count;++i)
{
Curr_Exc_delay[i]=curr_vDelay.at(i);
Curr_Exc_amp[i]=curr_vExcit.at(i);
Curr_Exc_dir[i]=curr_vDir.at(i);
}
Exc->setupCurrentExcitation( curr_vIndex, curr_vExcit, curr_vDelay, curr_vDir );
return true;
}
bool Operator::CalcPEC()
{
unsigned int pos[3];

43
FDTD/operator.h
View File

@ -21,8 +21,7 @@
#include "ContinuousStructure.h"
#include "tools/AdrOp.h"
#include "tools/constants.h"
#define FDTD_FLOAT float
#include "excitation.h"
class Operator_Extension;
@ -45,21 +44,6 @@ public:
inline virtual FDTD_FLOAT& GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[n][x][y][z]; }
inline virtual FDTD_FLOAT& GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[n][x][y][z]; }
//! Calculate a custom signal \return number of Nyquist timesteps defined by f0
virtual unsigned int CalcCustomExcitation(double f0, int nTS, string signal);
//! Calculate an excitation with center of f0 and the half bandwidth fc \return number of Nyquist timesteps
virtual unsigned int CalcGaussianPulsExcitation(double f0, double fc);
//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps \return number of Nyquist timesteps
virtual unsigned int CalcSinusExcitation(double f0, int nTS);
//! Calculate a dirac impuls excitation \return number of Nyquist timesteps
virtual unsigned int CalcDiracPulsExcitation();
//! Calculate a step excitation \return number of Nyquist timesteps
virtual unsigned int CalcStepExcitation();
//! Get the excitation timestep with the (first) max amplitude
virtual unsigned int GetMaxExcitationTimestep() const;
virtual void SetBoundaryCondition(int* BCs) {for (int n=0;n<6;++n) m_BC[n]=BCs[n];}
virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
virtual void ApplyMagneticBC(bool* dirs);
@ -69,10 +53,6 @@ public:
virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];}
void SetNyquistNum(unsigned int nyquist) {m_nyquistTS=nyquist;}
unsigned int GetNyquistNum() const {return m_nyquistTS;}
unsigned int CalcNyquistNum(double fmax);
void ShowStat() const;
void DumpOperator2File(string filename);
@ -102,6 +82,7 @@ protected:
virtual void Init();
virtual void Reset();
virtual void InitOperator();
virtual void InitExcitation();
struct Grid_Path
{
@ -122,7 +103,6 @@ protected:
//Calc timestep only internal use
virtual double CalcTimestep();
double dT; //FDTD timestep!
unsigned int m_nyquistTS;
//! Calc operator at certain pos
virtual void Calc_ECOperatorPos(int n, unsigned int* pos);
@ -152,24 +132,7 @@ public:
FDTD_FLOAT**** ii; //calc new current from old current
FDTD_FLOAT**** iv; //calc new current from old voltage
//Excitation time-signal
unsigned int ExciteLength;
FDTD_FLOAT* ExciteSignal_volt;
FDTD_FLOAT* ExciteSignal_curr;
//E-Field/voltage Excitation
unsigned int E_Exc_Count;
unsigned int* E_Exc_index[3];
unsigned short* E_Exc_dir;
FDTD_FLOAT* E_Exc_amp; //represented as edge-voltages!!
unsigned int* E_Exc_delay;
//H-Field/current Excitation
unsigned int Curr_Exc_Count;
unsigned int* Curr_Exc_index[3];
unsigned short* Curr_Exc_dir;
FDTD_FLOAT* Curr_Exc_amp; //represented as edge-currents!!
unsigned int* Curr_Exc_delay;
Excitation* Exc;
};
#endif // OPERATOR_H

5
FDTD/operator_sse.cpp
View File

@ -49,8 +49,11 @@ void Operator_sse::Reset()
Delete_N_3DArray_v4sf(f4_vi,numLines);
Delete_N_3DArray_v4sf(f4_iv,numLines);
Delete_N_3DArray_v4sf(f4_ii,numLines);
f4_vv = 0;
f4_vi = 0;
f4_iv = 0;
f4_ii = 0;
Operator::Reset();
// Init(); // FIXME this calls Operator::Init() twice...
}
void Operator_sse::InitOperator()

113
openEMS.pro
View File

@ -46,12 +46,13 @@ SOURCES += main.cpp \
FDTD/engine_multithread.cpp \
FDTD/operator_cylinder.cpp \
FDTD/engine_cylinder.cpp \
FDTD/engine_sse.cpp \
FDTD/operator_sse.cpp \
FDTD/engine_sse.cpp \
FDTD/operator_sse.cpp \
FDTD/operator_extension.cpp \
FDTD/engine_extension.cpp \
FDTD/engine_ext_mur_abc.cpp \
FDTD/operator_ext_mur_abc.cpp
FDTD/operator_ext_mur_abc.cpp \
FDTD/excitation.cpp
HEADERS += tools/ErrorMsg.h \
tools/AdrOp.h \
tools/constants.h \
@ -67,71 +68,73 @@ HEADERS += tools/ErrorMsg.h \
openems.h \
FDTD/engine_multithread.h \
FDTD/operator_cylinder.h \
FDTD/engine_sse.h \
FDTD/operator_sse.h \
FDTD/engine_sse.h \
FDTD/operator_sse.h \
FDTD/engine_cylinder.h \
FDTD/operator_extension.h \
FDTD/engine_extension.h \
FDTD/engine_ext_mur_abc.h \
FDTD/operator_ext_mur_abc.h
FDTD/operator_ext_mur_abc.h \
FDTD/excitation.h
QMAKE_CXXFLAGS_RELEASE = -O2 \
-g \
-march=native
-march=native -msse -msse2
QMAKE_CXXFLAGS_DEBUG = -O0 \
-g \
-march=native
-march=native -msse -msse2
#
# to use ABI2 target:
# qmake CONFIG+="ABI2 bits64" -o Makefile.ABI2-64 openEMS.pro
# make -fMakefile.ABI2-64
#
ABI2 {
CONFIG-=debug debug_and_release
CONFIG+=release
QMAKE_CFLAGS_RELEASE=-O2 -fabi-version=2
QMAKE_CXXFLAGS_RELEASE=-O2 -fabi-version=2
QMAKE_CC = apgcc
QMAKE_CXX = apg++
QMAKE_LINK = apg++
QMAKE_LINK_SHLIB = apg++
QMAKE_LFLAGS_RPATH =
QMAKE_LFLAGS = \'-Wl,-rpath,\$$ORIGIN/lib\'
ABI2 {
CONFIG -= debug \
debug_and_release
CONFIG += release
QMAKE_CFLAGS_RELEASE = -O2 \
-fabi-version=2
QMAKE_CXXFLAGS_RELEASE = -O2 \
-fabi-version=2
QMAKE_CC = apgcc
QMAKE_CXX = apg++
QMAKE_LINK = apg++
QMAKE_LINK_SHLIB = apg++
QMAKE_LFLAGS_RPATH =
QMAKE_LFLAGS = \'-Wl,-rpath,\$$ORIGIN/lib\'
}
bits64 {
QMAKE_CXXFLAGS_RELEASE+=-m64 -march=athlon64
QMAKE_LFLAGS_RELEASE+=-m64 -march=athlon64
OBJECTS_DIR = ABI2-64
LIBS = ../CSXCAD/ABI2-64/libCSXCAD.so
LIBS += ../fparser/ABI2-64/libfparser.so
LIBS += ../tinyxml/ABI2-64/libtinyxml.so
LIBS += ../boost-64/lib/libboost_thread.so
LIBS += ../hdf5-64/lib/libhdf5.so
LIBS += ../hdf5-64/lib/libhdf5_cpp.so -lpthread
INCLUDEPATH += ../hdf5-64/include
INCLUDEPATH += ../boost-64/include
bits64 {
QMAKE_CXXFLAGS_RELEASE += -m64 \
-march=athlon64
QMAKE_LFLAGS_RELEASE += -m64 \
-march=athlon64
OBJECTS_DIR = ABI2-64
LIBS = ../CSXCAD/ABI2-64/libCSXCAD.so
LIBS += ../fparser/ABI2-64/libfparser.so
LIBS += ../tinyxml/ABI2-64/libtinyxml.so
LIBS += ../boost-64/lib/libboost_thread.so
LIBS += ../hdf5-64/lib/libhdf5.so
LIBS += ../hdf5-64/lib/libhdf5_cpp.so \
-lpthread
INCLUDEPATH += ../hdf5-64/include
INCLUDEPATH += ../boost-64/include
}
bits32 {
QMAKE_CXXFLAGS_RELEASE+=-m32 -march=i686
QMAKE_LFLAGS_RELEASE+=-m32 -march=i686
OBJECTS_DIR = ABI2-32
LIBS = ../CSXCAD/ABI2-32/libCSXCAD.so
LIBS += ../fparser/ABI2-32/libfparser.so
LIBS += ../tinyxml/ABI2-32/libtinyxml.so
LIBS += ../boost-32/lib/libboost_thread.so
LIBS += ../hdf5-32/lib/libhdf5.so
LIBS += ../hdf5-32/lib/libhdf5_cpp.so
INCLUDEPATH += ../hdf5-32/include
INCLUDEPATH += ../boost-32/include
bits32 {
QMAKE_CXXFLAGS_RELEASE += -m32 \
-march=i686
QMAKE_LFLAGS_RELEASE += -m32 \
-march=i686
OBJECTS_DIR = ABI2-32
LIBS = ../CSXCAD/ABI2-32/libCSXCAD.so
LIBS += ../fparser/ABI2-32/libfparser.so
LIBS += ../tinyxml/ABI2-32/libtinyxml.so
LIBS += ../boost-32/lib/libboost_thread.so
LIBS += ../hdf5-32/lib/libhdf5.so
LIBS += ../hdf5-32/lib/libhdf5_cpp.so
INCLUDEPATH += ../hdf5-32/include
INCLUDEPATH += ../boost-32/include
}
ABI2 {
DESTDIR = $$OBJECTS_DIR
MOC_DIR = $$OBJECTS_DIR
UI_DIR = $$OBJECTS_DIR
RCC_DIR = $$OBJECTS_DIR
ABI2 {
DESTDIR = $$OBJECTS_DIR
MOC_DIR = $$OBJECTS_DIR
UI_DIR = $$OBJECTS_DIR
RCC_DIR = $$OBJECTS_DIR
}

56
openems.cpp
View File

@ -124,51 +124,6 @@ bool openEMS::parseCommandLineArgument( const char *argv )
return false;
}
void openEMS::SetupExcitation(TiXmlElement* Excite)
{
if (Excite==NULL)
{
cerr << "Can't read openEMS Excitation Settings... " << endl;
exit(-2);
}
int Excit_Type=0;
double f0=0;
double fc=0;
Excite->QueryIntAttribute("Type",&Excit_Type);
unsigned int Nyquist = 0;
switch (Excit_Type)
{
case 0:
Excite->QueryDoubleAttribute("f0",&f0);
Excite->QueryDoubleAttribute("fc",&fc);
Nyquist = FDTD_Op->CalcGaussianPulsExcitation(f0,fc);
break;
case 1:
Excite->QueryDoubleAttribute("f0",&f0);
Nyquist = FDTD_Op->CalcSinusExcitation(f0,NrTS);
break;
case 2:
Nyquist = FDTD_Op->CalcDiracPulsExcitation();
break;
case 3:
Nyquist = FDTD_Op->CalcStepExcitation();
break;
case 10:
Excite->QueryDoubleAttribute("f0",&f0);
Nyquist = FDTD_Op->CalcCustomExcitation(f0,NrTS,Excite->Attribute("Function"));
break;
}
if (!Nyquist)
{
cerr << "openEMS: excitation setup failed!!" << endl;
exit(2);
}
FDTD_Op->SetNyquistNum(Nyquist);
}
int openEMS::SetupFDTD(const char* file)
{
if (file==NULL) return -1;
@ -278,7 +233,8 @@ int openEMS::SetupFDTD(const char* file)
FDTD_Op->CalcECOperator();
SetupExcitation(FDTD_Opts->FirstChildElement("Excitation"));
if (!FDTD_Op->Exc->setupExcitation( FDTD_Opts->FirstChildElement("Excitation"), NrTS ))
exit(2);
if (DebugMat)
{
@ -318,7 +274,7 @@ int openEMS::SetupFDTD(const char* file)
//*************** setup processing ************//
cout << "Setting up processing..." << endl;
unsigned int Nyquist = FDTD_Op->GetNyquistNum();
unsigned int Nyquist = FDTD_Op->Exc->GetNyquistNum();
PA = new ProcessingArray(Nyquist);
double start[3];
@ -414,9 +370,9 @@ void openEMS::RunFDTD()
double maxE=0,currE=0;
//add all timesteps to end-crit field processing with max excite amplitude
unsigned int maxExcite = FDTD_Op->GetMaxExcitationTimestep();
for (unsigned int n=0;n<FDTD_Op->E_Exc_Count;++n)
ProcField->AddStep(FDTD_Op->E_Exc_delay[n]+maxExcite);
unsigned int maxExcite = FDTD_Op->Exc->GetMaxExcitationTimestep();
for (unsigned int n=0;n<FDTD_Op->Exc->E_Count;++n)
ProcField->AddStep(FDTD_Op->Exc->E_delay[n]+maxExcite);
double change=1;
int prevTS=0,currTS=0;

2
openems.h
View File

@ -45,8 +45,6 @@ public:
void DebugBox() {m_debugBox=true;}
protected:
void SetupExcitation(TiXmlElement* Excite);
bool CylinderCoords;
//! Number of Timesteps

2
tools/constants.h
View File

@ -18,6 +18,8 @@
#ifndef CONSTANTS_H
#define CONSTANTS_H
#define FDTD_FLOAT float
#define __EPS0__ 8.85418781762e-12
#define __MUE0__ 1.256637062e-6
#define __C0__ 299792458