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

Init style: applied astyle.sh 

astyle --style=allman --indent=tab --indent-cases {} \;

please apply this style prior to any future commit

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
pull/1/head
Thorsten Liebig 2010-12-06 13:04:37 +01:00
parent 70ae86d63f
commit 9d3761b1d0
57 changed files with 886 additions and 717 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -7,4 +7,5 @@ Makefile*
*~ *~
*.so* *.so*
*.pro.user* *.pro.user*
*.orig
openEMS openEMS

8
Common/engine_interface_base.h
View File

@ -40,18 +40,18 @@ public:
InterpolationType GetInterpolationType() {return m_InterpolType;} InterpolationType GetInterpolationType() {return m_InterpolType;}
//! Get the (interpolated) electric field at \p pos. \sa SetInterpolationType //! Get the (interpolated) electric field at \p pos. \sa SetInterpolationType
virtual double* GetEField(const unsigned int* pos, double* out) const {UNUSED(pos);return out;} virtual double* GetEField(const unsigned int* pos, double* out) const {UNUSED(pos); return out;}
//! Get the (interpolated) magnetic field at \p pos. \sa SetInterpolationType //! Get the (interpolated) magnetic field at \p pos. \sa SetInterpolationType
virtual double* GetHField(const unsigned int* pos, double* out) const {UNUSED(pos);return out;} virtual double* GetHField(const unsigned int* pos, double* out) const {UNUSED(pos); return out;}
//! Calculate the electric field integral along a given line //! Calculate the electric field integral along a given line
virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const {UNUSED(start);UNUSED(stop);return 0.0;} virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const {UNUSED(start); UNUSED(stop); return 0.0;}
//! Convert the interpolation type into a string. //! Convert the interpolation type into a string.
static std::string GetInterpolationNameByType(InterpolationType mode); static std::string GetInterpolationNameByType(InterpolationType mode);
//! Get the current simulation time //! Get the current simulation time
virtual double GetTime(bool dualTime=false) const {UNUSED(dualTime);return 0.0;} virtual double GetTime(bool dualTime=false) const {UNUSED(dualTime); return 0.0;}
//! Get the current number of timesteps //! Get the current number of timesteps
virtual unsigned int GetNumberOfTimesteps() const {return 0;} virtual unsigned int GetNumberOfTimesteps() const {return 0;}

6
Common/operator_base.cpp
View File

@ -39,15 +39,15 @@ std::string Operator_Base::GetDirName(int ny) const
void Operator_Base::Init() void Operator_Base::Init()
{ {
dT = 0; dT = 0;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
discLines[n]=NULL; discLines[n]=NULL;
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
m_BC[n]=0; m_BC[n]=0;
} }
void Operator_Base::Reset() void Operator_Base::Reset()
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
delete[] discLines[n]; delete[] discLines[n];
Init(); Init();
} }

16
Common/operator_base.h
View File

@ -44,32 +44,32 @@ public:
virtual double GetGridDelta() const {return 0;} virtual double GetGridDelta() const {return 0;}
//! Get the mesh delta times the grid delta for a 3D position (unit is meter) //! Get the mesh delta times the grid delta for a 3D position (unit is meter)
virtual double GetMeshDelta(int n, const unsigned int* pos, bool dualMesh=false) const {UNUSED(n);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetMeshDelta(int n, const unsigned int* pos, bool dualMesh=false) const {UNUSED(n); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
//! Get the disc line in \a n direction (in drawing units) //! Get the disc line in \a n direction (in drawing units)
virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const {UNUSED(n);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetDiscLine(int n, unsigned int pos, bool dualMesh=false) const {UNUSED(n); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
//! Get the node width for a given direction \a n and a given mesh position \a pos //! Get the node width for a given direction \a n and a given mesh position \a pos
virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
//! Get the node area for a given direction \a n and a given mesh position \a pos //! Get the node area for a given direction \a n and a given mesh position \a pos
virtual double GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
//! Get the length of an FDTD edge (unit is meter). //! Get the length of an FDTD edge (unit is meter).
virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos //! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
/*! /*!
This will return the area around an edge with a given direction, measured at the middle of the edge. This will return the area around an edge with a given direction, measured at the middle of the edge.
In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems. In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems.
*/ */
virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetEdgeArea(int ny, const unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
//! Snap the given coodinates to mesh indices //! Snap the given coodinates to mesh indices
virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL) {UNUSED(coord);UNUSED(uicoord);UNUSED(lower);UNUSED(inside);return false;}; virtual bool SnapToMesh(double* coord, unsigned int* uicoord, bool lower=false, bool* inside=NULL) {UNUSED(coord); UNUSED(uicoord); UNUSED(lower); UNUSED(inside); return false;};
//! Set the boundary conditions //! Set the boundary conditions
virtual void SetBoundaryCondition(int* BCs) {for (int n=0;n<6;++n) m_BC[n]=BCs[n];} virtual void SetBoundaryCondition(int* BCs) {for (int n=0; n<6; ++n) m_BC[n]=BCs[n];}
protected: protected:
Operator_Base(); Operator_Base();

4
Common/process_efield.cpp
View File

@ -65,7 +65,7 @@ void ProcessEField::Dump_FD_Data(vector<double_complex> value[3], double factor,
time_t rawTime; time_t rawTime;
time(&rawTime); time(&rawTime);
file << "%dump by openEMS @" << ctime(&rawTime) << "%frequency\treal_x\timag_x\treal_y\timag_y\treal_z\timag_z\n"; file << "%dump by openEMS @" << ctime(&rawTime) << "%frequency\treal_x\timag_x\treal_y\timag_y\treal_z\timag_z\n";
for (size_t n=0;n<value[0].size();++n) for (size_t n=0; n<value[0].size(); ++n)
{ {
file << m_FD_Samples.at(n) file << m_FD_Samples.at(n)
<< "\t" << 2.0 * std::real(value[0].at(n))*factor << "\t" << 2.0 * std::imag(value[0].at(n))*factor << "\t" << 2.0 * std::real(value[0].at(n))*factor << "\t" << 2.0 * std::imag(value[0].at(n))*factor
@ -110,7 +110,7 @@ int ProcessEField::Process()
{ {
FDTD_FLOAT field = Eng->GetVolt(pol,start) / Op->GetMeshDelta(pol,start); FDTD_FLOAT field = Eng->GetVolt(pol,start) / Op->GetMeshDelta(pol,start);
field *= m_weight; field *= m_weight;
for (size_t n=0;n<m_FD_Samples.size();++n) for (size_t n=0; n<m_FD_Samples.size(); ++n)
{ {
FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T ); FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T );
} }

5
Common/process_hfield.cpp
View File

@ -47,7 +47,8 @@ void ProcessHField::DefineStartStopCoord(double* dstart, double* dstop)
if (Op->SnapToMesh(dstop,stop,true,m_stop_inside)==false) if (Op->SnapToMesh(dstop,stop,true,m_stop_inside)==false)
cerr << "ProcessHField::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl; cerr << "ProcessHField::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
if (g_settings.showProbeDiscretization()) { if (g_settings.showProbeDiscretization())
{
cerr << m_Name << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], true ) << "," cerr << m_Name << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], true ) << ","
<< Op->GetDiscLine( 1, start[1], true ) << "," << Op->GetDiscLine( 2, start[2], true ) << ") -> (" << Op->GetDiscLine( 1, start[1], true ) << "," << Op->GetDiscLine( 2, start[2], true ) << ") -> ("
<< Op->GetDiscLine( 0, stop[0], true ) << ","<< Op->GetDiscLine( 1, stop[1], true ) << "," << Op->GetDiscLine( 0, stop[0], true ) << ","<< Op->GetDiscLine( 1, stop[1], true ) << ","
@ -93,7 +94,7 @@ int ProcessHField::Process()
{ {
FDTD_FLOAT field = Eng->GetCurr(pol,start) / Op->GetMeshDelta(pol,start,true); FDTD_FLOAT field = Eng->GetCurr(pol,start) / Op->GetMeshDelta(pol,start,true);
field *= m_weight; field *= m_weight;
for (size_t n=0;n<m_FD_Samples.size();++n) for (size_t n=0; n<m_FD_Samples.size(); ++n)
{ {
FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T ); FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T );
} }

26
Common/processcurrent.cpp
View File

@ -43,7 +43,7 @@ double ProcessCurrent::CalcIntegral()
int Dump_Dim = 0; int Dump_Dim = 0;
int NormDir = 0; int NormDir = 0;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if (start[n]>stop[n]) if (start[n]>stop[n])
{ {
@ -75,55 +75,55 @@ double ProcessCurrent::CalcIntegral()
case 0: case 0:
//y-current //y-current
if (m_stop_inside[0] && m_start_inside[2]) if (m_stop_inside[0] && m_start_inside[2])
for (unsigned int i=start[1]+1;i<=stop[1];++i) for (unsigned int i=start[1]+1; i<=stop[1]; ++i)
current+=Eng->GetCurr(1,stop[0],i,start[2]); current+=Eng->GetCurr(1,stop[0],i,start[2]);
//z-current //z-current
if (m_stop_inside[0] && m_stop_inside[1]) if (m_stop_inside[0] && m_stop_inside[1])
for (unsigned int i=start[2]+1;i<=stop[2];++i) for (unsigned int i=start[2]+1; i<=stop[2]; ++i)
current+=Eng->GetCurr(2,stop[0],stop[1],i); current+=Eng->GetCurr(2,stop[0],stop[1],i);
//y-current //y-current
if (m_start_inside[0] && m_stop_inside[2]) if (m_start_inside[0] && m_stop_inside[2])
for (unsigned int i=start[1]+1;i<=stop[1];++i) for (unsigned int i=start[1]+1; i<=stop[1]; ++i)
current-=Eng->GetCurr(1,start[0],i,stop[2]); current-=Eng->GetCurr(1,start[0],i,stop[2]);
//z-current //z-current
if (m_start_inside[0] && m_start_inside[1]) if (m_start_inside[0] && m_start_inside[1])
for (unsigned int i=start[2]+1;i<=stop[2];++i) for (unsigned int i=start[2]+1; i<=stop[2]; ++i)
current-=Eng->GetCurr(2,start[0],start[1],i); current-=Eng->GetCurr(2,start[0],start[1],i);
break; break;
case 1: case 1:
//z-current //z-current
if (m_start_inside[0] && m_start_inside[1]) if (m_start_inside[0] && m_start_inside[1])
for (unsigned int i=start[2]+1;i<=stop[2];++i) for (unsigned int i=start[2]+1; i<=stop[2]; ++i)
current+=Eng->GetCurr(2,start[0],start[1],i); current+=Eng->GetCurr(2,start[0],start[1],i);
//x-current //x-current
if (m_stop_inside[1] && m_stop_inside[2]) if (m_stop_inside[1] && m_stop_inside[2])
for (unsigned int i=start[0]+1;i<=stop[0];++i) for (unsigned int i=start[0]+1; i<=stop[0]; ++i)
current+=Eng->GetCurr(0,i,stop[1],stop[2]); current+=Eng->GetCurr(0,i,stop[1],stop[2]);
//z-current //z-current
if (m_stop_inside[0] && m_stop_inside[1]) if (m_stop_inside[0] && m_stop_inside[1])
for (unsigned int i=start[2]+1;i<=stop[2];++i) for (unsigned int i=start[2]+1; i<=stop[2]; ++i)
current-=Eng->GetCurr(2,stop[0],stop[1],i); current-=Eng->GetCurr(2,stop[0],stop[1],i);
//x-current //x-current
if (m_start_inside[1] && m_start_inside[2]) if (m_start_inside[1] && m_start_inside[2])
for (unsigned int i=start[0]+1;i<=stop[0];++i) for (unsigned int i=start[0]+1; i<=stop[0]; ++i)
current-=Eng->GetCurr(0,i,start[1],start[2]); current-=Eng->GetCurr(0,i,start[1],start[2]);
break; break;
case 2: case 2:
//x-current //x-current
if (m_start_inside[1] && m_start_inside[2]) if (m_start_inside[1] && m_start_inside[2])
for (unsigned int i=start[0]+1;i<=stop[0];++i) for (unsigned int i=start[0]+1; i<=stop[0]; ++i)
current+=Eng->GetCurr(0,i,start[1],start[2]); current+=Eng->GetCurr(0,i,start[1],start[2]);
//y-current //y-current
if (m_stop_inside[0] && m_start_inside[2]) if (m_stop_inside[0] && m_start_inside[2])
for (unsigned int i=start[1]+1;i<=stop[1];++i) for (unsigned int i=start[1]+1; i<=stop[1]; ++i)
current+=Eng->GetCurr(1,stop[0],i,start[2]); current+=Eng->GetCurr(1,stop[0],i,start[2]);
//x-current //x-current
if (m_stop_inside[1] && m_stop_inside[2]) if (m_stop_inside[1] && m_stop_inside[2])
for (unsigned int i=start[0]+1;i<=stop[0];++i) for (unsigned int i=start[0]+1; i<=stop[0]; ++i)
current-=Eng->GetCurr(0,i,stop[1],stop[2]); current-=Eng->GetCurr(0,i,stop[1],stop[2]);
//y-current //y-current
if (m_start_inside[0] && m_stop_inside[2]) if (m_start_inside[0] && m_stop_inside[2])
for (unsigned int i=start[1]+1;i<=stop[1];++i) for (unsigned int i=start[1]+1; i<=stop[1]; ++i)
current-=Eng->GetCurr(1,start[0],i,stop[2]); current-=Eng->GetCurr(1,start[0],i,stop[2]);
break; break;
default: default:

71
Common/processfields.cpp
View File

@ -29,7 +29,7 @@ ProcessFields::ProcessFields(Operator_Base* op) : Processing(op)
SetSubSampling(1); SetSubSampling(1);
SetPrecision(6); SetPrecision(6);
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
numLines[n]=0; numLines[n]=0;
discLines[n]=NULL; discLines[n]=NULL;
@ -38,7 +38,7 @@ ProcessFields::ProcessFields(Operator_Base* op) : Processing(op)
ProcessFields::~ProcessFields() ProcessFields::~ProcessFields()
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delete[] discLines[n]; delete[] discLines[n];
discLines[n]=NULL; discLines[n]=NULL;
@ -57,7 +57,7 @@ void ProcessFields::InitProcess()
H5::H5File* file = new H5::H5File( m_filename , H5F_ACC_TRUNC ); H5::H5File* file = new H5::H5File( m_filename , H5F_ACC_TRUNC );
H5::Group* group = new H5::Group( file->createGroup( "/Mesh" )); H5::Group* group = new H5::Group( file->createGroup( "/Mesh" ));
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
hsize_t dimsf[1]; // dataset dimensions hsize_t dimsf[1]; // dataset dimensions
dimsf[0] = numLines[n]; dimsf[0] = numLines[n];
@ -66,7 +66,7 @@ void ProcessFields::InitProcess()
H5::DataSet dataset = group->createDataSet( names[n].c_str(), datatype, dataspace ); H5::DataSet dataset = group->createDataSet( names[n].c_str(), datatype, dataspace );
//convert to float... //convert to float...
float* array = new float[numLines[n]]; float* array = new float[numLines[n]];
for (unsigned int i=0;i<numLines[n];++i) for (unsigned int i=0; i<numLines[n]; ++i)
{ {
#ifdef OUTPUT_IN_DRAWINGUNITS #ifdef OUTPUT_IN_DRAWINGUNITS
array[i] = Lines[n][i]; array[i] = Lines[n][i];
@ -117,7 +117,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
if (!Op->SnapToMesh(dstop,stop,dualMesh)) if (!Op->SnapToMesh(dstop,stop,dualMesh))
cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl; cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
// normalize order of start and stop // normalize order of start and stop
if (start[n]>stop[n]) if (start[n]>stop[n])
@ -129,14 +129,14 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
// construct new discLines // construct new discLines
lines.clear(); lines.clear();
for (unsigned int i=start[n];i<=stop[n];i+=subSample[n]) for (unsigned int i=start[n]; i<=stop[n]; i+=subSample[n])
{ {
lines.push_back(Op->GetDiscLine(n,i,dualMesh)); lines.push_back(Op->GetDiscLine(n,i,dualMesh));
} }
numLines[n] = lines.size(); numLines[n] = lines.size();
delete[] discLines[n]; delete[] discLines[n];
discLines[n] = new double[numLines[n]]; discLines[n] = new double[numLines[n]];
for (unsigned int i=0;i<numLines[n];++i) for (unsigned int i=0; i<numLines[n]; ++i)
discLines[n][i] = lines.at(i); discLines[n][i] = lines.at(i);
} }
} }
@ -146,7 +146,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
if (Op->SnapToMesh(dstop,stop)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl; if (Op->SnapToMesh(dstop,stop)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
//create mesh //create mesh
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if (start[n]>stop[n]) if (start[n]>stop[n])
{ {
@ -158,14 +158,14 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
// --stop[n]; // --stop[n];
// cerr << "start " << start[n] << "stop " << stop[n]; // cerr << "start " << start[n] << "stop " << stop[n];
lines.clear(); lines.clear();
for (unsigned int i=start[n];i<=stop[n];i+=subSample[n]) for (unsigned int i=start[n]; i<=stop[n]; i+=subSample[n])
{ {
lines.push_back(Op->GetDiscLine(n,i));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i])); lines.push_back(Op->GetDiscLine(n,i));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i]));
} }
numLines[n] = lines.size(); numLines[n] = lines.size();
delete[] discLines[n]; delete[] discLines[n];
discLines[n] = new double[numLines[n]]; discLines[n] = new double[numLines[n]];
for (unsigned int i=0;i<numLines[n];++i) for (unsigned int i=0; i<numLines[n]; ++i)
discLines[n][i] = lines.at(i); discLines[n][i] = lines.at(i);
} }
} }
@ -175,7 +175,7 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl; if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl;
//create dual mesh //create dual mesh
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
// cerr << "start " << start[n] << "stop " << stop[n]; // cerr << "start " << start[n] << "stop " << stop[n];
if (start[n]>stop[n]) if (start[n]>stop[n])
@ -186,19 +186,20 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
} }
++stop[n]; ++stop[n];
lines.clear(); lines.clear();
for (unsigned int i=start[n];i<stop[n];i+=subSample[n]) for (unsigned int i=start[n]; i<stop[n]; i+=subSample[n])
{ {
lines.push_back(Op->GetDiscLine(n,i,true));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i])); lines.push_back(Op->GetDiscLine(n,i,true));//0.5*(Op->discLines[n][i+1] + Op->discLines[n][i]));
} }
numLines[n] = lines.size(); numLines[n] = lines.size();
delete[] discLines[n]; delete[] discLines[n];
discLines[n] = new double[numLines[n]]; discLines[n] = new double[numLines[n]];
for (unsigned int i=0;i<numLines[n];++i) for (unsigned int i=0; i<numLines[n]; ++i)
discLines[n][i] = lines.at(i); discLines[n][i] = lines.at(i);
} }
} }
if (g_settings.showProbeDiscretization()) { if (g_settings.showProbeDiscretization())
{
// FIXME the information E-Field / H-Field and therefore which mesh to use is missing // FIXME the information E-Field / H-Field and therefore which mesh to use is missing
bool dualMesh = false; bool dualMesh = false;
cerr << m_filename << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], dualMesh ) << "," cerr << m_filename << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], dualMesh ) << ","
@ -222,11 +223,11 @@ double ProcessFields::CalcTotalEnergy() const
const Engine* Eng = EI_FDTD->GetFDTDEngine(); const Engine* Eng = EI_FDTD->GetFDTDEngine();
unsigned int pos[3]; unsigned int pos[3];
for (pos[0]=0;pos[0]<Op->GetNumberOfLines(0);++pos[0]) for (pos[0]=0; pos[0]<Op->GetNumberOfLines(0); ++pos[0])
{ {
for (pos[1]=0;pos[1]<Op->GetNumberOfLines(1);++pos[1]) for (pos[1]=0; pos[1]<Op->GetNumberOfLines(1); ++pos[1])
{ {
for (pos[2]=0;pos[2]<Op->GetNumberOfLines(2);++pos[2]) for (pos[2]=0; pos[2]<Op->GetNumberOfLines(2); ++pos[2])
{ {
energy+=fabs(Eng->GetVolt(0,pos[0],pos[1],pos[2]) * Eng->GetCurr(1,pos[0],pos[1],pos[2])); energy+=fabs(Eng->GetVolt(0,pos[0],pos[1],pos[2]) * Eng->GetCurr(1,pos[0],pos[1],pos[2]));
energy+=fabs(Eng->GetVolt(0,pos[0],pos[1],pos[2]) * Eng->GetCurr(2,pos[0],pos[1],pos[2])); energy+=fabs(Eng->GetVolt(0,pos[0],pos[1],pos[2]) * Eng->GetCurr(2,pos[0],pos[1],pos[2]));
@ -274,15 +275,15 @@ void ProcessFields::WriteVTKCartesianGridHeader(ofstream &file, double const* co
file << "DATASET RECTILINEAR_GRID " << endl; file << "DATASET RECTILINEAR_GRID " << endl;
file << "DIMENSIONS " << numLines[0] << " " << numLines[1] << " " << numLines[2] << endl; file << "DIMENSIONS " << numLines[0] << " " << numLines[1] << " " << numLines[2] << endl;
file << "X_COORDINATES " << numLines[0] << " " << __VTK_DATA_TYPE__ << endl; file << "X_COORDINATES " << numLines[0] << " " << __VTK_DATA_TYPE__ << endl;
for (unsigned int i=0;i<numLines[0];++i) for (unsigned int i=0; i<numLines[0]; ++i)
file << setprecision(precision) << discLines[0][i] * discLines_scaling << " "; file << setprecision(precision) << discLines[0][i] * discLines_scaling << " ";
file << endl; file << endl;
file << "Y_COORDINATES " << numLines[1] << " " << __VTK_DATA_TYPE__ << endl; file << "Y_COORDINATES " << numLines[1] << " " << __VTK_DATA_TYPE__ << endl;
for (unsigned int i=0;i<numLines[1];++i) for (unsigned int i=0; i<numLines[1]; ++i)
file << setprecision(precision) << discLines[1][i] * discLines_scaling << " "; file << setprecision(precision) << discLines[1][i] * discLines_scaling << " ";
file << endl; file << endl;
file << "Z_COORDINATES " << numLines[2] << " " << __VTK_DATA_TYPE__ << endl; file << "Z_COORDINATES " << numLines[2] << " " << __VTK_DATA_TYPE__ << endl;
for (unsigned int i=0;i<numLines[2];++i) for (unsigned int i=0; i<numLines[2]; ++i)
file << setprecision(precision) << discLines[2][i] * discLines_scaling << " "; file << setprecision(precision) << discLines[2][i] * discLines_scaling << " ";
file << endl << endl; file << endl << endl;
file << "POINT_DATA " << numLines[0]*numLines[1]*numLines[2] << endl; file << "POINT_DATA " << numLines[0]*numLines[1]*numLines[2] << endl;
@ -299,9 +300,9 @@ void ProcessFields::WriteVTKCylindricalGridHeader(ofstream &file, double const*
file << "DATASET STRUCTURED_GRID " << endl; file << "DATASET STRUCTURED_GRID " << endl;
file << "DIMENSIONS " << numLines[0] << " " << numLines[1] << " " << numLines[2] << endl; file << "DIMENSIONS " << numLines[0] << " " << numLines[1] << " " << numLines[2] << endl;
file << "POINTS " << numLines[0]*numLines[1]*numLines[2] << " " << __VTK_DATA_TYPE__ << endl; file << "POINTS " << numLines[0]*numLines[1]*numLines[2] << " " << __VTK_DATA_TYPE__ << endl;
for (unsigned int k=0;k<numLines[2];++k) for (unsigned int k=0; k<numLines[2]; ++k)
for (unsigned int j=0;j<numLines[1];++j) for (unsigned int j=0; j<numLines[1]; ++j)
for (unsigned int i=0;i<numLines[0];++i) for (unsigned int i=0; i<numLines[0]; ++i)
{ {
file << setprecision(precision) << discLines[0][i] * cos(discLines[1][j]) * discLines_scaling << " " file << setprecision(precision) << discLines[0][i] * cos(discLines[1][j]) * discLines_scaling << " "
<< discLines[0][i] * sin(discLines[1][j]) * discLines_scaling << " " << discLines[0][i] * sin(discLines[1][j]) * discLines_scaling << " "
@ -321,13 +322,13 @@ void ProcessFields::WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT
meshT = CARTESIAN_MESH; //dump field components as they are... meshT = CARTESIAN_MESH; //dump field components as they are...
unsigned int pos[3]; unsigned int pos[3];
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
double cos_a = cos(discLines[1][pos[1]]); //needed only for CYLINDRICAL_MESH double cos_a = cos(discLines[1][pos[1]]); //needed only for CYLINDRICAL_MESH
double sin_a = sin(discLines[1][pos[1]]); //needed only for CYLINDRICAL_MESH double sin_a = sin(discLines[1][pos[1]]); //needed only for CYLINDRICAL_MESH
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
switch (meshT) switch (meshT)
{ {
@ -365,7 +366,7 @@ bool ProcessFields::DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT
bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT, double discLines_scaling) bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT, double discLines_scaling)
{ {
WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT, discLines_scaling); WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT, discLines_scaling);
for (unsigned int n=0;n<numFields;++n) for (unsigned int n=0; n<numFields; ++n)
{ {
WriteVTKVectorArray(file, names[n], array[n], discLines, numLines, precision, meshT); WriteVTKVectorArray(file, names[n], array[n], discLines, numLines, precision, meshT);
file << endl; file << endl;
@ -379,11 +380,11 @@ void ProcessFields::WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT
file << "LOOKUP_TABLE default" << endl; file << "LOOKUP_TABLE default" << endl;
unsigned int pos[3]; unsigned int pos[3];
int count=0; int count=0;
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
file << setprecision(precision) << array[pos[0]][pos[1]][pos[2]] << " "; file << setprecision(precision) << array[pos[0]][pos[1]][pos[2]] << " ";
++count; ++count;
@ -404,7 +405,7 @@ bool ProcessFields::DumpScalarArray2VTK(ofstream &file, string name, FDTD_FLOAT
bool ProcessFields::DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT, double discLines_scaling) bool ProcessFields::DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT const* const* const* const* array, unsigned int numFields, double const* const* discLines, unsigned int const* numLines, unsigned int precision, string header_info, MeshType meshT, double discLines_scaling)
{ {
WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT, discLines_scaling); WriteVTKHeader(file, discLines, numLines, precision, header_info, meshT, discLines_scaling);
for (unsigned int n=0;n<numFields;++n) for (unsigned int n=0; n<numFields; ++n)
{ {
WriteVTKScalarArray(file, names[n], array[n], numLines, precision); WriteVTKScalarArray(file, names[n], array[n], numLines, precision);
file << endl; file << endl;
@ -443,13 +444,13 @@ bool ProcessFields::DumpVectorArray2HDF5(string filename, string name, FDTD_FLOA
// at least it is save in case FDTD_FLOAT was defined as double... // at least it is save in case FDTD_FLOAT was defined as double...
// why does hdf5 write the dimensions backwards??? or matlab??? // why does hdf5 write the dimensions backwards??? or matlab???
float hdf5array[3][numLines[2]][numLines[1]][numLines[0]]; float hdf5array[3][numLines[2]][numLines[1]][numLines[0]];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
for (unsigned int i=0;i<numLines[0];++i) for (unsigned int i=0; i<numLines[0]; ++i)
{ {
for (unsigned int j=0;j<numLines[1];++j) for (unsigned int j=0; j<numLines[1]; ++j)
{ {
for (unsigned int k=0;k<numLines[2];++k) for (unsigned int k=0; k<numLines[2]; ++k)
{ {
hdf5array[n][k][j][i] = array[n][i][j][k]; hdf5array[n][k][j][i] = array[n][i][j][k];
} }

17
Common/processfields_td.cpp
View File

@ -60,13 +60,13 @@ int ProcessFieldsTD::Process()
FDTD_FLOAT**** field = Create_N_3DArray<FDTD_FLOAT>(numLines); FDTD_FLOAT**** field = Create_N_3DArray<FDTD_FLOAT>(numLines);
if (m_DumpType==E_FIELD_DUMP) if (m_DumpType==E_FIELD_DUMP)
{ {
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
OpPos[0]=start[0]+pos[0]*subSample[0]; OpPos[0]=start[0]+pos[0]*subSample[0];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
OpPos[1]=start[1]+pos[1]*subSample[1]; OpPos[1]=start[1]+pos[1]*subSample[1];
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
OpPos[2]=start[2]+pos[2]*subSample[2]; OpPos[2]=start[2]+pos[2]*subSample[2];
m_Eng_Interface->GetEField(OpPos,out); m_Eng_Interface->GetEField(OpPos,out);
@ -80,13 +80,13 @@ int ProcessFieldsTD::Process()
if (m_DumpType==H_FIELD_DUMP) if (m_DumpType==H_FIELD_DUMP)
{ {
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
OpPos[0]=start[0]+pos[0]*subSample[0]; OpPos[0]=start[0]+pos[0]*subSample[0];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
OpPos[1]=start[1]+pos[1]*subSample[1]; OpPos[1]=start[1]+pos[1]*subSample[1];
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
OpPos[2]=start[2]+pos[2]*subSample[2]; OpPos[2]=start[2]+pos[2]*subSample[2];
m_Eng_Interface->GetHField(OpPos,out); m_Eng_Interface->GetHField(OpPos,out);
@ -101,7 +101,10 @@ int ProcessFieldsTD::Process()
if (m_fileType==VTK_FILETYPE) if (m_fileType==VTK_FILETYPE)
{ {
ofstream file(filename.c_str()); ofstream file(filename.c_str());
if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;}; if (file.is_open()==false)
{
cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl;
};
DumpVectorArray2VTK(file,GetFieldNameByType(m_DumpType),field,discLines,numLines,m_precision,string("Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString(), m_Mesh_Type, discLines_scaling); DumpVectorArray2VTK(file,GetFieldNameByType(m_DumpType),field,discLines,numLines,m_precision,string("Interpolation: ")+m_Eng_Interface->GetInterpolationTypeString(), m_Mesh_Type, discLines_scaling);
file.close(); file.close();
} }

31
Common/processing.cpp
View File

@ -113,7 +113,7 @@ void Processing::AddStep(unsigned int step)
void Processing::AddSteps(vector<unsigned int> steps) void Processing::AddSteps(vector<unsigned int> steps)
{ {
for (size_t n=0;n<steps.size();++n) for (size_t n=0; n<steps.size(); ++n)
{ {
AddStep(steps.at(n)); AddStep(steps.at(n));
} }
@ -143,7 +143,7 @@ void Processing::AddFrequency(double freq)
void Processing::AddFrequency(vector<double> *freqs) void Processing::AddFrequency(vector<double> *freqs)
{ {
for (size_t n=0;n<freqs->size();++n) for (size_t n=0; n<freqs->size(); ++n)
{ {
AddFrequency(freqs->at(n)); AddFrequency(freqs->at(n));
} }
@ -156,7 +156,8 @@ void Processing::DefineStartStopCoord(double* dstart, double* dstop)
if (Op->SnapToMesh(dstop,stop,m_dualMesh,m_stop_inside)==false) if (Op->SnapToMesh(dstop,stop,m_dualMesh,m_stop_inside)==false)
cerr << "Processing::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl; cerr << "Processing::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl;
if (g_settings.showProbeDiscretization()) { if (g_settings.showProbeDiscretization())
{
cerr << m_Name << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], m_dualMesh ) << "," cerr << m_Name << ": snapped coords: (" << Op->GetDiscLine( 0, start[0], m_dualMesh ) << ","
<< Op->GetDiscLine( 1, start[1], m_dualMesh ) << "," << Op->GetDiscLine( 2, start[2], m_dualMesh ) << ") -> (" << Op->GetDiscLine( 1, start[1], m_dualMesh ) << "," << Op->GetDiscLine( 2, start[2], m_dualMesh ) << ") -> ("
<< Op->GetDiscLine( 0, stop[0], m_dualMesh ) << ","<< Op->GetDiscLine( 1, stop[1], m_dualMesh ) << "," << Op->GetDiscLine( 0, stop[0], m_dualMesh ) << ","<< Op->GetDiscLine( 1, stop[1], m_dualMesh ) << ","
@ -191,15 +192,18 @@ void Processing::DumpBox2File( string vtkfilenameprefix, bool dualMesh ) const
// normalize coordinates // normalize coordinates
double s1[3], s2[3]; double s1[3], s2[3];
for (int i=0; i<3; i++) { for (int i=0; i<3; i++)
{
s1[i] = min(Op->GetDiscLine(i,start[i],dualMesh),Op->GetDiscLine(i,stop[i],dualMesh)); s1[i] = min(Op->GetDiscLine(i,start[i],dualMesh),Op->GetDiscLine(i,stop[i],dualMesh));
s2[i] = max(Op->GetDiscLine(i,start[i],dualMesh),Op->GetDiscLine(i,stop[i],dualMesh)); s2[i] = max(Op->GetDiscLine(i,start[i],dualMesh),Op->GetDiscLine(i,stop[i],dualMesh));
} }
// fix degenerate box/plane -> line (paraview display problem) // fix degenerate box/plane -> line (paraview display problem)
if (((s1[0] == s2[0]) && (s1[1] == s2[1])) || ((s1[0] == s2[0]) && (s1[2] == s2[2])) || ((s1[2] == s2[2]) && (s1[1] == s2[1]))) { if (((s1[0] == s2[0]) && (s1[1] == s2[1])) || ((s1[0] == s2[0]) && (s1[2] == s2[2])) || ((s1[2] == s2[2]) && (s1[1] == s2[1])))
{
// line are not displayed correctly -> enlarge // line are not displayed correctly -> enlarge
for (int i=0; i<3; i++) { for (int i=0; i<3; i++)
{
double delta = min( Op->GetMeshDelta( i, start,dualMesh ), Op->GetMeshDelta( i, stop,dualMesh ) ) / Op->GetGridDelta() / 4.0; double delta = min( Op->GetMeshDelta( i, start,dualMesh ), Op->GetMeshDelta( i, stop,dualMesh ) ) / Op->GetGridDelta() / 4.0;
s1[i] -= delta; s1[i] -= delta;
s2[i] += delta; s2[i] += delta;
@ -209,7 +213,8 @@ void Processing::DumpBox2File( string vtkfilenameprefix, bool dualMesh ) const
// rescale coordinates // rescale coordinates
#ifndef OUTPUT_IN_DRAWINGUNITS #ifndef OUTPUT_IN_DRAWINGUNITS
double scaling = Op->GetGridDelta(); double scaling = Op->GetGridDelta();
for (int i=0; i<3; i++) { for (int i=0; i<3; i++)
{
s1[i] *= scaling; s1[i] *= scaling;
s2[i] *= scaling; s2[i] *= scaling;
} }
@ -258,7 +263,7 @@ void Processing::Dump_FD_Data(vector<double_complex> value, double factor, strin
time_t rawTime; time_t rawTime;
time(&rawTime); time(&rawTime);
file << "%dump by openEMS @" << ctime(&rawTime) << "%frequency\treal\timag\n"; file << "%dump by openEMS @" << ctime(&rawTime) << "%frequency\treal\timag\n";
for (size_t n=0;n<value.size();++n) for (size_t n=0; n<value.size(); ++n)
{ {
file << m_FD_Samples.at(n) << "\t" << 2.0 * std::real(value.at(n))*factor << "\t" << 2.0 * std::imag(value.at(n))*factor << "\n"; file << m_FD_Samples.at(n) << "\t" << 2.0 * std::real(value.at(n))*factor << "\t" << 2.0 * std::imag(value.at(n))*factor << "\n";
} }
@ -272,7 +277,7 @@ void ProcessingArray::AddProcessing(Processing* proc)
void ProcessingArray::FlushNext() void ProcessingArray::FlushNext()
{ {
for (size_t i=0;i<ProcessArray.size();++i) for (size_t i=0; i<ProcessArray.size(); ++i)
{ {
ProcessArray.at(i)->FlushNext(); ProcessArray.at(i)->FlushNext();
} }
@ -280,7 +285,7 @@ void ProcessingArray::FlushNext()
void ProcessingArray::Reset() void ProcessingArray::Reset()
{ {
for (size_t i=0;i<ProcessArray.size();++i) for (size_t i=0; i<ProcessArray.size(); ++i)
{ {
ProcessArray.at(i)->Reset(); ProcessArray.at(i)->Reset();
} }
@ -288,7 +293,7 @@ void ProcessingArray::Reset()
void ProcessingArray::DeleteAll() void ProcessingArray::DeleteAll()
{ {
for (size_t i=0;i<ProcessArray.size();++i) for (size_t i=0; i<ProcessArray.size(); ++i)
{ {
delete ProcessArray.at(i); delete ProcessArray.at(i);
} }
@ -299,7 +304,7 @@ int ProcessingArray::Process()
{ {
int nextProcess=maxInterval; int nextProcess=maxInterval;
//this could be done nicely in parallel?? //this could be done nicely in parallel??
for (size_t i=0;i<ProcessArray.size();++i) for (size_t i=0; i<ProcessArray.size(); ++i)
{ {
int step = ProcessArray.at(i)->Process(); int step = ProcessArray.at(i)->Process();
if ((step>0) && (step<nextProcess)) if ((step>0) && (step<nextProcess))
@ -310,6 +315,6 @@ int ProcessingArray::Process()
void ProcessingArray::DumpBoxes2File( string vtkfilenameprefix ) const void ProcessingArray::DumpBoxes2File( string vtkfilenameprefix ) const
{ {
for (size_t i=0;i<ProcessArray.size();++i) for (size_t i=0; i<ProcessArray.size(); ++i)
ProcessArray.at(i)->DumpBox2File( vtkfilenameprefix ); ProcessArray.at(i)->DumpBox2File( vtkfilenameprefix );
} }

6
Common/processintegral.cpp
View File

@ -37,7 +37,7 @@ void ProcessIntegral::InitProcess()
m_filename = m_Name; m_filename = m_Name;
OpenFile(m_filename); OpenFile(m_filename);
FD_Values.clear(); FD_Values.clear();
for (size_t n=0;n<m_FD_Samples.size();++n) for (size_t n=0; n<m_FD_Samples.size(); ++n)
FD_Values.push_back(0); FD_Values.push_back(0);
} }
@ -64,7 +64,7 @@ int ProcessIntegral::Process()
{ {
TD_Values.push_back(integral); TD_Values.push_back(integral);
file << setprecision(m_precision) << time; file << setprecision(m_precision) << time;
for (int n=0;n<NrInt;++n) for (int n=0; n<NrInt; ++n)
file << "\t" << m_Results[n] * m_weight; file << "\t" << m_Results[n] * m_weight;
file << endl; file << endl;
} }
@ -75,7 +75,7 @@ int ProcessIntegral::Process()
if (m_Eng_Interface->GetNumberOfTimesteps()%m_FD_Interval==0) if (m_Eng_Interface->GetNumberOfTimesteps()%m_FD_Interval==0)
{ {
double T = time; double T = time;
for (size_t n=0;n<m_FD_Samples.size();++n) for (size_t n=0; n<m_FD_Samples.size(); ++n)
{ {
FD_Values.at(n) += (double)integral * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T ); FD_Values.at(n) += (double)integral * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T );
} }

30
Common/processmodematch.cpp
View File

@ -22,7 +22,7 @@
ProcessModeMatch::ProcessModeMatch(Operator_Base* op) : ProcessIntegral(op) ProcessModeMatch::ProcessModeMatch(Operator_Base* op) : ProcessIntegral(op)
{ {
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
m_ModeParser[n] = new CSFunctionParser(); m_ModeParser[n] = new CSFunctionParser();
m_ModeDist[n] = NULL; m_ModeDist[n] = NULL;
@ -35,7 +35,7 @@ ProcessModeMatch::ProcessModeMatch(Operator_Base* op) : ProcessIntegral(op)
ProcessModeMatch::~ProcessModeMatch() ProcessModeMatch::~ProcessModeMatch()
{ {
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
delete m_ModeParser[n]; delete m_ModeParser[n];
m_ModeParser[n] = NULL; m_ModeParser[n] = NULL;
@ -61,7 +61,7 @@ void ProcessModeMatch::InitProcess()
m_Eng_Interface->SetInterpolationType(Engine_Interface_Base::NODE_INTERPOLATE); m_Eng_Interface->SetInterpolationType(Engine_Interface_Base::NODE_INTERPOLATE);
int Dump_Dim=0; int Dump_Dim=0;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if (start[n]>stop[n]) if (start[n]>stop[n])
{ {
@ -96,7 +96,7 @@ void ProcessModeMatch::InitProcess()
return; return;
} }
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
int ny = (m_ny+n+1)%3; int ny = (m_ny+n+1)%3;
int res = m_ModeParser[n]->Parse(m_ModeFunction[ny], "x,y,z,rho,a,r,t"); int res = m_ModeParser[n]->Parse(m_ModeFunction[ny], "x,y,z,rho,a,r,t");
@ -109,7 +109,7 @@ void ProcessModeMatch::InitProcess()
} }
} }
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
m_ModeDist[n] = Create2DArray<double>(m_numLines); m_ModeDist[n] = Create2DArray<double>(m_numLines);
} }
@ -122,11 +122,11 @@ void ProcessModeMatch::InitProcess()
discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],m_dualMesh); discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],m_dualMesh);
double norm = 0; double norm = 0;
double area = 0; double area = 0;
for (unsigned int posP = 0;posP<m_numLines[0];++posP) for (unsigned int posP = 0; posP<m_numLines[0]; ++posP)
{ {
pos[nP] = start[nP] + posP; pos[nP] = start[nP] + posP;
discLine[nP] = Op->GetDiscLine(nP,pos[nP],m_dualMesh); discLine[nP] = Op->GetDiscLine(nP,pos[nP],m_dualMesh);
for (unsigned int posPP = 0;posPP<m_numLines[1];++posPP) for (unsigned int posPP = 0; posPP<m_numLines[1]; ++posPP)
{ {
pos[nPP] = start[nPP] + posPP; pos[nPP] = start[nPP] + posPP;
discLine[nPP] = Op->GetDiscLine(nPP,pos[nPP],m_dualMesh); discLine[nPP] = Op->GetDiscLine(nPP,pos[nPP],m_dualMesh);
@ -149,7 +149,7 @@ void ProcessModeMatch::InitProcess()
var[6] = asin(1)-atan(var[2]/var[3]); //theta (t) var[6] = asin(1)-atan(var[2]/var[3]); //theta (t)
} }
area = Op->GetNodeArea(m_ny,pos,m_dualMesh); area = Op->GetNodeArea(m_ny,pos,m_dualMesh);
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
m_ModeDist[n][posP][posPP] = m_ModeParser[n]->Eval(var); //calc mode template m_ModeDist[n][posP][posPP] = m_ModeParser[n]->Eval(var); //calc mode template
if ((isnan(m_ModeDist[n][posP][posPP])) || (isinf(m_ModeDist[n][posP][posPP]))) if ((isnan(m_ModeDist[n][posP][posPP])) || (isinf(m_ModeDist[n][posP][posPP])))
@ -163,10 +163,10 @@ void ProcessModeMatch::InitProcess()
norm = sqrt(norm); norm = sqrt(norm);
// cerr << norm << endl; // cerr << norm << endl;
// normalize template function... // normalize template function...
for (unsigned int posP = 0;posP<m_numLines[0];++posP) for (unsigned int posP = 0; posP<m_numLines[0]; ++posP)
for (unsigned int posPP = 0;posPP<m_numLines[1];++posPP) for (unsigned int posPP = 0; posPP<m_numLines[1]; ++posPP)
{ {
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
m_ModeDist[n][posP][posPP] /= norm; m_ModeDist[n][posP][posPP] /= norm;
} }
@ -179,7 +179,7 @@ void ProcessModeMatch::InitProcess()
void ProcessModeMatch::Reset() void ProcessModeMatch::Reset()
{ {
ProcessIntegral::Reset(); ProcessIntegral::Reset();
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
Delete2DArray<double>(m_ModeDist[n],m_numLines); Delete2DArray<double>(m_ModeDist[n],m_numLines);
m_ModeDist[n] = NULL; m_ModeDist[n] = NULL;
@ -215,10 +215,10 @@ double* ProcessModeMatch::CalcMultipleIntegrals()
double out[3]={0,0,0}; double out[3]={0,0,0};
for (unsigned int posP = 0;posP<m_numLines[0];++posP) for (unsigned int posP = 0; posP<m_numLines[0]; ++posP)
{ {
pos[nP] = start[nP] + posP; pos[nP] = start[nP] + posP;
for (unsigned int posPP = 0;posPP<m_numLines[1];++posPP) for (unsigned int posPP = 0; posPP<m_numLines[1]; ++posPP)
{ {
pos[nPP] = start[nPP] + posPP; pos[nPP] = start[nPP] + posPP;
area = Op->GetNodeArea(m_ny,pos,m_dualMesh); area = Op->GetNodeArea(m_ny,pos,m_dualMesh);
@ -227,7 +227,7 @@ double* ProcessModeMatch::CalcMultipleIntegrals()
if (m_ModeFieldType==1) if (m_ModeFieldType==1)
m_Eng_Interface->GetHField(pos,out); m_Eng_Interface->GetHField(pos,out);
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
field = out[(m_ny+n+1)%3]; field = out[(m_ny+n+1)%3];
value += field * m_ModeDist[n][posP][posPP] * area; value += field * m_ModeDist[n][posP][posPP] * area;

36
FDTD/engine.cpp
View File

@ -35,7 +35,7 @@ Engine::Engine(const Operator* op)
m_type = BASIC; m_type = BASIC;
numTS = 0; numTS = 0;
Op = op; Op = op;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
numLines[n] = Op->GetOriginalNumLines(n); numLines[n] = Op->GetOriginalNumLines(n);
volt=NULL; volt=NULL;
curr=NULL; curr=NULL;
@ -61,7 +61,7 @@ void Engine::Init()
void Engine::InitExtensions() void Engine::InitExtensions()
{ {
for (size_t n=0;n<Op->GetNumberOfExtentions();++n) for (size_t n=0; n<Op->GetNumberOfExtentions(); ++n)
{ {
Operator_Extension* op_ext = Op->GetExtension(n); Operator_Extension* op_ext = Op->GetExtension(n);
Engine_Extension* eng_ext = op_ext->CreateEngineExtention(); Engine_Extension* eng_ext = op_ext->CreateEngineExtention();
@ -75,7 +75,7 @@ void Engine::InitExtensions()
void Engine::ClearExtensions() void Engine::ClearExtensions()
{ {
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
delete m_Eng_exts.at(n); delete m_Eng_exts.at(n);
m_Eng_exts.clear(); m_Eng_exts.clear();
} }
@ -111,13 +111,13 @@ void Engine::UpdateVoltages(unsigned int startX, unsigned int numX)
pos[0] = startX; pos[0] = startX;
//voltage updates //voltage updates
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
shift[0]=pos[0]; shift[0]=pos[0];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
shift[1]=pos[1]; shift[1]=pos[1];
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
shift[2]=pos[2]; shift[2]=pos[2];
//do the updates here //do the updates here
@ -144,7 +144,7 @@ void Engine::ApplyVoltageExcite()
unsigned int ny; unsigned int ny;
unsigned int pos[3]; unsigned int pos[3];
//soft voltage excitation here (E-field excite) //soft voltage excitation here (E-field excite)
for (unsigned int n=0;n<Op->Exc->Volt_Count;++n) for (unsigned int n=0; n<Op->Exc->Volt_Count; ++n)
{ {
exc_pos = (int)numTS - (int)Op->Exc->Volt_delay[n]; exc_pos = (int)numTS - (int)Op->Exc->Volt_delay[n];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length); exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
@ -165,11 +165,11 @@ void Engine::UpdateCurrents(unsigned int startX, unsigned int numX)
{ {
unsigned int pos[3]; unsigned int pos[3];
pos[0] = startX; pos[0] = startX;
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
for (pos[1]=0;pos[1]<numLines[1]-1;++pos[1]) for (pos[1]=0; pos[1]<numLines[1]-1; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2]-1;++pos[2]) for (pos[2]=0; pos[2]<numLines[2]-1; ++pos[2])
{ {
//do the updates here //do the updates here
//for x //for x
@ -195,7 +195,7 @@ void Engine::ApplyCurrentExcite()
unsigned int ny; unsigned int ny;
unsigned int pos[3]; unsigned int pos[3];
//soft current excitation here (H-field excite) //soft current excitation here (H-field excite)
for (unsigned int n=0;n<Op->Exc->Curr_Count;++n) for (unsigned int n=0; n<Op->Exc->Curr_Count; ++n)
{ {
exc_pos = (int)numTS - (int)Op->Exc->Curr_delay[n]; exc_pos = (int)numTS - (int)Op->Exc->Curr_delay[n];
exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length); exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length);
@ -215,7 +215,7 @@ void Engine::ApplyCurrentExcite()
void Engine::DoPreVoltageUpdates() void Engine::DoPreVoltageUpdates()
{ {
//execute extensions in reverse order -> highest priority gets access to the voltages last //execute extensions in reverse order -> highest priority gets access to the voltages last
for (int n=m_Eng_exts.size()-1;n>=0;--n) for (int n=m_Eng_exts.size()-1; n>=0; --n)
m_Eng_exts.at(n)->DoPreVoltageUpdates(); m_Eng_exts.at(n)->DoPreVoltageUpdates();
} }
@ -223,41 +223,41 @@ void Engine::DoPreVoltageUpdates()
void Engine::DoPostVoltageUpdates() void Engine::DoPostVoltageUpdates()
{ {
//execute extensions in normal order -> highest priority gets access to the voltages first //execute extensions in normal order -> highest priority gets access to the voltages first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
m_Eng_exts.at(n)->DoPostVoltageUpdates(); m_Eng_exts.at(n)->DoPostVoltageUpdates();
} }
void Engine::Apply2Voltages() void Engine::Apply2Voltages()
{ {
//execute extensions in normal order -> highest priority gets access to the voltages first //execute extensions in normal order -> highest priority gets access to the voltages first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
m_Eng_exts.at(n)->Apply2Voltages(); m_Eng_exts.at(n)->Apply2Voltages();
} }
void Engine::DoPreCurrentUpdates() void Engine::DoPreCurrentUpdates()
{ {
//execute extensions in reverse order -> highest priority gets access to the currents last //execute extensions in reverse order -> highest priority gets access to the currents last
for (int n=m_Eng_exts.size()-1;n>=0;--n) for (int n=m_Eng_exts.size()-1; n>=0; --n)
m_Eng_exts.at(n)->DoPreCurrentUpdates(); m_Eng_exts.at(n)->DoPreCurrentUpdates();
} }
void Engine::DoPostCurrentUpdates() void Engine::DoPostCurrentUpdates()
{ {
//execute extensions in normal order -> highest priority gets access to the currents first //execute extensions in normal order -> highest priority gets access to the currents first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
m_Eng_exts.at(n)->DoPostCurrentUpdates(); m_Eng_exts.at(n)->DoPostCurrentUpdates();
} }
void Engine::Apply2Current() void Engine::Apply2Current()
{ {
//execute extensions in normal order -> highest priority gets access to the currents first //execute extensions in normal order -> highest priority gets access to the currents first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
m_Eng_exts.at(n)->Apply2Current(); m_Eng_exts.at(n)->Apply2Current();
} }
bool Engine::IterateTS(unsigned int iterTS) bool Engine::IterateTS(unsigned int iterTS)
{ {
for (unsigned int iter=0;iter<iterTS;++iter) for (unsigned int iter=0; iter<iterTS; ++iter)
{ {
//voltage updates with extensions //voltage updates with extensions
DoPreVoltageUpdates(); DoPreVoltageUpdates();

5
FDTD/engine.h
View File

@ -21,8 +21,9 @@
#include <fstream> #include <fstream>
#include "operator.h" #include "operator.h"
namespace NS_Engine_Multithread { namespace NS_Engine_Multithread
class thread; // evil hack to access numTS from multithreading context {
class thread; // evil hack to access numTS from multithreading context
} }
class Engine_Extension; class Engine_Extension;

12
FDTD/engine_cylindermultigrid.cpp
View File

@ -47,7 +47,7 @@ Engine_CylinderMultiGrid::Engine_CylinderMultiGrid(const Operator_CylinderMultiG
m_InnerEng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync); m_InnerEng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync);
// if already has a base extension, switch places ... seems to be faster... // if already has a base extension, switch places ... seems to be faster...
for (size_t n=0;n<m_InnerEngine->m_Eng_exts.size();++n) for (size_t n=0; n<m_InnerEngine->m_Eng_exts.size(); ++n)
{ {
Engine_Ext_CylinderMultiGrid* eng_mg = dynamic_cast<Engine_Ext_CylinderMultiGrid*>(m_InnerEngine->m_Eng_exts.at(n)); Engine_Ext_CylinderMultiGrid* eng_mg = dynamic_cast<Engine_Ext_CylinderMultiGrid*>(m_InnerEngine->m_Eng_exts.at(n));
if (eng_mg) if (eng_mg)
@ -114,9 +114,9 @@ bool Engine_CylinderMultiGrid::IterateTS(unsigned int iterTS)
m_stopBarrier->wait(); //tell base and child to wait for another start event... m_stopBarrier->wait(); //tell base and child to wait for another start event...
//interpolate child data to base mesh... //interpolate child data to base mesh...
for (unsigned int n=0;n<Op_CMG->m_Split_Pos-1;++n) for (unsigned int n=0; n<Op_CMG->m_Split_Pos-1; ++n)
InterpolVoltChild2Base(n); InterpolVoltChild2Base(n);
for (unsigned int n=0;n<Op_CMG->m_Split_Pos-2;++n) for (unsigned int n=0; n<Op_CMG->m_Split_Pos-2; ++n)
InterpolCurrChild2Base(n); InterpolCurrChild2Base(n);
return true; return true;
@ -131,7 +131,7 @@ void Engine_CylinderMultiGrid::InterpolVoltChild2Base(unsigned int rzPlane)
pos[0] = rzPlane; pos[0] = rzPlane;
bool isOdd, isEven; bool isOdd, isEven;
f4vector half, one_eighth, three_eighth; f4vector half, one_eighth, three_eighth;
for (int n=0;n<4;++n) for (int n=0; n<4; ++n)
{ {
half.f[n]=0.5; half.f[n]=0.5;
one_eighth.f[n] = 1.0/8.0; one_eighth.f[n] = 1.0/8.0;
@ -195,7 +195,7 @@ void Engine_CylinderMultiGrid::InterpolCurrChild2Base(unsigned int rzPlane)
pos[0] = rzPlane; pos[0] = rzPlane;
bool isOdd, isEven; bool isOdd, isEven;
f4vector one_fourth, three_fourth; f4vector one_fourth, three_fourth;
for (int n=0;n<4;++n) for (int n=0; n<4; ++n)
{ {
one_fourth.f[n] = 1.0/4.0; one_fourth.f[n] = 1.0/4.0;
three_fourth.f[n] = 3.0/4.0; three_fourth.f[n] = 3.0/4.0;
@ -263,7 +263,7 @@ void Engine_CylinderMultiGrid_Thread::operator()()
{ {
m_startBarrier->wait(); //wait for Base engine to start the iterations... m_startBarrier->wait(); //wait for Base engine to start the iterations...
while(*m_numTS>0) //m_numTS==0 request to terminate this thread... while (*m_numTS>0) //m_numTS==0 request to terminate this thread...
{ {
if (m_isBase) if (m_isBase)
m_Eng->Engine_Multithread::IterateTS(*m_numTS); m_Eng->Engine_Multithread::IterateTS(*m_numTS);

16
FDTD/engine_interface_fdtd.cpp
View File

@ -42,7 +42,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
{ {
default: default:
case NO_INTERPOLATION: case NO_INTERPOLATION:
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delta = m_Op->GetEdgeLength(n,pos,false); delta = m_Op->GetEdgeLength(n,pos,false);
if (delta) if (delta)
@ -52,7 +52,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
} }
break; break;
case NODE_INTERPOLATE: case NODE_INTERPOLATE:
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delta = m_Op->GetMeshDelta(n,iPos); delta = m_Op->GetMeshDelta(n,iPos);
out[n] = m_Eng->GetVolt(n,iPos); out[n] = m_Eng->GetVolt(n,iPos);
@ -74,7 +74,7 @@ double* Engine_Interface_FDTD::GetEField(const unsigned int* pos, double* out) c
} }
break; break;
case CELL_INTERPOLATE: case CELL_INTERPOLATE:
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
nP = (n+1)%3; nP = (n+1)%3;
nPP = (n+2)%3; nPP = (n+2)%3;
@ -120,7 +120,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
out[2] = m_Eng->GetCurr(2,pos) / m_Op->GetEdgeLength(2,pos,true); out[2] = m_Eng->GetCurr(2,pos) / m_Op->GetEdgeLength(2,pos,true);
break; break;
case NODE_INTERPOLATE: case NODE_INTERPOLATE:
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
nP = (n+1)%3; nP = (n+1)%3;
nPP = (n+2)%3; nPP = (n+2)%3;
@ -141,7 +141,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
} }
break; break;
case CELL_INTERPOLATE: case CELL_INTERPOLATE:
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delta = m_Op->GetMeshDelta(n,iPos,true); delta = m_Op->GetMeshDelta(n,iPos,true);
out[n] = m_Eng->GetCurr(n,iPos); out[n] = m_Eng->GetCurr(n,iPos);
@ -165,18 +165,18 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const
{ {
double result=0; double result=0;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if (start[n]<stop[n]) if (start[n]<stop[n])
{ {
unsigned int pos[3]={start[0],start[1],start[2]}; unsigned int pos[3]={start[0],start[1],start[2]};
for (;pos[n]<stop[n];++pos[n]) for (; pos[n]<stop[n]; ++pos[n])
result += m_Eng->GetVolt(n,pos[0],pos[1],pos[2]); result += m_Eng->GetVolt(n,pos[0],pos[1],pos[2]);
} }
else else
{ {
unsigned int pos[3]={stop[0],stop[1],stop[2]}; unsigned int pos[3]={stop[0],stop[1],stop[2]};
for (;pos[n]<start[n];++pos[n]) for (; pos[n]<start[n]; ++pos[n])
result -= m_Eng->GetVolt(n,pos[0],pos[1],pos[2]); result -= m_Eng->GetVolt(n,pos[0],pos[1],pos[2]);
} }
} }

41
FDTD/engine_multithread.cpp
View File

@ -18,9 +18,9 @@
//#define ENABLE_DEBUG_TIME //#define ENABLE_DEBUG_TIME
#ifdef ENABLE_DEBUG_TIME #ifdef ENABLE_DEBUG_TIME
#define DEBUG_TIME(x) x; #define DEBUG_TIME(x) x;
#else #else
#define DEBUG_TIME(x) ; #define DEBUG_TIME(x) ;
#endif #endif
@ -58,10 +58,12 @@ Engine_Multithread::~Engine_Multithread()
#ifdef ENABLE_DEBUG_TIME #ifdef ENABLE_DEBUG_TIME
NS_Engine_Multithread::DBG().cout() << "Engine_Multithread::~Engine_Multithread()" << endl; NS_Engine_Multithread::DBG().cout() << "Engine_Multithread::~Engine_Multithread()" << endl;
std::map<boost::thread::id, std::vector<double> >::iterator it; std::map<boost::thread::id, std::vector<double> >::iterator it;
for (it=m_timer_list.begin(); it!=m_timer_list.end(); it++) { for (it=m_timer_list.begin(); it!=m_timer_list.end(); it++)
{
NS_Engine_Multithread::DBG().cout() << "*** DEBUG Thread: " << it->first << std::endl; NS_Engine_Multithread::DBG().cout() << "*** DEBUG Thread: " << it->first << std::endl;
std::vector<double>::iterator it2; std::vector<double>::iterator it2;
for (it2=it->second.begin(); it2<it->second.end();) { for (it2=it->second.begin(); it2<it->second.end();)
{
NS_Engine_Multithread::DBG().cout() << "after voltage update, before barrier1: " << fixed << setprecision(6) << *(it2++) << std::endl; NS_Engine_Multithread::DBG().cout() << "after voltage update, before barrier1: " << fixed << setprecision(6) << *(it2++) << std::endl;
NS_Engine_Multithread::DBG().cout() << "after barrier1, before barrier2: " << fixed << setprecision(6) << *(it2++) << std::endl; NS_Engine_Multithread::DBG().cout() << "after barrier1, before barrier2: " << fixed << setprecision(6) << *(it2++) << std::endl;
NS_Engine_Multithread::DBG().cout() << "after barrier2, before current update: " << fixed << setprecision(6) << *(it2++) << std::endl; NS_Engine_Multithread::DBG().cout() << "after barrier2, before current update: " << fixed << setprecision(6) << *(it2++) << std::endl;
@ -117,7 +119,7 @@ void Engine_Multithread::Init()
m_thread_group.add_thread( t ); m_thread_group.add_thread( t );
} }
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
m_Eng_exts.at(n)->SetNumberOfThreads(m_numThreads); m_Eng_exts.at(n)->SetNumberOfThreads(m_numThreads);
} }
@ -134,9 +136,12 @@ void Engine_Multithread::Reset()
m_stopThreads = true; m_stopThreads = true;
m_stopBarrier->wait(); // wait for the threads to finish m_stopBarrier->wait(); // wait for the threads to finish
m_thread_group.join_all(); // wait for termination m_thread_group.join_all(); // wait for termination
delete m_IterateBarrier; m_IterateBarrier = 0; delete m_IterateBarrier;
delete m_startBarrier; m_startBarrier = 0; m_IterateBarrier = 0;
delete m_stopBarrier; m_stopBarrier = 0; delete m_startBarrier;
m_startBarrier = 0;
delete m_stopBarrier;
m_stopBarrier = 0;
} }
Engine_SSE_Compressed::Reset(); Engine_SSE_Compressed::Reset();
@ -158,7 +163,7 @@ bool Engine_Multithread::IterateTS(unsigned int iterTS)
void Engine_Multithread::DoPreVoltageUpdates(int threadID) void Engine_Multithread::DoPreVoltageUpdates(int threadID)
{ {
//execute extensions in reverse order -> highest priority gets access to the voltages last //execute extensions in reverse order -> highest priority gets access to the voltages last
for (int n=m_Eng_exts.size()-1;n>=0;--n) for (int n=m_Eng_exts.size()-1; n>=0; --n)
{ {
m_Eng_exts.at(n)->DoPreVoltageUpdates(threadID); m_Eng_exts.at(n)->DoPreVoltageUpdates(threadID);
m_IterateBarrier->wait(); m_IterateBarrier->wait();
@ -169,7 +174,7 @@ void Engine_Multithread::DoPreVoltageUpdates(int threadID)
void Engine_Multithread::DoPostVoltageUpdates(int threadID) void Engine_Multithread::DoPostVoltageUpdates(int threadID)
{ {
//execute extensions in normal order -> highest priority gets access to the voltages first //execute extensions in normal order -> highest priority gets access to the voltages first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
{ {
m_Eng_exts.at(n)->DoPostVoltageUpdates(threadID); m_Eng_exts.at(n)->DoPostVoltageUpdates(threadID);
m_IterateBarrier->wait(); m_IterateBarrier->wait();
@ -179,7 +184,7 @@ void Engine_Multithread::DoPostVoltageUpdates(int threadID)
void Engine_Multithread::Apply2Voltages(int threadID) void Engine_Multithread::Apply2Voltages(int threadID)
{ {
//execute extensions in normal order -> highest priority gets access to the voltages first //execute extensions in normal order -> highest priority gets access to the voltages first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
{ {
m_Eng_exts.at(n)->Apply2Voltages(threadID); m_Eng_exts.at(n)->Apply2Voltages(threadID);
m_IterateBarrier->wait(); m_IterateBarrier->wait();
@ -189,7 +194,7 @@ void Engine_Multithread::Apply2Voltages(int threadID)
void Engine_Multithread::DoPreCurrentUpdates(int threadID) void Engine_Multithread::DoPreCurrentUpdates(int threadID)
{ {
//execute extensions in reverse order -> highest priority gets access to the currents last //execute extensions in reverse order -> highest priority gets access to the currents last
for (int n=m_Eng_exts.size()-1;n>=0;--n) for (int n=m_Eng_exts.size()-1; n>=0; --n)
{ {
m_Eng_exts.at(n)->DoPreCurrentUpdates(threadID); m_Eng_exts.at(n)->DoPreCurrentUpdates(threadID);
m_IterateBarrier->wait(); m_IterateBarrier->wait();
@ -199,7 +204,7 @@ void Engine_Multithread::DoPreCurrentUpdates(int threadID)
void Engine_Multithread::DoPostCurrentUpdates(int threadID) void Engine_Multithread::DoPostCurrentUpdates(int threadID)
{ {
//execute extensions in normal order -> highest priority gets access to the currents first //execute extensions in normal order -> highest priority gets access to the currents first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
{ {
m_Eng_exts.at(n)->DoPostCurrentUpdates(threadID); m_Eng_exts.at(n)->DoPostCurrentUpdates(threadID);
m_IterateBarrier->wait(); m_IterateBarrier->wait();
@ -209,7 +214,7 @@ void Engine_Multithread::DoPostCurrentUpdates(int threadID)
void Engine_Multithread::Apply2Current(int threadID) void Engine_Multithread::Apply2Current(int threadID)
{ {
//execute extensions in normal order -> highest priority gets access to the currents first //execute extensions in normal order -> highest priority gets access to the currents first
for (size_t n=0;n<m_Eng_exts.size();++n) for (size_t n=0; n<m_Eng_exts.size(); ++n)
{ {
m_Eng_exts.at(n)->Apply2Current(threadID); m_Eng_exts.at(n)->Apply2Current(threadID);
m_IterateBarrier->wait(); m_IterateBarrier->wait();
@ -219,7 +224,8 @@ void Engine_Multithread::Apply2Current(int threadID)
// //
// ************************************************************************************************************************* // *************************************************************************************************************************
// //
namespace NS_Engine_Multithread { namespace NS_Engine_Multithread
{
thread::thread( Engine_Multithread* ptr, unsigned int start, unsigned int stop, unsigned int stop_h, unsigned int threadID ) thread::thread( Engine_Multithread* ptr, unsigned int start, unsigned int stop, unsigned int stop_h, unsigned int threadID )
{ {
@ -236,7 +242,8 @@ void thread::operator()()
//DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") started." << endl; //DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") started." << endl;
while (!m_enginePtr->m_stopThreads) { while (!m_enginePtr->m_stopThreads)
{
// wait for start // wait for start
//DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") waiting..." << endl; //DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") waiting..." << endl;
m_enginePtr->m_startBarrier->wait(); m_enginePtr->m_startBarrier->wait();
@ -244,7 +251,7 @@ void thread::operator()()
DEBUG_TIME( Timer timer1 ); DEBUG_TIME( Timer timer1 );
for (unsigned int iter=0;iter<m_enginePtr->m_iterTS;++iter) for (unsigned int iter=0; iter<m_enginePtr->m_iterTS; ++iter)
{ {
// pre voltage stuff... // pre voltage stuff...
m_enginePtr->DoPreVoltageUpdates(m_threadID); m_enginePtr->DoPreVoltageUpdates(m_threadID);

30
FDTD/engine_multithread.h
View File

@ -34,34 +34,38 @@
class Engine_Multithread; class Engine_Multithread;
namespace NS_Engine_Multithread { namespace NS_Engine_Multithread
{
class DBG { // debug class DBG // debug
public: {
public:
DBG() {} DBG() {}
~DBG() { std::cout << os.str();} ~DBG() { std::cout << os.str();}
std::ostringstream& cout() {return os;} std::ostringstream& cout() {return os;}
protected: protected:
std::ostringstream os; std::ostringstream os;
}; };
class Timer { //debug class Timer //debug
public: {
public:
Timer() {gettimeofday(&t1,NULL);} Timer() {gettimeofday(&t1,NULL);}
double elapsed() {gettimeofday(&t2,NULL); return (t2.tv_sec-t1.tv_sec) + (t2.tv_usec-t1.tv_usec)*1e-6;} double elapsed() {gettimeofday(&t2,NULL); return (t2.tv_sec-t1.tv_sec) + (t2.tv_usec-t1.tv_usec)*1e-6;}
protected: protected:
timeval t1,t2; timeval t1,t2;
}; };
class thread { class thread
public: {
public:
thread( Engine_Multithread* ptr, unsigned int start, unsigned int stop, unsigned int stop_h, unsigned int threadID ); thread( Engine_Multithread* ptr, unsigned int start, unsigned int stop, unsigned int stop_h, unsigned int threadID );
void operator()(); void operator()();
protected: protected:
unsigned int m_start, m_stop, m_stop_h, m_threadID; unsigned int m_start, m_stop, m_stop_h, m_threadID;
Engine_Multithread *m_enginePtr; Engine_Multithread *m_enginePtr;
}; };
} // namespace } // namespace

12
FDTD/engine_sse.cpp
View File

@ -79,13 +79,13 @@ void Engine_sse::UpdateVoltages(unsigned int startX, unsigned int numX)
f4vector temp; f4vector temp;
pos[0] = startX; pos[0] = startX;
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
shift[0]=pos[0]; shift[0]=pos[0];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
shift[1]=pos[1]; shift[1]=pos[1];
for (pos[2]=1;pos[2]<numVectors;++pos[2]) for (pos[2]=1; pos[2]<numVectors; ++pos[2])
{ {
// x-polarization // x-polarization
f4_volt[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[0][pos[0]][pos[1]][pos[2]].v; f4_volt[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[0][pos[0]][pos[1]][pos[2]].v;
@ -131,11 +131,11 @@ void Engine_sse::UpdateCurrents(unsigned int startX, unsigned int numX)
f4vector temp; f4vector temp;
pos[0] = startX; pos[0] = startX;
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
for (pos[1]=0;pos[1]<numLines[1]-1;++pos[1]) for (pos[1]=0; pos[1]<numLines[1]-1; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numVectors-1;++pos[2]) for (pos[2]=0; pos[2]<numVectors-1; ++pos[2])
{ {
// x-pol // x-pol
f4_curr[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[0][pos[0]][pos[1]][pos[2]].v; f4_curr[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[0][pos[0]][pos[1]][pos[2]].v;

14
FDTD/engine_sse_compressed.cpp
View File

@ -17,7 +17,7 @@
#include "engine_sse_compressed.h" #include "engine_sse_compressed.h"
#ifdef __SSE2__ #ifdef __SSE2__
#include <emmintrin.h> #include <emmintrin.h>
#endif #endif
Engine_SSE_Compressed* Engine_SSE_Compressed::New(const Operator_SSE_Compressed* op) Engine_SSE_Compressed* Engine_SSE_Compressed::New(const Operator_SSE_Compressed* op)
@ -45,13 +45,13 @@ void Engine_SSE_Compressed::UpdateVoltages(unsigned int startX, unsigned int num
pos[0] = startX; pos[0] = startX;
unsigned int index=0; unsigned int index=0;
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
shift[0]=pos[0]; shift[0]=pos[0];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
shift[1]=pos[1]; shift[1]=pos[1];
for (pos[2]=1;pos[2]<numVectors;++pos[2]) for (pos[2]=1; pos[2]<numVectors; ++pos[2])
{ {
index = Op->m_Op_index[pos[0]][pos[1]][pos[2]]; index = Op->m_Op_index[pos[0]][pos[1]][pos[2]];
// x-polarization // x-polarization
@ -108,11 +108,11 @@ void Engine_SSE_Compressed::UpdateCurrents(unsigned int startX, unsigned int num
pos[0] = startX; pos[0] = startX;
unsigned int index; unsigned int index;
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
for (pos[1]=0;pos[1]<numLines[1]-1;++pos[1]) for (pos[1]=0; pos[1]<numLines[1]-1; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numVectors-1;++pos[2]) for (pos[2]=0; pos[2]<numVectors-1; ++pos[2])
{ {
index = Op->m_Op_index[pos[0]][pos[1]][pos[2]]; index = Op->m_Op_index[pos[0]][pos[1]][pos[2]];
// x-pol // x-pol

19
FDTD/excitation.cpp
View File

@ -34,7 +34,8 @@ Excitation::Excitation( double timestep )
Curr_dir = 0; Curr_dir = 0;
Curr_Count = 0; Curr_Count = 0;
for (int n=0;n<3;++n) { for (int n=0; n<3; ++n)
{
Volt_index[n] = 0; Volt_index[n] = 0;
Curr_index[n] = 0; Curr_index[n] = 0;
Volt_Count_Dir[n] = 0; Volt_Count_Dir[n] = 0;
@ -55,7 +56,8 @@ Excitation::~Excitation()
delete[] Curr_delay; delete[] Curr_delay;
delete[] Curr_dir; delete[] Curr_dir;
delete[] Curr_amp; delete[] Curr_amp;
for (int n=0;n<3;++n) { for (int n=0; n<3; ++n)
{
delete[] Volt_index[n]; delete[] Volt_index[n];
delete[] Curr_index[n]; delete[] Curr_index[n];
} }
@ -64,7 +66,8 @@ Excitation::~Excitation()
bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS ) bool Excitation::setupExcitation( TiXmlElement* Excite, unsigned int maxTS )
{ {
if (!Excite) { if (!Excite)
{
cerr << "Can't read openEMS excitation settings... " << endl; cerr << "Can't read openEMS excitation settings... " << endl;
return false; return false;
} }
@ -113,7 +116,7 @@ unsigned int Excitation::GetMaxExcitationTimestep() const
{ {
FDTD_FLOAT maxAmp=0; FDTD_FLOAT maxAmp=0;
unsigned int maxStep=0; unsigned int maxStep=0;
for (unsigned int n=1;n<Length+1;++n) for (unsigned int n=1; n<Length+1; ++n)
{ {
if (fabs(Signal_volt[n])>maxAmp) if (fabs(Signal_volt[n])>maxAmp)
{ {
@ -141,7 +144,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
Signal_curr = new FDTD_FLOAT[Length+1]; Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0; Signal_volt[0]=0.0;
Signal_curr[0]=0.0; Signal_curr[0]=0.0;
for (unsigned int n=1;n<Length+1;++n) for (unsigned int n=1; n<Length+1; ++n)
{ {
double t = (n-1)*dT; 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)); 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));
@ -210,7 +213,7 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
exit(1); exit(1);
} }
double vars[1]; double vars[1];
for (unsigned int n=1;n<Length+1;++n) for (unsigned int n=1; n<Length+1; ++n)
{ {
vars[0] = (n-1)*dT; vars[0] = (n-1)*dT;
Signal_volt[n] = fParse.Eval(vars); Signal_volt[n] = fParse.Eval(vars);
@ -234,7 +237,7 @@ void Excitation::CalcSinusExcitation(double f0, int nTS)
Signal_curr = new FDTD_FLOAT[Length+1]; Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0; Signal_volt[0]=0.0;
Signal_curr[0]=0.0; Signal_curr[0]=0.0;
for (unsigned int n=1;n<Length+1;++n) for (unsigned int n=1; n<Length+1; ++n)
{ {
double t = (n-1)*dT; double t = (n-1)*dT;
Signal_volt[n] = sin(2.0*PI*f0*t); Signal_volt[n] = sin(2.0*PI*f0*t);
@ -297,7 +300,7 @@ void Excitation::setupCurrentExcitation( vector<unsigned int> const curr_vIndex[
// cerr << "Excitation::setupCurrentExcitation(): Number of current excitation points: " << Curr_Count << endl; // cerr << "Excitation::setupCurrentExcitation(): Number of current excitation points: " << Curr_Count << endl;
// if (Curr_Count==0) // if (Curr_Count==0)
// cerr << "No H-Field/current excitation found!" << endl; // cerr << "No H-Field/current excitation found!" << endl;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
for (unsigned int i=0; i<Curr_Count; i++) for (unsigned int i=0; i<Curr_Count; i++)
Curr_index[n][i] = curr_vIndex[n].at(i); Curr_index[n][i] = curr_vIndex[n].at(i);
for (unsigned int i=0; i<Curr_Count; i++) for (unsigned int i=0; i<Curr_Count; i++)

18
FDTD/extensions/engine_ext_cylinder.cpp
View File

@ -26,7 +26,7 @@ Engine_Ext_Cylinder::Engine_Ext_Cylinder(Operator_Ext_Cylinder* op_ext) : Engine
CC_closedAlpha = op_ext->CC_closedAlpha; CC_closedAlpha = op_ext->CC_closedAlpha;
CC_R0_included = op_ext->CC_R0_included; CC_R0_included = op_ext->CC_R0_included;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
numLines[n] = op_ext->m_Op->GetOriginalNumLines(n); numLines[n] = op_ext->m_Op->GetOriginalNumLines(n);
//this cylindrical extension should be executed first? //this cylindrical extension should be executed first?
@ -41,17 +41,17 @@ void Engine_Ext_Cylinder::Apply2Voltages()
{ {
unsigned int pos[3]; unsigned int pos[3];
pos[0] = 0; pos[0] = 0;
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2])*cyl_Op->vv_R0[pos[2]]); m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2])*cyl_Op->vv_R0[pos[2]]);
for (pos[1]=0;pos[1]<numLines[1]-1;++pos[1]) for (pos[1]=0; pos[1]<numLines[1]-1; ++pos[1])
{ {
m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2]) + cyl_Op->vi_R0[pos[2]] * m_Eng->GetCurr(1,0,pos[1],pos[2]) ); m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2]) + cyl_Op->vi_R0[pos[2]] * m_Eng->GetCurr(1,0,pos[1],pos[2]) );
} }
} }
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
m_Eng->SetVolt(1,0,pos[1],pos[2], 0); //no voltage in alpha-direction at r=0 m_Eng->SetVolt(1,0,pos[1],pos[2], 0); //no voltage in alpha-direction at r=0
m_Eng->SetVolt(2,0,pos[1],pos[2], m_Eng->GetVolt(2,0,0,pos[2]) ); m_Eng->SetVolt(2,0,pos[1],pos[2], m_Eng->GetVolt(2,0,0,pos[2]) );
@ -63,9 +63,9 @@ void Engine_Ext_Cylinder::Apply2Voltages()
unsigned int pos[3]; unsigned int pos[3];
// copy voltages from last alpha-plane to first // copy voltages from last alpha-plane to first
unsigned int last_A_Line = numLines[1]-1; unsigned int last_A_Line = numLines[1]-1;
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
m_Eng->SetVolt(0,pos[0],0,pos[2], m_Eng->GetVolt(0,pos[0],last_A_Line,pos[2]) ); m_Eng->SetVolt(0,pos[0],0,pos[2], m_Eng->GetVolt(0,pos[0],last_A_Line,pos[2]) );
m_Eng->SetVolt(1,pos[0],0,pos[2], m_Eng->GetVolt(1,pos[0],last_A_Line,pos[2]) ); m_Eng->SetVolt(1,pos[0],0,pos[2], m_Eng->GetVolt(1,pos[0],last_A_Line,pos[2]) );
@ -81,10 +81,10 @@ void Engine_Ext_Cylinder::Apply2Current()
//close alpha //close alpha
unsigned int pos[3]; unsigned int pos[3];
// copy currents from first alpha-plane to last // copy currents from first alpha-plane to last
for (pos[0]=0;pos[0]<numLines[0]-1;++pos[0]) for (pos[0]=0; pos[0]<numLines[0]-1; ++pos[0])
{ {
unsigned int last_A_Line = numLines[1]-1; unsigned int last_A_Line = numLines[1]-1;
for (pos[2]=0;pos[2]<numLines[2]-1;++pos[2]) for (pos[2]=0; pos[2]<numLines[2]-1; ++pos[2])
{ {
m_Eng->SetCurr(0,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(0,pos[0],0,pos[2]) ); m_Eng->SetCurr(0,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(0,pos[0],0,pos[2]) );
m_Eng->SetCurr(1,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(1,pos[0],0,pos[2]) ); m_Eng->SetCurr(1,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(1,pos[0],0,pos[2]) );

6
FDTD/extensions/engine_ext_cylindermultigrid.cpp
View File

@ -96,7 +96,11 @@ void Engine_Ext_CylinderMultiGrid::SyncVoltages()
unsigned int pos[3]; unsigned int pos[3];
pos[0] = m_Eng_MG->Op_CMG->GetSplitPos()-1; pos[0] = m_Eng_MG->Op_CMG->GetSplitPos()-1;
unsigned int pos1_half = 0; unsigned int pos1_half = 0;
f4vector v_null; v_null.f[0] = 0;v_null.f[1] = 0;v_null.f[2] = 0;v_null.f[3] = 0; f4vector v_null;
v_null.f[0] = 0;
v_null.f[1] = 0;
v_null.f[2] = 0;
v_null.f[3] = 0;
for (pos[1]=0; pos[1]<numLines[1]-1; pos[1]+=2) for (pos[1]=0; pos[1]<numLines[1]-1; pos[1]+=2)
{ {
pos1_half = pos[1]/2; pos1_half = pos[1]/2;

20
FDTD/extensions/engine_ext_dispersive.cpp
View File

@ -22,12 +22,12 @@
Engine_Ext_Dispersive::Engine_Ext_Dispersive(Operator_Ext_Dispersive* op_ext_disp) : Engine_Extension(op_ext_disp) Engine_Ext_Dispersive::Engine_Ext_Dispersive(Operator_Ext_Dispersive* op_ext_disp) : Engine_Extension(op_ext_disp)
{ {
m_Op_Ext_Disp = op_ext_disp; m_Op_Ext_Disp = op_ext_disp;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if (m_Op_Ext_Disp->m_curr_ADE_On==true) if (m_Op_Ext_Disp->m_curr_ADE_On==true)
{ {
curr_ADE[n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count]; curr_ADE[n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count];
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
curr_ADE[n][i]=0.0; curr_ADE[n][i]=0.0;
} }
else else
@ -35,7 +35,7 @@ Engine_Ext_Dispersive::Engine_Ext_Dispersive(Operator_Ext_Dispersive* op_ext_dis
if (m_Op_Ext_Disp->m_volt_ADE_On==true) if (m_Op_Ext_Disp->m_volt_ADE_On==true)
{ {
volt_ADE[n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count]; volt_ADE[n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count];
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
volt_ADE[n][i]=0.0; volt_ADE[n][i]=0.0;
} }
else else
@ -45,7 +45,7 @@ Engine_Ext_Dispersive::Engine_Ext_Dispersive(Operator_Ext_Dispersive* op_ext_dis
Engine_Ext_Dispersive::~Engine_Ext_Dispersive() Engine_Ext_Dispersive::~Engine_Ext_Dispersive()
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delete[] curr_ADE[n]; delete[] curr_ADE[n];
curr_ADE[n] = NULL; curr_ADE[n] = NULL;
@ -65,7 +65,7 @@ void Engine_Ext_Dispersive::Apply2Voltages()
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
{ {
m_Eng->Engine::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]); m_Eng->Engine::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
m_Eng->Engine::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]); m_Eng->Engine::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
@ -76,7 +76,7 @@ void Engine_Ext_Dispersive::Apply2Voltages()
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*)m_Eng; Engine_sse* eng_sse = (Engine_sse*)m_Eng;
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
{ {
eng_sse->Engine_sse::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]); eng_sse->Engine_sse::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
eng_sse->Engine_sse::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]); eng_sse->Engine_sse::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
@ -85,7 +85,7 @@ void Engine_Ext_Dispersive::Apply2Voltages()
break; break;
} }
default: default:
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
{ {
m_Eng->SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]); m_Eng->SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
m_Eng->SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]); m_Eng->SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
@ -106,7 +106,7 @@ void Engine_Ext_Dispersive::Apply2Current()
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
{ {
m_Eng->Engine::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]); m_Eng->Engine::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
m_Eng->Engine::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]); m_Eng->Engine::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
@ -117,7 +117,7 @@ void Engine_Ext_Dispersive::Apply2Current()
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*)m_Eng; Engine_sse* eng_sse = (Engine_sse*)m_Eng;
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
{ {
eng_sse->Engine_sse::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]); eng_sse->Engine_sse::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
eng_sse->Engine_sse::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]); eng_sse->Engine_sse::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
@ -126,7 +126,7 @@ void Engine_Ext_Dispersive::Apply2Current()
break; break;
} }
default: default:
for (unsigned int i=0;i<m_Op_Ext_Disp->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
{ {
m_Eng->SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]); m_Eng->SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
m_Eng->SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]); m_Eng->SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);

12
FDTD/extensions/engine_ext_lorentzmaterial.cpp
View File

@ -40,7 +40,7 @@ void Engine_Ext_LorentzMaterial::DoPreVoltageUpdates()
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int i=0;i<m_Op_Ext_Lor->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
{ {
volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i]; volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i];
volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
@ -56,7 +56,7 @@ void Engine_Ext_LorentzMaterial::DoPreVoltageUpdates()
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*)m_Eng; Engine_sse* eng_sse = (Engine_sse*)m_Eng;
for (unsigned int i=0;i<m_Op_Ext_Lor->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
{ {
volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i]; volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i];
volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
@ -70,7 +70,7 @@ void Engine_Ext_LorentzMaterial::DoPreVoltageUpdates()
break; break;
} }
default: default:
for (unsigned int i=0;i<m_Op_Ext_Lor->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
{ {
volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i]; volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i];
volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]); volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
@ -96,7 +96,7 @@ void Engine_Ext_LorentzMaterial::DoPreCurrentUpdates()
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int i=0;i<m_Op_Ext_Lor->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
{ {
curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i]; curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i];
curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
@ -112,7 +112,7 @@ void Engine_Ext_LorentzMaterial::DoPreCurrentUpdates()
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*)m_Eng; Engine_sse* eng_sse = (Engine_sse*)m_Eng;
for (unsigned int i=0;i<m_Op_Ext_Lor->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
{ {
curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i]; curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i];
curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
@ -126,7 +126,7 @@ void Engine_Ext_LorentzMaterial::DoPreCurrentUpdates()
break; break;
} }
default: default:
for (unsigned int i=0;i<m_Op_Ext_Lor->m_LM_Count;++i) for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
{ {
curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i]; curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i];
curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]); curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);

40
FDTD/extensions/engine_ext_mur_abc.cpp
View File

@ -41,7 +41,7 @@ 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... //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; int maxDelay=-1;
for (unsigned int n=0;n<m_Op_mur->m_Op->Exc->Volt_Count;++n) for (unsigned int n=0; n<m_Op_mur->m_Op->Exc->Volt_Count; ++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 ( ((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) )
{ {
@ -75,7 +75,7 @@ void Engine_Ext_Mur_ABC::SetNumberOfThreads(int nrThread)
m_numX = AssignJobs2Threads(m_numLines[0],m_NrThreads,false); m_numX = AssignJobs2Threads(m_numLines[0],m_NrThreads,false);
m_start.resize(m_NrThreads,0); m_start.resize(m_NrThreads,0);
m_start.at(0)=0; m_start.at(0)=0;
for (size_t n=1;n<m_numX.size();++n) for (size_t n=1; n<m_numX.size(); ++n)
m_start.at(n) = m_start.at(n-1) + m_numX.at(n-1); m_start.at(n) = m_start.at(n-1) + m_numX.at(n-1);
} }
@ -96,11 +96,11 @@ void Engine_Ext_Mur_ABC::DoPreVoltageUpdates(int threadID)
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
pos_shift[m_nyP] = pos[m_nyP]; pos_shift[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
pos_shift[m_nyPP] = pos[m_nyPP]; pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->Engine::GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyP,pos); m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->Engine::GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyP,pos);
@ -112,10 +112,10 @@ void Engine_Ext_Mur_ABC::DoPreVoltageUpdates(int threadID)
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*) m_Eng; Engine_sse* eng_sse = (Engine_sse*) m_Eng;
for (pos[m_nyP]=0;pos[m_nyP]<m_numLines[0];++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_numLines[0]; ++pos[m_nyP])
{ {
pos_shift[m_nyP] = pos[m_nyP]; pos_shift[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
pos_shift[m_nyPP] = pos[m_nyPP]; pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = eng_sse->Engine_sse::GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyP,pos); m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = eng_sse->Engine_sse::GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyP,pos);
@ -125,11 +125,11 @@ void Engine_Ext_Mur_ABC::DoPreVoltageUpdates(int threadID)
break; break;
} }
default: default:
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
pos_shift[m_nyP] = pos[m_nyP]; pos_shift[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
pos_shift[m_nyPP] = pos[m_nyPP]; pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos); m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] = m_Eng->GetVolt(m_nyP,pos_shift) - m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos);
@ -156,11 +156,11 @@ void Engine_Ext_Mur_ABC::DoPostVoltageUpdates(int threadID)
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
pos_shift[m_nyP] = pos[m_nyP]; pos_shift[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
pos_shift[m_nyPP] = pos[m_nyPP]; pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyP,pos_shift); m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->Engine::GetVolt(m_nyP,pos_shift);
@ -173,11 +173,11 @@ void Engine_Ext_Mur_ABC::DoPostVoltageUpdates(int threadID)
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*) m_Eng; Engine_sse* eng_sse = (Engine_sse*) m_Eng;
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
pos_shift[m_nyP] = pos[m_nyP]; pos_shift[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
pos_shift[m_nyPP] = pos[m_nyPP]; pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyP,pos_shift); m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * eng_sse->Engine_sse::GetVolt(m_nyP,pos_shift);
@ -188,11 +188,11 @@ void Engine_Ext_Mur_ABC::DoPostVoltageUpdates(int threadID)
} }
default: default:
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
pos_shift[m_nyP] = pos[m_nyP]; pos_shift[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
pos_shift[m_nyPP] = pos[m_nyPP]; pos_shift[m_nyPP] = pos[m_nyPP];
m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos_shift); m_volt_nyP[pos[m_nyP]][pos[m_nyPP]] += m_Op_mur->m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] * m_Eng->GetVolt(m_nyP,pos_shift);
@ -217,10 +217,10 @@ void Engine_Ext_Mur_ABC::Apply2Voltages(int threadID)
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
m_Eng->Engine::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); m_Eng->Engine::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]);
m_Eng->Engine::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); m_Eng->Engine::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]);
@ -232,10 +232,10 @@ void Engine_Ext_Mur_ABC::Apply2Voltages(int threadID)
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*) m_Eng; Engine_sse* eng_sse = (Engine_sse*) m_Eng;
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
eng_sse->Engine_sse::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); eng_sse->Engine_sse::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]);
eng_sse->Engine_sse::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); eng_sse->Engine_sse::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]);
@ -245,10 +245,10 @@ void Engine_Ext_Mur_ABC::Apply2Voltages(int threadID)
} }
default: default:
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
pos[m_nyP]=lineX+m_start.at(threadID); pos[m_nyP]=lineX+m_start.at(threadID);
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
m_Eng->SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); m_Eng->SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]);
m_Eng->SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); m_Eng->SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]);

32
FDTD/extensions/engine_ext_pml_sf.cpp
View File

@ -50,13 +50,13 @@ void Engine_Ext_PML_SF::UpdateVoltages(unsigned int startX, unsigned int numX)
pos[0] = startX; pos[0] = startX;
//voltage updates //voltage updates
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
shift[0]=pos[0]; shift[0]=pos[0];
for (pos[1]=0;pos[1]<m_Op_PML_SF->m_numLines[1];++pos[1]) for (pos[1]=0; pos[1]<m_Op_PML_SF->m_numLines[1]; ++pos[1])
{ {
shift[1]=pos[1]; shift[1]=pos[1];
for (pos[2]=0;pos[2]<m_Op_PML_SF->m_numLines[2];++pos[2]) for (pos[2]=0; pos[2]<m_Op_PML_SF->m_numLines[2]; ++pos[2])
{ {
shift[2]=pos[2]; shift[2]=pos[2];
//do the updates here //do the updates here
@ -90,11 +90,11 @@ void Engine_Ext_PML_SF::UpdateCurrents(unsigned int startX, unsigned int numX)
{ {
unsigned int pos[3]; unsigned int pos[3];
pos[0] = startX; pos[0] = startX;
for (unsigned int posX=0;posX<numX;++posX) for (unsigned int posX=0; posX<numX; ++posX)
{ {
for (pos[1]=0;pos[1]<m_Op_PML_SF->m_numLines[1]-1;++pos[1]) for (pos[1]=0; pos[1]<m_Op_PML_SF->m_numLines[1]-1; ++pos[1])
{ {
for (pos[2]=0;pos[2]<m_Op_PML_SF->m_numLines[2]-1;++pos[2]) for (pos[2]=0; pos[2]<m_Op_PML_SF->m_numLines[2]-1; ++pos[2])
{ {
//do the updates here //do the updates here
//for x //for x
@ -150,10 +150,10 @@ void Engine_Ext_PML_SF_Plane::Apply2Voltages()
// copy currents data from main engine to pml engine (lowest main line to highest pml line) // copy currents data from main engine to pml engine (lowest main line to highest pml line)
if (m_Op_PML_SF_PL->m_top==false) if (m_Op_PML_SF_PL->m_top==false)
{ {
for (pos[m_nyP]=0;pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1;++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1; ++pos[m_nyP])
{ {
pml_pos[m_nyP] = pos[m_nyP]; pml_pos[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1;++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1; ++pos[m_nyPP])
{ {
pml_pos[m_nyPP] = pos[m_nyPP]; pml_pos[m_nyPP] = pos[m_nyPP];
curr[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetCurr(0,pos); curr[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetCurr(0,pos);
@ -174,10 +174,10 @@ void Engine_Ext_PML_SF_Plane::Apply2Voltages()
// copy voltage data from pml engine to main engine (highest pml line to lowest main line) // copy voltage data from pml engine to main engine (highest pml line to lowest main line)
pos[m_ny] = 0; pos[m_ny] = 0;
pml_pos[m_ny] = m_Op_PML_SF_PL->m_numLines[m_ny]-1; pml_pos[m_ny] = m_Op_PML_SF_PL->m_numLines[m_ny]-1;
for (pos[m_nyP]=0;pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP];++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]; ++pos[m_nyP])
{ {
pml_pos[m_nyP] = pos[m_nyP]; pml_pos[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]; ++pos[m_nyPP])
{ {
pml_pos[m_nyPP] = pos[m_nyPP]; pml_pos[m_nyPP] = pos[m_nyPP];
m_Eng->SetVolt(0,pos, volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] ); m_Eng->SetVolt(0,pos, volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] );
@ -201,10 +201,10 @@ void Engine_Ext_PML_SF_Plane::Apply2Current()
// copy voltage data from main engine to pml engine (highest main line to lowest pml line) // copy voltage data from main engine to pml engine (highest main line to lowest pml line)
if (m_Op_PML_SF_PL->m_top) if (m_Op_PML_SF_PL->m_top)
{ {
for (pos[m_nyP]=0;pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP];++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]; ++pos[m_nyP])
{ {
pml_pos[m_nyP] = pos[m_nyP]; pml_pos[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]; ++pos[m_nyPP])
{ {
pml_pos[m_nyPP] = pos[m_nyPP]; pml_pos[m_nyPP] = pos[m_nyPP];
volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(0,pos); volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(0,pos);
@ -223,10 +223,10 @@ void Engine_Ext_PML_SF_Plane::Apply2Current()
// this is necessary to catch the voltage excitation on the lowest main voltage line... // this is necessary to catch the voltage excitation on the lowest main voltage line...
if (m_Op_PML_SF_PL->m_top==false) if (m_Op_PML_SF_PL->m_top==false)
{ {
for (pos[m_nyP]=0;pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1;++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1; ++pos[m_nyP])
{ {
pml_pos[m_nyP] = pos[m_nyP]; pml_pos[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1;++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1; ++pos[m_nyPP])
{ {
pml_pos[m_nyPP] = pos[m_nyPP]; pml_pos[m_nyPP] = pos[m_nyPP];
volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(0,pos); volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] = m_Eng->GetVolt(0,pos);
@ -246,10 +246,10 @@ void Engine_Ext_PML_SF_Plane::Apply2Current()
// copy currents data from pml engine to main engine (lowest pml line to highest main line) // copy currents data from pml engine to main engine (lowest pml line to highest main line)
if (m_Op_PML_SF_PL->m_top) if (m_Op_PML_SF_PL->m_top)
{ {
for (pos[m_nyP]=0;pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1;++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_Op_PML_SF_PL->m_numLines[m_nyP]-1; ++pos[m_nyP])
{ {
pml_pos[m_nyP] = pos[m_nyP]; pml_pos[m_nyP] = pos[m_nyP];
for (pos[m_nyPP]=0;pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1;++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_Op_PML_SF_PL->m_numLines[m_nyPP]-1; ++pos[m_nyPP])
{ {
pml_pos[m_nyPP] = pos[m_nyPP]; pml_pos[m_nyPP] = pos[m_nyPP];

80
FDTD/extensions/engine_ext_upml.cpp
View File

@ -50,7 +50,7 @@ void Engine_Ext_UPML::SetNumberOfThreads(int nrThread)
m_numX = AssignJobs2Threads(m_Op_UPML->m_numLines[0],m_NrThreads,false); m_numX = AssignJobs2Threads(m_Op_UPML->m_numLines[0],m_NrThreads,false);
m_start.resize(m_NrThreads,0); m_start.resize(m_NrThreads,0);
m_start.at(0)=0; m_start.at(0)=0;
for (size_t n=1;n<m_numX.size();++n) for (size_t n=1; n<m_numX.size(); ++n)
m_start.at(n) = m_start.at(n-1) + m_numX.at(n-1); m_start.at(n) = m_start.at(n-1) + m_numX.at(n-1);
} }
@ -70,14 +70,14 @@ void Engine_Ext_UPML::DoPreVoltageUpdates(int threadID)
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -103,14 +103,14 @@ void Engine_Ext_UPML::DoPreVoltageUpdates(int threadID)
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*) m_Eng; Engine_sse* eng_sse = (Engine_sse*) m_Eng;
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -135,14 +135,14 @@ void Engine_Ext_UPML::DoPreVoltageUpdates(int threadID)
} }
default: default:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -184,14 +184,14 @@ void Engine_Ext_UPML::DoPostVoltageUpdates(int threadID)
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -214,14 +214,14 @@ void Engine_Ext_UPML::DoPostVoltageUpdates(int threadID)
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*) m_Eng; Engine_sse* eng_sse = (Engine_sse*) m_Eng;
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -243,14 +243,14 @@ void Engine_Ext_UPML::DoPostVoltageUpdates(int threadID)
} }
default: default:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -264,7 +264,8 @@ void Engine_Ext_UPML::DoPostVoltageUpdates(int threadID)
f_help = volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; f_help = volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]];
volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(2,pos); volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetVolt(2,pos);
m_Eng->SetVolt(2,pos, f_help + m_Op_UPML->vvfn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); } m_Eng->SetVolt(2,pos, f_help + m_Op_UPML->vvfn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * volt_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]);
}
} }
} }
break; break;
@ -288,14 +289,14 @@ void Engine_Ext_UPML::DoPreCurrentUpdates(int threadID)
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -321,14 +322,14 @@ void Engine_Ext_UPML::DoPreCurrentUpdates(int threadID)
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*) m_Eng; Engine_sse* eng_sse = (Engine_sse*) m_Eng;
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -354,14 +355,14 @@ void Engine_Ext_UPML::DoPreCurrentUpdates(int threadID)
} }
default: default:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -402,14 +403,14 @@ void Engine_Ext_UPML::DoPostCurrentUpdates(int threadID)
{ {
case Engine::BASIC: case Engine::BASIC:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -432,14 +433,14 @@ void Engine_Ext_UPML::DoPostCurrentUpdates(int threadID)
case Engine::SSE: case Engine::SSE:
{ {
Engine_sse* eng_sse = (Engine_sse*) m_Eng; Engine_sse* eng_sse = (Engine_sse*) m_Eng;
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -461,14 +462,14 @@ void Engine_Ext_UPML::DoPostCurrentUpdates(int threadID)
} }
default: default:
{ {
for (unsigned int lineX=0;lineX<m_numX.at(threadID);++lineX) for (unsigned int lineX=0; lineX<m_numX.at(threadID); ++lineX)
{ {
loc_pos[0]=lineX+m_start.at(threadID); loc_pos[0]=lineX+m_start.at(threadID);
pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0]; pos[0] = loc_pos[0] + m_Op_UPML->m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_Op_UPML->m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_Op_UPML->m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1]; pos[1] = loc_pos[1] + m_Op_UPML->m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_Op_UPML->m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_Op_UPML->m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2]; pos[2] = loc_pos[2] + m_Op_UPML->m_StartPos[2];
@ -482,7 +483,8 @@ void Engine_Ext_UPML::DoPostCurrentUpdates(int threadID)
f_help = curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]; f_help = curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]];
curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(2,pos); curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] = m_Eng->GetCurr(2,pos);
m_Eng->SetCurr(2,pos, f_help + m_Op_UPML->iifn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]); } m_Eng->SetCurr(2,pos, f_help + m_Op_UPML->iifn[2][loc_pos[0]][loc_pos[1]][loc_pos[2]] * curr_flux[2][loc_pos[0]][loc_pos[1]][loc_pos[2]]);
}
} }
} }
break; break;

16
FDTD/extensions/operator_ext_cylinder.cpp
View File

@ -32,14 +32,18 @@ Operator_Ext_Cylinder::Operator_Ext_Cylinder(Operator_Cylinder* op) : Operator_E
Operator_Ext_Cylinder::~Operator_Ext_Cylinder() Operator_Ext_Cylinder::~Operator_Ext_Cylinder()
{ {
delete[] vv_R0;vv_R0=NULL; delete[] vv_R0;
delete[] vi_R0;vi_R0=NULL; vv_R0=NULL;
delete[] vi_R0;
vi_R0=NULL;
} }
bool Operator_Ext_Cylinder::BuildExtension() bool Operator_Ext_Cylinder::BuildExtension()
{ {
delete[] vv_R0;vv_R0=NULL; delete[] vv_R0;
delete[] vi_R0;vi_R0=NULL; vv_R0=NULL;
delete[] vi_R0;
vi_R0=NULL;
//if r=0 is not included -> obviously no special treatment for r=0 //if r=0 is not included -> obviously no special treatment for r=0
//if alpha direction is not closed, PEC-BC at r=0 necessary and already set... //if alpha direction is not closed, PEC-BC at r=0 necessary and already set...
@ -53,11 +57,11 @@ bool Operator_Ext_Cylinder::BuildExtension()
double inEC[4]; double inEC[4];
double dT = m_Op->GetTimestep(); double dT = m_Op->GetTimestep();
pos[0]=0; pos[0]=0;
for (pos[2]=0;pos[2]<m_Op->GetOriginalNumLines(2);++pos[2]) for (pos[2]=0; pos[2]<m_Op->GetOriginalNumLines(2); ++pos[2])
{ {
double C=0; double C=0;
double G=0; double G=0;
for (pos[1]=0;pos[1]<m_Op->GetOriginalNumLines(1)-1;++pos[1]) for (pos[1]=0; pos[1]<m_Op->GetOriginalNumLines(1)-1; ++pos[1])
{ {
m_Op_Cyl->Calc_ECPos(2,pos,inEC); m_Op_Cyl->Calc_ECPos(2,pos,inEC);
C+=inEC[0]*0.5; C+=inEC[0]*0.5;

2
FDTD/extensions/operator_ext_cylinder.h
View File

@ -35,7 +35,7 @@ public:
virtual Engine_Extension* CreateEngineExtention(); virtual Engine_Extension* CreateEngineExtention();
virtual bool IsCylinderCoordsSave() const {return true;} virtual bool IsCylinderCoordsSave() const {return true;}
virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child);return true;} virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;}
virtual std::string GetExtensionName() const {return std::string("Extension for the Cylinder-Coords Operator");} virtual std::string GetExtensionName() const {return std::string("Extension for the Cylinder-Coords Operator");}

24
FDTD/extensions/operator_ext_lorentzmaterial.cpp
View File

@ -20,7 +20,7 @@
Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Operator_Ext_Dispersive(op) Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Operator_Ext_Dispersive(op)
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
v_int_ADE[n] = NULL; v_int_ADE[n] = NULL;
v_ext_ADE[n] = NULL; v_ext_ADE[n] = NULL;
@ -31,7 +31,7 @@ Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Opera
Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial() Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial()
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delete[] v_int_ADE[n]; delete[] v_int_ADE[n];
v_int_ADE[n] = NULL; v_int_ADE[n] = NULL;
@ -61,16 +61,16 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
vector<double> i_ext[3]; vector<double> i_ext[3];
vector<unsigned int> v_pos[3]; vector<unsigned int> v_pos[3];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
unsigned int index = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]); unsigned int index = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]);
//calc epsilon lorentz material //calc epsilon lorentz material
b_pos_on = false; b_pos_on = false;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
L_D[n]=0; L_D[n]=0;
coord[0] = m_Op->GetDiscLine(0,pos[0]); coord[0] = m_Op->GetDiscLine(0,pos[0]);
@ -91,7 +91,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
} }
} }
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
C_D[n]=0; C_D[n]=0;
coord[0] = m_Op->GetDiscLine(0,pos[0],true); coord[0] = m_Op->GetDiscLine(0,pos[0],true);
@ -114,7 +114,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
if (b_pos_on) //this position has active lorentz material if (b_pos_on) //this position has active lorentz material
{ {
for (unsigned int n=0;n<3;++n) for (unsigned int n=0; n<3; ++n)
{ {
v_pos[n].push_back(pos[n]); v_pos[n].push_back(pos[n]);
v_int[n].push_back(1); v_int[n].push_back(1);
@ -137,17 +137,17 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
//copy all vectors into the array's //copy all vectors into the array's
m_LM_Count = v_pos[0].size(); m_LM_Count = v_pos[0].size();
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
m_LM_pos[n] = new unsigned int[m_LM_Count]; m_LM_pos[n] = new unsigned int[m_LM_Count];
for (unsigned int i=0;i<m_LM_Count;++i) for (unsigned int i=0; i<m_LM_Count; ++i)
m_LM_pos[n][i] = v_pos[n].at(i); m_LM_pos[n][i] = v_pos[n].at(i);
if (m_volt_ADE_On) if (m_volt_ADE_On)
{ {
v_int_ADE[n] = new FDTD_FLOAT[m_LM_Count]; v_int_ADE[n] = new FDTD_FLOAT[m_LM_Count];
v_ext_ADE[n] = new FDTD_FLOAT[m_LM_Count]; v_ext_ADE[n] = new FDTD_FLOAT[m_LM_Count];
for (unsigned int i=0;i<m_LM_Count;++i) for (unsigned int i=0; i<m_LM_Count; ++i)
{ {
v_int_ADE[n][i] = v_int[n].at(i); v_int_ADE[n][i] = v_int[n].at(i);
v_ext_ADE[n][i] = v_ext[n].at(i); v_ext_ADE[n][i] = v_ext[n].at(i);
@ -158,7 +158,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
i_int_ADE[n] = new FDTD_FLOAT[m_LM_Count]; i_int_ADE[n] = new FDTD_FLOAT[m_LM_Count];
i_ext_ADE[n] = new FDTD_FLOAT[m_LM_Count]; i_ext_ADE[n] = new FDTD_FLOAT[m_LM_Count];
for (unsigned int i=0;i<m_LM_Count;++i) for (unsigned int i=0; i<m_LM_Count; ++i)
{ {
i_int_ADE[n][i] = i_int[n].at(i); i_int_ADE[n][i] = i_int[n].at(i);
i_ext_ADE[n][i] = i_ext[n].at(i); i_ext_ADE[n][i] = i_ext[n].at(i);

4
FDTD/extensions/operator_ext_mur_abc.cpp
View File

@ -138,10 +138,10 @@ bool Operator_Ext_Mur_ABC::BuildExtension()
else else
coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) - delta/2 / m_Op->GetGridDelta(); coord[m_ny] = m_Op->GetDiscLine(m_ny,pos[m_ny]) - delta/2 / m_Op->GetGridDelta();
for (pos[m_nyP]=0;pos[m_nyP]<m_numLines[0];++pos[m_nyP]) for (pos[m_nyP]=0; pos[m_nyP]<m_numLines[0]; ++pos[m_nyP])
{ {
coord[m_nyP] = m_Op->GetDiscLine(m_nyP,pos[m_nyP]); coord[m_nyP] = m_Op->GetDiscLine(m_nyP,pos[m_nyP]);
for (pos[m_nyPP]=0;pos[m_nyPP]<m_numLines[1];++pos[m_nyPP]) for (pos[m_nyPP]=0; pos[m_nyPP]<m_numLines[1]; ++pos[m_nyPP])
{ {
coord[m_nyPP] = m_Op->GetDiscLine(m_nyPP,pos[m_nyPP]); coord[m_nyPP] = m_Op->GetDiscLine(m_nyPP,pos[m_nyPP]);
CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, CSProperties::MATERIAL, false); CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, CSProperties::MATERIAL, false);

36
FDTD/extensions/operator_ext_pml_sf.cpp
View File

@ -23,7 +23,7 @@
bool Build_Split_Field_PML(Operator* op, int BC[6], int size[6], string gradFunc) bool Build_Split_Field_PML(Operator* op, int BC[6], int size[6], string gradFunc)
{ {
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
{ {
if (BC[n]==3) //split field PML if (BC[n]==3) //split field PML
{ {
@ -61,12 +61,16 @@ Operator_Ext_PML_SF::Operator_Ext_PML_SF(Operator* op) : Operator_Extension(op)
m_numLines[1]=0; m_numLines[1]=0;
m_numLines[2]=0; m_numLines[2]=0;
vv[0] = NULL; vv[1] = NULL; vv[0] = NULL;
vi[0] = NULL; vi[1] = NULL; vv[1] = NULL;
ii[0] = NULL; ii[1] = NULL; vi[0] = NULL;
iv[0] = NULL; iv[1] = NULL; vi[1] = NULL;
ii[0] = NULL;
ii[1] = NULL;
iv[0] = NULL;
iv[1] = NULL;
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
m_BC[n]=0; m_BC[n]=0;
m_GradingFunction = new FunctionParser(); m_GradingFunction = new FunctionParser();
@ -131,7 +135,7 @@ bool Operator_Ext_PML_SF::SetGradingFunction(string func)
m_GradFunc = func; m_GradFunc = func;
int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N"); int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N");
if(res < 0) return true; if (res < 0) return true;
cerr << "Operator_Ext_PML_SF::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl; cerr << "Operator_Ext_PML_SF::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl;
cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n"; cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n";
@ -152,13 +156,13 @@ bool Operator_Ext_PML_SF::BuildExtension()
double inEC[4]; double inEC[4];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
for (pos[0]=0;pos[0]<m_numLines[0];++pos[0]) for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
{ {
for (pos[1]=0;pos[1]<m_numLines[1];++pos[1]) for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<m_numLines[2];++pos[2]) for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
{ {
Calc_ECPos(0,n,pos,inEC); Calc_ECPos(0,n,pos,inEC);
if (inEC[0]>0) if (inEC[0]>0)
@ -334,7 +338,7 @@ void Operator_Ext_PML_SF_Plane::ApplyBC()
if (m_top==false) if (m_top==false)
PEC[2*m_ny+1] = 0; PEC[2*m_ny+1] = 0;
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
{ {
PMC[n] = (m_BC[n] == 1); PMC[n] = (m_BC[n] == 1);
if (n/2 == m_ny) if (n/2 == m_ny)
@ -343,15 +347,15 @@ void Operator_Ext_PML_SF_Plane::ApplyBC()
//apply BC //apply BC
unsigned int pos[3] = {0,0,0}; unsigned int pos[3] = {0,0,0};
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
int nP = (n+1)%3; int nP = (n+1)%3;
int nPP = (n+2)%3; int nPP = (n+2)%3;
for (pos[nP]=0;pos[nP]<m_numLines[nP];++pos[nP]) for (pos[nP]=0; pos[nP]<m_numLines[nP]; ++pos[nP])
{ {
for (pos[nPP]=0;pos[nPP]<m_numLines[nPP];++pos[nPP]) for (pos[nPP]=0; pos[nPP]<m_numLines[nPP]; ++pos[nPP])
{ {
for (int m=0;m<2;++m) for (int m=0; m<2; ++m)
{ {
pos[n]=0; pos[n]=0;
GetVV(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n]; GetVV(m,nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!PEC[2*n];

10
FDTD/extensions/operator_ext_pml_sf.h
View File

@ -34,18 +34,18 @@ class Operator_Ext_PML_SF : public Operator_Extension
public: public:
~Operator_Ext_PML_SF(); ~Operator_Ext_PML_SF();
virtual void SetBoundaryCondition(int* BCs) {for (int n=0;n<6;++n) m_BC[n]=BCs[n];} virtual void SetBoundaryCondition(int* BCs) {for (int n=0; n<6; ++n) m_BC[n]=BCs[n];}
inline virtual FDTD_FLOAT& GetVV(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[nP][n][x][y][z]; } inline virtual FDTD_FLOAT& GetVV(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[nP][n][x][y][z]; }
inline virtual FDTD_FLOAT& GetVI(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[nP][n][x][y][z]; } inline virtual FDTD_FLOAT& GetVI(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[nP][n][x][y][z]; }
inline virtual FDTD_FLOAT& GetII(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[nP][n][x][y][z]; } inline virtual FDTD_FLOAT& GetII(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[nP][n][x][y][z]; }
inline virtual FDTD_FLOAT& GetIV(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[nP][n][x][y][z]; } inline virtual FDTD_FLOAT& GetIV(unsigned int nP, unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[nP][n][x][y][z]; }
virtual double GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetEdgeArea(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
virtual double GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny);UNUSED(pos);UNUSED(dualMesh);return 0.0;} virtual double GetEdgeLength(int ny, unsigned int pos[3], bool dualMesh = false) const {UNUSED(ny); UNUSED(pos); UNUSED(dualMesh); return 0.0;}
//! This will resturn the pml parameter grading //! This will resturn the pml parameter grading
virtual double GetKappaGraded(double depth, double Zm) const {UNUSED(depth);UNUSED(Zm);return 0.0;} virtual double GetKappaGraded(double depth, double Zm) const {UNUSED(depth); UNUSED(Zm); return 0.0;}
//! Set the grading function for the pml //! Set the grading function for the pml
/*! /*!
@ -73,7 +73,7 @@ protected:
virtual void ApplyBC() {}; virtual void ApplyBC() {};
virtual bool Calc_ECPos(int nP, int n, unsigned int* pos, double* inEC) const {UNUSED(n);UNUSED(nP);UNUSED(pos);UNUSED(inEC);return true;}; virtual bool Calc_ECPos(int nP, int n, unsigned int* pos, double* inEC) const {UNUSED(n); UNUSED(nP); UNUSED(pos); UNUSED(inEC); return true;};
unsigned int m_numLines[3]; unsigned int m_numLines[3];
bool m_SetupDone; bool m_SetupDone;

24
FDTD/extensions/operator_ext_upml.cpp
View File

@ -27,12 +27,12 @@ Operator_Ext_UPML::Operator_Ext_UPML(Operator* op) : Operator_Extension(op)
//default grading function //default grading function
SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z "); SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z ");
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
{ {
m_BC[n]=0; m_BC[n]=0;
m_Size[n]=0; m_Size[n]=0;
} }
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
m_StartPos[n]=0; m_StartPos[n]=0;
m_numLines[n]=0; m_numLines[n]=0;
@ -55,7 +55,7 @@ Operator_Ext_UPML::~Operator_Ext_UPML()
void Operator_Ext_UPML::SetBoundaryCondition(const int* BCs, const unsigned int size[6]) void Operator_Ext_UPML::SetBoundaryCondition(const int* BCs, const unsigned int size[6])
{ {
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
{ {
m_BC[n]=BCs[n]; m_BC[n]=BCs[n];
m_Size[n]=size[n]; m_Size[n]=size[n];
@ -64,7 +64,7 @@ void Operator_Ext_UPML::SetBoundaryCondition(const int* BCs, const unsigned int
void Operator_Ext_UPML::SetRange(const unsigned int start[3], const unsigned int stop[3]) void Operator_Ext_UPML::SetRange(const unsigned int start[3], const unsigned int stop[3])
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
m_StartPos[n]=start[n]; m_StartPos[n]=start[n];
m_numLines[n]=stop[n]-start[n]+1; m_numLines[n]=stop[n]-start[n]+1;
@ -77,7 +77,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
unsigned int size[6]={ui_size[0],ui_size[1],ui_size[2],ui_size[3],ui_size[4],ui_size[5]}; unsigned int size[6]={ui_size[0],ui_size[1],ui_size[2],ui_size[3],ui_size[4],ui_size[5]};
//check if mesh is large enough to support the pml //check if mesh is large enough to support the pml
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetOriginalNumLines(n) ) if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetOriginalNumLines(n) )
{ {
cerr << "Operator_Ext_UPML::Create_UPML: Warning: Not enough lines in direction: " << n << ", resetting to PEC" << endl; cerr << "Operator_Ext_UPML::Create_UPML: Warning: Not enough lines in direction: " << n << ", resetting to PEC" << endl;
@ -199,7 +199,7 @@ bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsi
{ {
Operator_Cylinder* op_child = op_cyl_MG->GetInnerOperator(); Operator_Cylinder* op_child = op_cyl_MG->GetInnerOperator();
op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op_child); op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op_child);
for (int n=0;n<2;++n) for (int n=0; n<2; ++n)
{ {
start[n]=0; start[n]=0;
stop[n]=op_child->GetOriginalNumLines(n)-1; stop[n]=op_child->GetOriginalNumLines(n)-1;
@ -258,7 +258,7 @@ bool Operator_Ext_UPML::SetGradingFunction(string func)
m_GradFunc = func; m_GradFunc = func;
int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N"); int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N");
if(res < 0) return true; if (res < 0) return true;
cerr << "Operator_Ext_UPML::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl; cerr << "Operator_Ext_UPML::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl;
cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n"; cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n";
@ -269,7 +269,7 @@ void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm,
{ {
double depth=0; double depth=0;
double width=0; double width=0;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if ((pos[n] <= m_Size[2*n]) && (m_BC[2*n]==3)) //lower pml in n-dir if ((pos[n] <= m_Size[2*n]) && (m_BC[2*n]==3)) //lower pml in n-dir
{ {
@ -355,16 +355,16 @@ bool Operator_Ext_UPML::BuildExtension()
double eff_Mat[4]; double eff_Mat[4];
double dT = m_Op->GetTimestep(); double dT = m_Op->GetTimestep();
for (loc_pos[0]=0;loc_pos[0]<m_numLines[0];++loc_pos[0]) for (loc_pos[0]=0; loc_pos[0]<m_numLines[0]; ++loc_pos[0])
{ {
pos[0] = loc_pos[0] + m_StartPos[0]; pos[0] = loc_pos[0] + m_StartPos[0];
for (loc_pos[1]=0;loc_pos[1]<m_numLines[1];++loc_pos[1]) for (loc_pos[1]=0; loc_pos[1]<m_numLines[1]; ++loc_pos[1])
{ {
pos[1] = loc_pos[1] + m_StartPos[1]; pos[1] = loc_pos[1] + m_StartPos[1];
for (loc_pos[2]=0;loc_pos[2]<m_numLines[2];++loc_pos[2]) for (loc_pos[2]=0; loc_pos[2]<m_numLines[2]; ++loc_pos[2])
{ {
pos[2] = loc_pos[2] + m_StartPos[2]; pos[2] = loc_pos[2] + m_StartPos[2];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
m_Op->Calc_EffMatPos(0,pos,eff_Mat); m_Op->Calc_EffMatPos(0,pos,eff_Mat);
CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i); CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i);

2
FDTD/extensions/operator_extension.h
View File

@ -53,7 +53,7 @@ public:
virtual bool IsCylinderCoordsSave() const {return false;} virtual bool IsCylinderCoordsSave() const {return false;}
//! The cylindrical multi grid operator will check whether the extension is save to use. Default is false. Derive this method to override. //! The cylindrical multi grid operator will check whether the extension is save to use. Default is false. Derive this method to override.
virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child);return false;} virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return false;}
virtual std::string GetExtensionName() const {return std::string("Abstract Operator Extension Base Class");} virtual std::string GetExtensionName() const {return std::string("Abstract Operator Extension Base Class");}

178
FDTD/operator.cpp
View File

@ -39,7 +39,7 @@ Operator::Operator() : Operator_Base()
Operator::~Operator() Operator::~Operator()
{ {
for (size_t n=0;n<m_Op_exts.size();++n) for (size_t n=0; n<m_Op_exts.size(); ++n)
delete m_Op_exts.at(n); delete m_Op_exts.at(n);
m_Op_exts.clear(); m_Op_exts.clear();
Reset(); Reset();
@ -65,7 +65,7 @@ void Operator::Init()
MainOp=NULL; MainOp=NULL;
DualOp=NULL; DualOp=NULL;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
EC_C[n]=NULL; EC_C[n]=NULL;
EC_G[n]=NULL; EC_G[n]=NULL;
@ -84,7 +84,7 @@ void Operator::Reset()
Delete_N_3DArray(ii,numLines); Delete_N_3DArray(ii,numLines);
delete MainOp; delete MainOp;
delete DualOp; delete DualOp;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delete[] EC_C[n]; delete[] EC_C[n];
delete[] EC_G[n]; delete[] EC_G[n];
@ -128,7 +128,7 @@ bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, boo
{ {
bool ok=true; bool ok=true;
unsigned int numLines[3]; unsigned int numLines[3];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
numLines[n] = GetNumberOfLines(n); numLines[n] = GetNumberOfLines(n);
if (inside) //set defaults if (inside) //set defaults
@ -136,7 +136,8 @@ bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, boo
uicoord[n]=0; uicoord[n]=0;
if (dcoord[n]<discLines[n][0]) if (dcoord[n]<discLines[n][0])
{ {
ok=false;uicoord[n]=0; ok=false;
uicoord[n]=0;
if (inside) inside[n] = false; if (inside) inside[n] = false;
} }
else if (dcoord[n]==discLines[n][0]) else if (dcoord[n]==discLines[n][0])
@ -154,7 +155,7 @@ bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower, boo
if (lower) uicoord[n]=numLines[n]-2; if (lower) uicoord[n]=numLines[n]-2;
} }
else else
for (unsigned int i=1;i<numLines[n];++i) for (unsigned int i=1; i<numLines[n]; ++i)
{ {
if (dcoord[n]<discLines[n][i]) if (dcoord[n]<discLines[n][i])
{ {
@ -201,7 +202,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
while (minFoot<stopFoot) while (minFoot<stopFoot)
{ {
minDist=1e300; minDist=1e300;
for (int n=0;n<3;++n) //check all 6 surrounding points for (int n=0; n<3; ++n) //check all 6 surrounding points
{ {
P[0] = discLines[0][currPos[0]]; P[0] = discLines[0][currPos[0]];
P[1] = discLines[1][currPos[1]]; P[1] = discLines[1][currPos[1]];
@ -251,7 +252,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
} }
//close missing edges, if currPos is not equal to uiStopPos //close missing edges, if currPos is not equal to uiStopPos
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if (currPos[n]>uiStop[n]) if (currPos[n]>uiStop[n])
{ {
@ -315,7 +316,7 @@ void Operator::ShowExtStat() const
{ {
if (m_Op_exts.size()==0) return; if (m_Op_exts.size()==0) return;
cout << "-----------------------------------" << endl; cout << "-----------------------------------" << endl;
for (size_t n=0;n<m_Op_exts.size();++n) for (size_t n=0; n<m_Op_exts.size(); ++n)
m_Op_exts.at(n)->ShowStat(cout); m_Op_exts.at(n)->ShowStat(cout);
cout << "-----------------------------------" << endl; cout << "-----------------------------------" << endl;
} }
@ -340,7 +341,7 @@ void Operator::DumpOperator2File(string filename)
FDTD_FLOAT**** exc = Create_N_3DArray<FDTD_FLOAT>(numLines); FDTD_FLOAT**** exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
if (Exc) if (Exc)
{ {
for (unsigned int n=0;n<Exc->Volt_Count;++n) for (unsigned int n=0; n<Exc->Volt_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]; exc[Exc->Volt_dir[n]][Exc->Volt_index[0][n]][Exc->Volt_index[1][n]][Exc->Volt_index[2][n]] = Exc->Volt_amp[n];
} }
@ -400,9 +401,12 @@ void Operator::DumpPEC2File( string filename )
double scaling = 1; double scaling = 1;
#endif #endif
for (pos[0]=0; pos[0]<numLines[0]-1; pos[0]++) { for (pos[0]=0; pos[0]<numLines[0]-1; pos[0]++)
for (pos[1]=0; pos[1]<numLines[1]-1; pos[1]++) { {
for (pos[2]=0; pos[2]<numLines[2]-1; pos[2]++) { for (pos[1]=0; pos[1]<numLines[1]-1; pos[1]++)
{
for (pos[2]=0; pos[2]<numLines[2]-1; pos[2]++)
{
if ((pos[1] != 0) && (pos[2] != 0)) if ((pos[1] != 0) && (pos[2] != 0))
{ {
// PEC surrounds the computational area; do not output this // PEC surrounds the computational area; do not output this
@ -484,13 +488,13 @@ void Operator::DumpMaterial2File(string filename)
FDTD_FLOAT**** sigma = Create_N_3DArray<FDTD_FLOAT>(numLines); FDTD_FLOAT**** sigma = Create_N_3DArray<FDTD_FLOAT>(numLines);
unsigned int pos[3]; unsigned int pos[3];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
double inMat[4]; double inMat[4];
Calc_EffMatPos(n, pos, inMat); Calc_EffMatPos(n, pos, inMat);
@ -525,15 +529,25 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo)
CSX = geo; CSX = geo;
CSRectGrid* grid=CSX->GetGrid(); CSRectGrid* grid=CSX->GetGrid();
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true); discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
if (n==1) if (n==1)
if (numLines[n]<3) {cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl; Reset(); return false;} if (numLines[n]<3)
{
cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl;
Reset();
return false;
}
} }
MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]); MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
MainOp->SetGrid(discLines[0],discLines[1],discLines[2]); MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
if (grid->GetDeltaUnit()<=0) {cerr << "CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!!" << endl; Reset(); return false;} if (grid->GetDeltaUnit()<=0)
{
cerr << "CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!!" << endl;
Reset();
return false;
}
else gridDelta=grid->GetDeltaUnit(); else gridDelta=grid->GetDeltaUnit();
MainOp->SetGridDelta(1); MainOp->SetGridDelta(1);
MainOp->AddCellAdrOp(); MainOp->AddCellAdrOp();
@ -612,13 +626,13 @@ int Operator::CalcECOperator( DebugFlags debugFlags )
unsigned int pos[3]; unsigned int pos[3];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
Calc_ECOperatorPos(n,pos); Calc_ECOperatorPos(n,pos);
} }
@ -629,7 +643,7 @@ int Operator::CalcECOperator( DebugFlags debugFlags )
//Apply PEC to all boundary's //Apply PEC to all boundary's
bool PEC[6]={1,1,1,1,1,1}; bool PEC[6]={1,1,1,1,1,1};
//make an exception for BC == -1 //make an exception for BC == -1
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
if ((m_BC[n]==-1)) if ((m_BC[n]==-1))
PEC[n] = false; PEC[n] = false;
ApplyElectricBC(PEC); ApplyElectricBC(PEC);
@ -637,7 +651,7 @@ int Operator::CalcECOperator( DebugFlags debugFlags )
CalcPEC(); CalcPEC();
bool PMC[6]; bool PMC[6];
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
PMC[n] = m_BC[n]==1; PMC[n] = m_BC[n]==1;
ApplyMagneticBC(PMC); ApplyMagneticBC(PMC);
@ -646,7 +660,7 @@ int Operator::CalcECOperator( DebugFlags debugFlags )
if (CalcFieldExcitation()==false) return -1; if (CalcFieldExcitation()==false) return -1;
//all information available for extension... create now... //all information available for extension... create now...
for (size_t n=0;n<m_Op_exts.size();++n) for (size_t n=0; n<m_Op_exts.size(); ++n)
m_Op_exts.at(n)->BuildExtension(); m_Op_exts.at(n)->BuildExtension();
if (debugFlags & debugMaterial) if (debugFlags & debugMaterial)
@ -657,12 +671,16 @@ int Operator::CalcECOperator( DebugFlags debugFlags )
DumpPEC2File( "PEC_dump.vtk" ); DumpPEC2File( "PEC_dump.vtk" );
//cleanup //cleanup
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
delete[] EC_C[n];EC_C[n]=NULL; delete[] EC_C[n];
delete[] EC_G[n];EC_G[n]=NULL; EC_C[n]=NULL;
delete[] EC_L[n];EC_L[n]=NULL; delete[] EC_G[n];
delete[] EC_R[n];EC_R[n]=NULL; EC_G[n]=NULL;
delete[] EC_L[n];
EC_L[n]=NULL;
delete[] EC_R[n];
EC_R[n]=NULL;
} }
return 0; return 0;
@ -674,13 +692,13 @@ void Operator::ApplyElectricBC(bool* dirs)
return; return;
unsigned int pos[3]; unsigned int pos[3];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
int nP = (n+1)%3; int nP = (n+1)%3;
int nPP = (n+2)%3; int nPP = (n+2)%3;
for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP]) for (pos[nP]=0; pos[nP]<numLines[nP]; ++pos[nP])
{ {
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP]) for (pos[nPP]=0; pos[nPP]<numLines[nPP]; ++pos[nPP])
{ {
if (dirs[2*n]) if (dirs[2*n])
{ {
@ -715,13 +733,13 @@ void Operator::ApplyMagneticBC(bool* dirs)
return; return;
unsigned int pos[3]; unsigned int pos[3];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
int nP = (n+1)%3; int nP = (n+1)%3;
int nPP = (n+2)%3; int nPP = (n+2)%3;
for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP]) for (pos[nP]=0; pos[nP]<numLines[nP]; ++pos[nP])
{ {
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP]) for (pos[nPP]=0; pos[nPP]<numLines[nPP]; ++pos[nPP])
{ {
if (dirs[2*n]) if (dirs[2*n])
{ {
@ -912,7 +930,9 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
EffMat[1]/=area; EffMat[1]/=area;
//******************************* mu,sigma averaging *****************************// //******************************* mu,sigma averaging *****************************//
loc_pos[0]=pos[0];loc_pos[1]=pos[1];loc_pos[2]=pos[2]; loc_pos[0]=pos[0];
loc_pos[1]=pos[1];
loc_pos[2]=pos[2];
double length=0; double length=0;
//shift down //shift down
shiftCoord[n] = coord[n]-delta_M*0.25; shiftCoord[n] = coord[n]-delta_M*0.25;
@ -963,7 +983,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
EffMat[2] = length * __MUE0__ / EffMat[2]; EffMat[2] = length * __MUE0__ / EffMat[2];
if (EffMat[3]) EffMat[3]=length / EffMat[3]; if (EffMat[3]) EffMat[3]=length / EffMat[3];
for (int n=0;n<4;++n) for (int n=0; n<4; ++n)
if (isnan(EffMat[n]) || isinf(EffMat[n])) if (isnan(EffMat[n]) || isinf(EffMat[n]))
{ {
cerr << "Operator::Calc_EffMatPos: An effective material parameter is not a valid result, this should NOT have happend... exit..." << endl; cerr << "Operator::Calc_EffMatPos: An effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
@ -975,7 +995,7 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat) c
void Operator::Init_EC() void Operator::Init_EC()
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
//init x-cell-array //init x-cell-array
delete[] EC_C[n]; delete[] EC_C[n];
@ -986,7 +1006,7 @@ void Operator::Init_EC()
EC_G[n] = new double[MainOp->GetSize()]; EC_G[n] = new double[MainOp->GetSize()];
EC_L[n] = new double[MainOp->GetSize()]; EC_L[n] = new double[MainOp->GetSize()];
EC_R[n] = new double[MainOp->GetSize()]; EC_R[n] = new double[MainOp->GetSize()];
for (unsigned int i=0;i<MainOp->GetSize();i++) //init all for (unsigned int i=0; i<MainOp->GetSize(); i++) //init all
{ {
EC_C[n][i]=0; EC_C[n][i]=0;
EC_G[n][i]=0; EC_G[n][i]=0;
@ -998,17 +1018,21 @@ void Operator::Init_EC()
bool Operator::Calc_EC() bool Operator::Calc_EC()
{ {
if (CSX==NULL) {cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl; return false;} if (CSX==NULL)
{
cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl;
return false;
}
unsigned int ipos; unsigned int ipos;
unsigned int pos[3]; unsigned int pos[3];
double inEC[4]; double inEC[4];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
Calc_ECPos(n,pos,inEC); Calc_ECPos(n,pos,inEC);
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
@ -1046,16 +1070,16 @@ double Operator::CalcTimestep_Var1()
unsigned int ipos_PM; unsigned int ipos_PM;
unsigned int ipos_PPM; unsigned int ipos_PPM;
MainOp->SetReflection2Cell(); MainOp->SetReflection2Cell();
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
int nP = (n+1)%3; int nP = (n+1)%3;
int nPP = (n+2)%3; int nPP = (n+2)%3;
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
ipos_PM = MainOp->Shift(nP,-1); ipos_PM = MainOp->Shift(nP,-1);
@ -1082,7 +1106,7 @@ double min(double* val, unsigned int count)
if (count==0) if (count==0)
return 0.0; return 0.0;
double min = val[0]; double min = val[0];
for (unsigned int n=1;n<count;++n) for (unsigned int n=1; n<count; ++n)
if (val[n]<min) if (val[n]<min)
min = val[n]; min = val[n];
return min; return min;
@ -1101,16 +1125,16 @@ double Operator::CalcTimestep_Var3()
double wqp=0,wt1=0,wt2=0; double wqp=0,wt1=0,wt2=0;
double wt_4[4]={0,0,0,0}; double wt_4[4]={0,0,0,0};
MainOp->SetReflection2Cell(); MainOp->SetReflection2Cell();
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
int nP = (n+1)%3; int nP = (n+1)%3;
int nPP = (n+2)%3; int nPP = (n+2)%3;
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
MainOp->ResetShift(); MainOp->ResetShift();
ipos = MainOp->SetPos(pos[0],pos[1],pos[2]); ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
@ -1191,24 +1215,24 @@ bool Operator::CalcFieldExcitation()
CSProperties* prop=NULL; CSProperties* prop=NULL;
int priority=0; int priority=0;
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
delta[2]=fabs(MainOp->GetIndexDelta(2,pos[2])); delta[2]=fabs(MainOp->GetIndexDelta(2,pos[2]));
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
delta[1]=fabs(MainOp->GetIndexDelta(1,pos[1])); delta[1]=fabs(MainOp->GetIndexDelta(1,pos[1]));
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
delta[0]=fabs(MainOp->GetIndexDelta(0,pos[0])); delta[0]=fabs(MainOp->GetIndexDelta(0,pos[0]));
//electric field excite //electric field excite
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
volt_coord[0] = discLines[0][pos[0]]; volt_coord[0] = discLines[0][pos[0]];
volt_coord[1] = discLines[1][pos[1]]; volt_coord[1] = discLines[1][pos[1]];
volt_coord[2] = discLines[2][pos[2]]; volt_coord[2] = discLines[2][pos[2]];
volt_coord[n]+=delta[n]*0.5; volt_coord[n]+=delta[n]*0.5;
for (size_t p=0;p<vec_prop.size();++p) for (size_t p=0; p<vec_prop.size(); ++p)
{ {
prop = vec_prop.at(p); prop = vec_prop.at(p);
elec = prop->ToElectrode(); elec = prop->ToElectrode();
@ -1239,7 +1263,7 @@ bool Operator::CalcFieldExcitation()
} }
//magnetic field excite //magnetic field excite
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
if ((pos[0]>=numLines[0]-1) || (pos[1]>=numLines[1]-1) || (pos[2]>=numLines[2]-1)) 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 continue; //skip the last H-Line which is outside the FDTD-domain
@ -1250,7 +1274,7 @@ bool Operator::CalcFieldExcitation()
curr_coord[2] = discLines[2][pos[2]]; curr_coord[2] = discLines[2][pos[2]];
curr_coord[nP] +=delta[nP]*0.5; curr_coord[nP] +=delta[nP]*0.5;
curr_coord[nPP] +=delta[nPP]*0.5; curr_coord[nPP] +=delta[nPP]*0.5;
for (size_t p=0;p<vec_prop.size();++p) for (size_t p=0; p<vec_prop.size(); ++p)
{ {
prop = vec_prop.at(p); prop = vec_prop.at(p);
elec = prop->ToElectrode(); elec = prop->ToElectrode();
@ -1289,24 +1313,24 @@ bool Operator::CalcFieldExcitation()
double p2[3]; double p2[3];
double deltaN=0.0; double deltaN=0.0;
struct Grid_Path path; struct Grid_Path path;
for (size_t p=0;p<vec_prop.size();++p) for (size_t p=0; p<vec_prop.size(); ++p)
{ {
prop = vec_prop.at(p); prop = vec_prop.at(p);
elec = prop->ToElectrode(); elec = prop->ToElectrode();
for (size_t n=0;n<prop->GetQtyPrimitives();++n) for (size_t n=0; n<prop->GetQtyPrimitives(); ++n)
{ {
CSPrimitives* prim = prop->GetPrimitive(n); CSPrimitives* prim = prop->GetPrimitive(n);
CSPrimCurve* curv = prim->ToCurve(); CSPrimCurve* curv = prim->ToCurve();
if (curv) if (curv)
{ {
for (size_t i=1;i<curv->GetNumberOfPoints();++i) for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
{ {
curv->GetPoint(i-1,p1); curv->GetPoint(i-1,p1);
curv->GetPoint(i,p2); curv->GetPoint(i,p2);
path = FindPath(p1,p2); path = FindPath(p1,p2);
if (path.dir.size()>0) if (path.dir.size()>0)
prim->SetPrimitiveUsed(true); prim->SetPrimitiveUsed(true);
for (size_t t=0;t<path.dir.size();++t) for (size_t t=0; t<path.dir.size(); ++t)
{ {
n = path.dir.at(t); n = path.dir.at(t);
pos[0] = path.posPath[0].at(t); pos[0] = path.posPath[0].at(t);
@ -1358,7 +1382,9 @@ bool Operator::CalcFieldExcitation()
bool Operator::CalcPEC() bool Operator::CalcPEC()
{ {
m_Nr_PEC[0]=0; m_Nr_PEC[1]=0; m_Nr_PEC[2]=0; m_Nr_PEC[0]=0;
m_Nr_PEC[1]=0;
m_Nr_PEC[2]=0;
CalcPEC_Range(0,numLines[0]-1,m_Nr_PEC); CalcPEC_Range(0,numLines[0]-1,m_Nr_PEC);
@ -1372,13 +1398,13 @@ void Operator::CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned i
double coord[3]; double coord[3];
double delta; double delta;
unsigned int pos[3]; unsigned int pos[3];
for (pos[0]=startX;pos[0]<=stopX;++pos[0]) for (pos[0]=startX; pos[0]<=stopX; ++pos[0])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
coord[0] = discLines[0][pos[0]]; coord[0] = discLines[0][pos[0]];
coord[1] = discLines[1][pos[1]]; coord[1] = discLines[1][pos[1]];
@ -1409,23 +1435,23 @@ void Operator::CalcPEC_Curves()
double p2[3]; double p2[3];
struct Grid_Path path; struct Grid_Path path;
vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::METAL); vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::METAL);
for (size_t p=0;p<vec_prop.size();++p) for (size_t p=0; p<vec_prop.size(); ++p)
{ {
CSProperties* prop = vec_prop.at(p); CSProperties* prop = vec_prop.at(p);
for (size_t n=0;n<prop->GetQtyPrimitives();++n) for (size_t n=0; n<prop->GetQtyPrimitives(); ++n)
{ {
CSPrimitives* prim = prop->GetPrimitive(n); CSPrimitives* prim = prop->GetPrimitive(n);
CSPrimCurve* curv = prim->ToCurve(); CSPrimCurve* curv = prim->ToCurve();
if (curv) if (curv)
{ {
for (size_t i=1;i<curv->GetNumberOfPoints();++i) for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
{ {
curv->GetPoint(i-1,p1); curv->GetPoint(i-1,p1);
curv->GetPoint(i,p2); curv->GetPoint(i,p2);
path = FindPath(p1,p2); path = FindPath(p1,p2);
if (path.dir.size()>0) if (path.dir.size()>0)
prim->SetPrimitiveUsed(true); prim->SetPrimitiveUsed(true);
for (size_t t=0;t<path.dir.size();++t) for (size_t t=0; t<path.dir.size(); ++t)
{ {
// cerr << path.dir.at(t) << " " << path.posPath[0].at(t) << " " << path.posPath[1].at(t) << " " << path.posPath[2].at(t) << endl; // cerr << path.dir.at(t) << " " << path.posPath[0].at(t) << " " << path.posPath[1].at(t) << " " << path.posPath[2].at(t) << endl;
SetVV(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t), 0 ); SetVV(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t), 0 );

22
FDTD/operator_cylinder.cpp
View File

@ -183,14 +183,22 @@ bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo)
} }
} }
else if (minmaxA>2*PI) else if (minmaxA>2*PI)
{cerr << "Operator_Cylinder::SetGeometryCSX: Alpha Max-Min must not be larger than 2*PI!!!" << endl; Reset(); return false;} {
cerr << "Operator_Cylinder::SetGeometryCSX: Alpha Max-Min must not be larger than 2*PI!!!" << endl;
Reset();
return false;
}
else else
{ {
CC_closedAlpha=false; CC_closedAlpha=false;
} }
if (discLines[0][0]<0) if (discLines[0][0]<0)
{cerr << "Operator_Cylinder::SetGeometryCSX: r<0 not allowed in Cylinder Coordinates!!!" << endl; Reset(); return false;} {
cerr << "Operator_Cylinder::SetGeometryCSX: r<0 not allowed in Cylinder Coordinates!!!" << endl;
Reset();
return false;
}
else if (discLines[0][0]==0.0) else if (discLines[0][0]==0.0)
{ {
cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl; cout << "Operator_Cylinder::SetGeometryCSX: r=0 included..." << endl;
@ -208,15 +216,16 @@ void Operator_Cylinder::ApplyElectricBC(bool* dirs)
if (dirs==NULL) return; if (dirs==NULL) return;
if (CC_closedAlpha) if (CC_closedAlpha)
{ {
dirs[2]=0;dirs[3]=0; //no PEC in alpha directions... dirs[2]=0;
dirs[3]=0; //no PEC in alpha directions...
} }
if (CC_R0_included) if (CC_R0_included)
{ {
// E in alpha direction ( aka volt[1][x][y][z] ) is not defined for r==0 --> always zero... // E in alpha direction ( aka volt[1][x][y][z] ) is not defined for r==0 --> always zero...
unsigned int pos[3] = {0,0,0}; unsigned int pos[3] = {0,0,0};
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
SetVV(1,pos[0],pos[1],pos[2], 0 ); SetVV(1,pos[0],pos[1],pos[2], 0 );
SetVI(1,pos[0],pos[1],pos[2], 0 ); SetVI(1,pos[0],pos[1],pos[2], 0 );
@ -231,7 +240,8 @@ void Operator_Cylinder::ApplyMagneticBC(bool* dirs)
if (dirs==NULL) return; if (dirs==NULL) return;
if (CC_closedAlpha) if (CC_closedAlpha)
{ {
dirs[2]=0;dirs[3]=0; //no PMC in alpha directions... dirs[2]=0;
dirs[3]=0; //no PMC in alpha directions...
} }
if (CC_R0_included) if (CC_R0_included)
{ {

6
FDTD/operator_cylindermultigrid.cpp
View File

@ -67,7 +67,7 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
//check if mesh is homogenous in alpha-direction //check if mesh is homogenous in alpha-direction
double diff=discLines[1][1]-discLines[1][0]; double diff=discLines[1][1]-discLines[1][0];
for (unsigned int n=2;n<numLines[1];++n) for (unsigned int n=2; n<numLines[1]; ++n)
{ {
if ( fabs((discLines[1][n]-discLines[1][n-1]) - diff)/diff > 1e-10) if ( fabs((discLines[1][n]-discLines[1][n-1]) - diff)/diff > 1e-10)
{ {
@ -77,7 +77,7 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
} }
m_Split_Pos = 0; m_Split_Pos = 0;
for (unsigned int n=0;n<numLines[0];++n) for (unsigned int n=0; n<numLines[0]; ++n)
{ {
if (m_Split_Rad < discLines[0][n]) if (m_Split_Rad < discLines[0][n])
{ {
@ -97,7 +97,7 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo)
grid->ClearLines(0); grid->ClearLines(0);
grid->ClearLines(1); grid->ClearLines(1);
for (unsigned int n=0; n<m_Split_Pos ;++n) for (unsigned int n=0; n<m_Split_Pos ; ++n)
grid->AddDiscLine(0,discLines[0][n]); grid->AddDiscLine(0,discLines[0][n]);
for (unsigned int n=0; n<numLines[1]; n+=2) for (unsigned int n=0; n<numLines[1]; n+=2)
grid->AddDiscLine(1,discLines[1][n]); grid->AddDiscLine(1,discLines[1][n]);

48
FDTD/operator_multithread.cpp
View File

@ -68,11 +68,15 @@ void Operator_Multithread::Reset()
m_thread_group.join_all(); m_thread_group.join_all();
delete m_CalcEC_Start;m_CalcEC_Start=NULL; delete m_CalcEC_Start;
delete m_CalcEC_Stop;m_CalcEC_Stop=NULL; m_CalcEC_Start=NULL;
delete m_CalcEC_Stop;
m_CalcEC_Stop=NULL;
delete m_CalcPEC_Start;m_CalcPEC_Start=NULL; delete m_CalcPEC_Start;
delete m_CalcPEC_Stop;m_CalcPEC_Stop=NULL; m_CalcPEC_Start=NULL;
delete m_CalcPEC_Stop;
m_CalcPEC_Stop=NULL;
} }
void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector<unsigned int> &start, vector<unsigned int> &stop) const
@ -106,11 +110,15 @@ int Operator_Multithread::CalcECOperator( DebugFlags debugFlags )
cout << "Multithreaded operator using " << m_numThreads << " threads." << std::endl; cout << "Multithreaded operator using " << m_numThreads << " threads." << std::endl;
m_thread_group.join_all(); m_thread_group.join_all();
delete m_CalcEC_Start;m_CalcEC_Start = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller delete m_CalcEC_Start;
delete m_CalcEC_Stop;m_CalcEC_Stop = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller m_CalcEC_Start = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
delete m_CalcEC_Stop;
m_CalcEC_Stop = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
delete m_CalcPEC_Start;m_CalcPEC_Start = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller delete m_CalcPEC_Start;
delete m_CalcPEC_Stop;m_CalcPEC_Stop = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller m_CalcPEC_Start = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
delete m_CalcPEC_Stop;
m_CalcPEC_Stop = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
for (unsigned int n=0; n<m_numThreads; n++) for (unsigned int n=0; n<m_numThreads; n++)
{ {
@ -123,7 +131,11 @@ int Operator_Multithread::CalcECOperator( DebugFlags debugFlags )
bool Operator_Multithread::Calc_EC() bool Operator_Multithread::Calc_EC()
{ {
if (CSX==NULL) {cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl; return false;} if (CSX==NULL)
{
cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl;
return false;
}
MainOp->SetPos(0,0,0); MainOp->SetPos(0,0,0);
@ -136,7 +148,9 @@ bool Operator_Multithread::Calc_EC()
bool Operator_Multithread::CalcPEC() bool Operator_Multithread::CalcPEC()
{ {
m_Nr_PEC[0]=0; m_Nr_PEC[1]=0; m_Nr_PEC[2]=0; m_Nr_PEC[0]=0;
m_Nr_PEC[1]=0;
m_Nr_PEC[2]=0;
m_Nr_PEC_thread = new unsigned int[m_numThreads][3]; m_Nr_PEC_thread = new unsigned int[m_numThreads][3];
@ -144,8 +158,8 @@ bool Operator_Multithread::CalcPEC()
m_CalcPEC_Stop->wait(); m_CalcPEC_Stop->wait();
for (unsigned int t=0;t<m_numThreads;++t) for (unsigned int t=0; t<m_numThreads; ++t)
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
m_Nr_PEC[n]+=m_Nr_PEC_thread[t][n]; m_Nr_PEC[n]+=m_Nr_PEC_thread[t][n];
CalcPEC_Curves(); CalcPEC_Curves();
@ -171,13 +185,13 @@ void Operator_Thread::operator()()
unsigned int ipos; unsigned int ipos;
unsigned int pos[3]; unsigned int pos[3];
double inEC[4]; double inEC[4];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
for (pos[0]=m_start;pos[0]<=m_stop;++pos[0]) for (pos[0]=m_start; pos[0]<=m_stop; ++pos[0])
{ {
for (pos[1]=0;pos[1]<m_OpPtr->numLines[1];++pos[1]) for (pos[1]=0; pos[1]<m_OpPtr->numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<m_OpPtr->numLines[2];++pos[2]) for (pos[2]=0; pos[2]<m_OpPtr->numLines[2]; ++pos[2])
{ {
m_OpPtr->Calc_ECPos(n,pos,inEC); m_OpPtr->Calc_ECPos(n,pos,inEC);
ipos = m_OpPtr->MainOp->GetPos(pos[0],pos[1],pos[2]); ipos = m_OpPtr->MainOp->GetPos(pos[0],pos[1],pos[2]);
@ -193,7 +207,7 @@ void Operator_Thread::operator()()
//************** calculate EC (Calc_EC) ***********************// //************** calculate EC (Calc_EC) ***********************//
m_OpPtr->m_CalcPEC_Start->wait(); m_OpPtr->m_CalcPEC_Start->wait();
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
m_OpPtr->m_Nr_PEC_thread[m_threadID][n] = 0; m_OpPtr->m_Nr_PEC_thread[m_threadID][n] = 0;
m_OpPtr->CalcPEC_Range(m_start,m_stop,m_OpPtr->m_Nr_PEC_thread[m_threadID]); m_OpPtr->CalcPEC_Range(m_start,m_stop,m_OpPtr->m_Nr_PEC_thread[m_threadID]);

6
FDTD/operator_sse_compressed.cpp
View File

@ -120,13 +120,13 @@ void Operator_SSE_Compressed::ShowStat() const
INLINE int equal(f4vector v1, f4vector v2) INLINE int equal(f4vector v1, f4vector v2)
{ {
#if defined(__SSE__) #if defined(__SSE__)
#if (__GNUC__ == 4) && (__GNUC_MINOR__ < 4) #if (__GNUC__ == 4) && (__GNUC_MINOR__ < 4)
v4si compare = __builtin_ia32_cmpeqps( v1.v, v2.v ); v4si compare = __builtin_ia32_cmpeqps( v1.v, v2.v );
return __builtin_ia32_movmskps( (v4sf)compare ) == 0x0f; return __builtin_ia32_movmskps( (v4sf)compare ) == 0x0f;
#else #else
v4sf compare = __builtin_ia32_cmpeqps( v1.v, v2.v ); v4sf compare = __builtin_ia32_cmpeqps( v1.v, v2.v );
return __builtin_ia32_movmskps( compare ) == 0x0f; return __builtin_ia32_movmskps( compare ) == 0x0f;
#endif #endif
#else #else
return ( return (
v1.f[0] == v2.f[0] && v1.f[0] == v2.f[0] &&

2
main.cpp
View File

@ -70,7 +70,7 @@ int main(int argc, char *argv[])
if (argc>=3) if (argc>=3)
{ {
for (int n=2;n<argc;++n) for (int n=2; n<argc; ++n)
{ {
if (g_settings.parseCommandLineArgument(argv[n])) if (g_settings.parseCommandLineArgument(argv[n]))
continue; continue;

42
openems.cpp
View File

@ -77,9 +77,12 @@ openEMS::~openEMS()
void openEMS::Reset() void openEMS::Reset()
{ {
if (PA) PA->DeleteAll(); if (PA) PA->DeleteAll();
delete PA; PA=0; delete PA;
delete FDTD_Eng; FDTD_Eng=0; PA=0;
delete FDTD_Op; FDTD_Op=0; delete FDTD_Eng;
FDTD_Eng=0;
delete FDTD_Op;
FDTD_Op=0;
} }
//! \brief processes a command line argument //! \brief processes a command line argument
@ -203,7 +206,7 @@ bool openEMS::SetupBoundaryConditions(TiXmlElement* BC)
string bound_names[] = {"xmin","xmax","ymin","ymax","zmin","zmax"}; string bound_names[] = {"xmin","xmax","ymin","ymax","zmin","zmax"};
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
{ {
EC = BC->QueryIntAttribute(bound_names[n].c_str(),&bounds[n]); EC = BC->QueryIntAttribute(bound_names[n].c_str(),&bounds[n]);
if (EC==TIXML_SUCCESS) if (EC==TIXML_SUCCESS)
@ -219,7 +222,7 @@ bool openEMS::SetupBoundaryConditions(TiXmlElement* BC)
bounds[n] = 1; bounds[n] = 1;
else if (s_bc=="MUR") else if (s_bc=="MUR")
bounds[n] = 2; bounds[n] = 2;
else if(strncmp(s_bc.c_str(),"PML_=",4)==0) else if (strncmp(s_bc.c_str(),"PML_=",4)==0)
{ {
bounds[n] = 3; bounds[n] = 3;
pml_size[n] = atoi(s_bc.c_str()+4); pml_size[n] = atoi(s_bc.c_str()+4);
@ -240,7 +243,7 @@ bool openEMS::SetupBoundaryConditions(TiXmlElement* BC)
if (BC->QueryDoubleAttribute("MUR_PhaseVelocity",&mur_v_ph) != TIXML_SUCCESS) if (BC->QueryDoubleAttribute("MUR_PhaseVelocity",&mur_v_ph) != TIXML_SUCCESS)
mur_v_ph = -1; mur_v_ph = -1;
string mur_v_ph_names[6] = {"MUR_PhaseVelocity_xmin", "MUR_PhaseVelocity_xmax", "MUR_PhaseVelocity_ymin", "MUR_PhaseVelocity_ymax", "MUR_PhaseVelocity_zmin", "MUR_PhaseVelocity_zmax"}; string mur_v_ph_names[6] = {"MUR_PhaseVelocity_xmin", "MUR_PhaseVelocity_xmax", "MUR_PhaseVelocity_ymin", "MUR_PhaseVelocity_ymax", "MUR_PhaseVelocity_zmin", "MUR_PhaseVelocity_zmax"};
for (int n=0;n<6;++n) for (int n=0; n<6; ++n)
{ {
if (bounds[n]==2) //Mur-ABC if (bounds[n]==2) //Mur-ABC
{ {
@ -412,7 +415,8 @@ int openEMS::SetupFDTD(const char* file)
cout << "Creation time for operator: " << CalcDiffTime(OpDoneTime,startTime) << " s" << endl; cout << "Creation time for operator: " << CalcDiffTime(OpDoneTime,startTime) << " s" << endl;
if (m_no_simulation) { if (m_no_simulation)
{
// simulation was disabled (to generate debug output only) // simulation was disabled (to generate debug output only)
return 1; return 1;
} }
@ -429,7 +433,7 @@ int openEMS::SetupFDTD(const char* file)
double start[3]; double start[3];
double stop[3]; double stop[3];
vector<CSProperties*> Probes = CSX.GetPropertyByType(CSProperties::PROBEBOX); vector<CSProperties*> Probes = CSX.GetPropertyByType(CSProperties::PROBEBOX);
for (size_t i=0;i<Probes.size();++i) for (size_t i=0; i<Probes.size(); ++i)
{ {
//only looking for one prim atm //only looking for one prim atm
CSPrimitives* prim = Probes.at(i)->GetPrimitive(0); CSPrimitives* prim = Probes.at(i)->GetPrimitive(0);
@ -438,8 +442,12 @@ int openEMS::SetupFDTD(const char* file)
bool acc; bool acc;
double bnd[6] = {0,0,0,0,0,0}; double bnd[6] = {0,0,0,0,0,0};
acc = prim->GetBoundBox(bnd,true); acc = prim->GetBoundBox(bnd,true);
start[0]= bnd[0];start[1]=bnd[2];start[2]=bnd[4]; start[0]= bnd[0];
stop[0] = bnd[1];stop[1] =bnd[3];stop[2] =bnd[5]; start[1]=bnd[2];
start[2]=bnd[4];
stop[0] = bnd[1];
stop[1] =bnd[3];
stop[2] =bnd[5];
CSPropProbeBox* pb = Probes.at(i)->ToProbeBox(); CSPropProbeBox* pb = Probes.at(i)->ToProbeBox();
Processing* proc = NULL; Processing* proc = NULL;
if (pb) if (pb)
@ -490,7 +498,7 @@ int openEMS::SetupFDTD(const char* file)
} }
vector<CSProperties*> DumpProps = CSX.GetPropertyByType(CSProperties::DUMPBOX); vector<CSProperties*> DumpProps = CSX.GetPropertyByType(CSProperties::DUMPBOX);
for (size_t i=0;i<DumpProps.size();++i) for (size_t i=0; i<DumpProps.size(); ++i)
{ {
ProcessFieldsTD* ProcTD = new ProcessFieldsTD(FDTD_Op); ProcessFieldsTD* ProcTD = new ProcessFieldsTD(FDTD_Op);
ProcTD->SetEnable(Enable_Dumps); ProcTD->SetEnable(Enable_Dumps);
@ -505,8 +513,12 @@ int openEMS::SetupFDTD(const char* file)
bool acc; bool acc;
double bnd[6] = {0,0,0,0,0,0}; double bnd[6] = {0,0,0,0,0,0};
acc = prim->GetBoundBox(bnd,true); acc = prim->GetBoundBox(bnd,true);
start[0]= bnd[0];start[1]=bnd[2];start[2]=bnd[4]; start[0]= bnd[0];
stop[0] = bnd[1];stop[1] =bnd[3];stop[2] =bnd[5]; start[1]=bnd[2];
start[2]=bnd[4];
stop[0] = bnd[1];
stop[1] =bnd[3];
stop[2] =bnd[5];
CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox(); CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox();
if (db) if (db)
{ {
@ -516,7 +528,7 @@ int openEMS::SetupFDTD(const char* file)
ProcTD->SetFileType((ProcessFields::FileType)db->GetFileType()); ProcTD->SetFileType((ProcessFields::FileType)db->GetFileType());
if (CylinderCoords) if (CylinderCoords)
ProcTD->SetMeshType(Processing::CYLINDRICAL_MESH); ProcTD->SetMeshType(Processing::CYLINDRICAL_MESH);
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
ProcTD->SetSubSampling(db->GetSubSampling(n),n); ProcTD->SetSubSampling(db->GetSubSampling(n),n);
ProcTD->SetFilePattern(db->GetName()); ProcTD->SetFilePattern(db->GetName());
ProcTD->SetFileName(db->GetName()); ProcTD->SetFileName(db->GetName());
@ -587,7 +599,7 @@ void openEMS::RunFDTD()
//add all timesteps to end-crit field processing with max excite amplitude //add all timesteps to end-crit field processing with max excite amplitude
unsigned int maxExcite = FDTD_Op->Exc->GetMaxExcitationTimestep(); unsigned int maxExcite = FDTD_Op->Exc->GetMaxExcitationTimestep();
for (unsigned int n=0;n<FDTD_Op->Exc->Volt_Count;++n) for (unsigned int n=0; n<FDTD_Op->Exc->Volt_Count; ++n)
ProcField->AddStep(FDTD_Op->Exc->Volt_delay[n]+maxExcite); ProcField->AddStep(FDTD_Op->Exc->Volt_delay[n]+maxExcite);
double change=1; double change=1;

70
tools/AdrOp.cpp
View File

@ -20,7 +20,12 @@
AdrOp::AdrOp(unsigned int muiImax, unsigned int muiJmax, unsigned int muiKmax, unsigned int muiLmax) AdrOp::AdrOp(unsigned int muiImax, unsigned int muiJmax, unsigned int muiKmax, unsigned int muiLmax)
{ {
//error-handling... //error-handling...
error = new ErrorMsg(9); if (error==NULL) { fprintf(stderr,"Memory allocation failed!! exiting..."); exit(1); } error = new ErrorMsg(9);
if (error==NULL)
{
fprintf(stderr,"Memory allocation failed!! exiting...");
exit(1);
}
error->SetMsg(1,"Adress Operator: Memory allocation failed!! exiting..."); error->SetMsg(1,"Adress Operator: Memory allocation failed!! exiting...");
error->SetMsg(2,"Adress Operator: Invalid Adress requested!! exiting..."); error->SetMsg(2,"Adress Operator: Invalid Adress requested!! exiting...");
error->SetMsg(3,"Adress Operator: Invalid Position set!! exiting..."); error->SetMsg(3,"Adress Operator: Invalid Position set!! exiting...");
@ -58,7 +63,10 @@ AdrOp::AdrOp(unsigned int muiImax, unsigned int muiJmax, unsigned int muiKmax, u
reflect=false; reflect=false;
uiTypeOffset=0; uiTypeOffset=0;
clCellAdr=NULL; clCellAdr=NULL;
dGrid[0]=NULL;dGrid[1]=NULL;dGrid[2]=NULL;dGrid[3]=NULL; dGrid[0]=NULL;
dGrid[1]=NULL;
dGrid[2]=NULL;
dGrid[3]=NULL;
dDeltaUnit=1; dDeltaUnit=1;
bDebug=false; bDebug=false;
} }
@ -78,12 +86,12 @@ AdrOp::AdrOp(AdrOp* origOP)
uiJpos=origOP->uiJpos; uiJpos=origOP->uiJpos;
uiKpos=origOP->uiKpos; uiKpos=origOP->uiKpos;
uiLpos=origOP->uiLpos; uiLpos=origOP->uiLpos;
for (int ii=0;ii<4;++ii) dGrid[ii]=origOP->dGrid[ii]; for (int ii=0; ii<4; ++ii) dGrid[ii]=origOP->dGrid[ii];
dDeltaUnit=origOP->dDeltaUnit; dDeltaUnit=origOP->dDeltaUnit;
iIshift=origOP->iIshift; iIshift=origOP->iIshift;
iJshift=origOP->iJshift; iJshift=origOP->iJshift;
iKshift=origOP->iKshift; iKshift=origOP->iKshift;
for (int ii=0;ii<3;++ii) iCellShift[ii]=origOP->iCellShift[ii]; for (int ii=0; ii<3; ++ii) iCellShift[ii]=origOP->iCellShift[ii];
i=origOP->i; i=origOP->i;
j=origOP->j; j=origOP->j;
k=origOP->k; k=origOP->k;
@ -100,8 +108,10 @@ AdrOp::AdrOp(AdrOp* origOP)
AdrOp::~AdrOp() AdrOp::~AdrOp()
{ {
// cerr << "\n------Adress Operator deconstructed-----\n" << endl; // cerr << "\n------Adress Operator deconstructed-----\n" << endl;
delete error; error=NULL; delete error;
delete clCellAdr; clCellAdr=NULL; error=NULL;
delete clCellAdr;
clCellAdr=NULL;
} }
unsigned int AdrOp::SetPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos) unsigned int AdrOp::SetPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos)
@ -173,14 +183,14 @@ unsigned int AdrOp::GetPos(int muiIrel, int muiJrel, int muiKrel, int /*muiLrel*
if (bPosSet==false) error->Error(6); if (bPosSet==false) error->Error(6);
if (reflect) if (reflect)
{ {
#if EXPENSE_LOG==1 #if EXPENSE_LOG==1
if (muiIrel+(int)uiIpos<0) ADRESSEXPENSE(2,1,0,0,1,0) if (muiIrel+(int)uiIpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
if (muiIrel+(int)uiIpos>(int)uiImax-1) ADRESSEXPENSE(4,1,0,0,1,0) if (muiIrel+(int)uiIpos>(int)uiImax-1) ADRESSEXPENSE(4,1,0,0,1,0)
if (muiJrel+(int)uiJpos<0) ADRESSEXPENSE(2,1,0,0,1,0) if (muiJrel+(int)uiJpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
if (muiJrel+(int)uiJpos>(int)uiJmax-1) ADRESSEXPENSE(4,1,0,0,1,0) if (muiJrel+(int)uiJpos>(int)uiJmax-1) ADRESSEXPENSE(4,1,0,0,1,0)
if (muiKrel+(int)uiKpos<0) ADRESSEXPENSE(2,1,0,0,1,0) if (muiKrel+(int)uiKpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
if (muiKrel+(int)uiKpos>(int)uiKmax-1) ADRESSEXPENSE(4,1,0,0,1,0) if (muiKrel+(int)uiKpos>(int)uiKmax-1) ADRESSEXPENSE(4,1,0,0,1,0)
#endif #endif
if (muiIrel+(int)uiIpos<0) muiIrel=-2*uiIpos-muiIrel-uiTypeOffset; if (muiIrel+(int)uiIpos<0) muiIrel=-2*uiIpos-muiIrel-uiTypeOffset;
if (muiIrel+(int)uiIpos>(int)uiImax-1) muiIrel=2*(uiImax-1-uiIpos)-muiIrel+uiTypeOffset; if (muiIrel+(int)uiIpos>(int)uiImax-1) muiIrel=2*(uiImax-1-uiIpos)-muiIrel+uiTypeOffset;
@ -253,22 +263,34 @@ unsigned int AdrOp::GetPosFromNode(int ny, unsigned int uiNode)
double AdrOp::GetNodeVolume(unsigned int uiNode) double AdrOp::GetNodeVolume(unsigned int uiNode)
{ {
for (unsigned int n=0;n<uiDimension;n++) if (dGrid[n]==NULL) error->Error(9); for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
double dVol=1; double dVol=1;
unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax}; unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
unsigned int uiPos[4]={GetPosFromNode(0,uiNode),GetPosFromNode(1,uiNode),GetPosFromNode(2,uiNode),GetPosFromNode(3,uiNode)}; unsigned int uiPos[4]={GetPosFromNode(0,uiNode),GetPosFromNode(1,uiNode),GetPosFromNode(2,uiNode),GetPosFromNode(3,uiNode)};
for (unsigned int n=0;n<uiDimension;n++) for (unsigned int n=0; n<uiDimension; n++)
{ {
if ((uiPos[n]>0) && (uiPos[n]<uiMax[n]-1)) { dVol*=0.5*dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]-1]); ADRESSEXPENSE(4,0,1,3,0,4) } if ((uiPos[n]>0) && (uiPos[n]<uiMax[n]-1))
else if ((uiPos[n]==0) && (uiPos[n]<uiMax[n]-1)) { dVol*=dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]]); ADRESSEXPENSE(3,0,1,2,0,4) } {
else if ((uiPos[n]>0) && (uiPos[n]==uiMax[n]-1)) { dVol*=dDeltaUnit*(dGrid[n][uiPos[n]]-dGrid[n][uiPos[n]-1]); ADRESSEXPENSE(3,0,1,2,0,4) } dVol*=0.5*dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]-1]);
ADRESSEXPENSE(4,0,1,3,0,4)
}
else if ((uiPos[n]==0) && (uiPos[n]<uiMax[n]-1))
{
dVol*=dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]]);
ADRESSEXPENSE(3,0,1,2,0,4)
}
else if ((uiPos[n]>0) && (uiPos[n]==uiMax[n]-1))
{
dVol*=dDeltaUnit*(dGrid[n][uiPos[n]]-dGrid[n][uiPos[n]-1]);
ADRESSEXPENSE(3,0,1,2,0,4)
}
} }
return dVol; return dVol;
} }
double AdrOp::GetIndexWidth(int ny, int index) double AdrOp::GetIndexWidth(int ny, int index)
{ {
for (unsigned int n=0;n<uiDimension;n++) if (dGrid[n]==NULL) error->Error(9); for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
double width=0; double width=0;
while (ny<0) ny+=uiDimension; while (ny<0) ny+=uiDimension;
ny=ny%uiDimension; ny=ny%uiDimension;
@ -282,7 +304,7 @@ double AdrOp::GetIndexWidth(int ny, int index)
double AdrOp::GetIndexCoord(int ny, int index) double AdrOp::GetIndexCoord(int ny, int index)
{ {
for (unsigned int n=0;n<uiDimension;n++) if (dGrid[n]==NULL) error->Error(9); for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
while (ny<0) ny+=uiDimension; while (ny<0) ny+=uiDimension;
ny=ny%uiDimension; ny=ny%uiDimension;
unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax}; unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
@ -430,7 +452,8 @@ AdrOp* AdrOp::AddCellAdrOp()
AdrOp* AdrOp::DeleteCellAdrOp() AdrOp* AdrOp::DeleteCellAdrOp()
{ {
delete clCellAdr; clCellAdr=NULL; delete clCellAdr;
clCellAdr=NULL;
return NULL; return NULL;
} }
@ -459,9 +482,9 @@ unsigned int AdrOp::GetCellPos(bool incShift)
{ {
if (bPosSet==false) error->Error(6); if (bPosSet==false) error->Error(6);
if (clCellAdr==NULL) error->Error(7); if (clCellAdr==NULL) error->Error(7);
#if EXPENSE_LOG==1 #if EXPENSE_LOG==1
if (incShift) ADRESSEXPENSE(3,0,0,0,0,2) if (incShift) ADRESSEXPENSE(3,0,0,0,0,2)
#endif #endif
if (incShift) return clCellAdr->GetPos(uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]); if (incShift) return clCellAdr->GetPos(uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]);
else return clCellAdr->GetPos(uiIpos,uiJpos,uiKpos); else return clCellAdr->GetPos(uiIpos,uiJpos,uiKpos);
} }
@ -475,14 +498,14 @@ unsigned int AdrOp::GetCellPos(int i, int j, int k)
double AdrOp::GetShiftCellVolume(int ny, int step) double AdrOp::GetShiftCellVolume(int ny, int step)
{ {
for (unsigned int n=0;n<uiDimension;n++) if (dGrid[n]==NULL) error->Error(9); for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
int uiMax[4]={uiImax-1,uiJmax-1,uiKmax-1,uiLmax-1}; int uiMax[4]={uiImax-1,uiJmax-1,uiKmax-1,uiLmax-1};
while (ny<0) ny+=uiDimension; while (ny<0) ny+=uiDimension;
ny=ny%uiDimension; ny=ny%uiDimension;
iCellShift[ny]=step; iCellShift[ny]=step;
int uiPos[4]={uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]}; int uiPos[4]={uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]};
double dVol=1; double dVol=1;
for (unsigned int n=0;n<uiDimension;++n) for (unsigned int n=0; n<uiDimension; ++n)
{ {
if (uiMax[n]>0) if (uiMax[n]>0)
{ {
@ -517,7 +540,12 @@ unsigned int deltaAdrOp::GetAdr(int pos)
if (uiMax==1) return 0; if (uiMax==1) return 0;
if (pos<0) pos=pos*-1; if (pos<0) pos=pos*-1;
else if (pos>(int)uiMax-1) pos=2*(uiMax-1)-pos+1; else if (pos>(int)uiMax-1) pos=2*(uiMax-1)-pos+1;
if ((pos<0) || (pos>(int)uiMax-1)) {fprintf(stderr," Error exiting... "); getchar(); exit(-1);} if ((pos<0) || (pos>(int)uiMax-1))
{
fprintf(stderr," Error exiting... ");
getchar();
exit(-1);
}
return pos; return pos;
} }

3
tools/AdrOp.h
View File

@ -34,7 +34,8 @@
using namespace std; using namespace std;
class AdrOp{ class AdrOp
{
public: public:
///Constructor, define dimension/size here ///Constructor, define dimension/size here
AdrOp(unsigned int muiImax, unsigned int muiYmax, unsigned int muiKmax=0, unsigned int muiLmax=0); AdrOp(unsigned int muiImax, unsigned int muiYmax, unsigned int muiKmax=0, unsigned int muiLmax=0);

25
tools/ErrorMsg.cpp
View File

@ -23,20 +23,35 @@
ErrorMsg::ErrorMsg(unsigned int NoMessage) ErrorMsg::ErrorMsg(unsigned int NoMessage)
{ {
NoMsg=NoMessage; NoMsg=NoMessage;
if (NoMsg>0) Msg = new char*[NoMsg]; if (Msg==NULL) { fprintf(stderr,"Memory allocation failed!! exiting..."); exit(1); } if (NoMsg>0) Msg = new char*[NoMsg];
for (unsigned int i=0;i<NoMsg;i++) Msg[i]=NULL; if (Msg==NULL)
{
fprintf(stderr,"Memory allocation failed!! exiting...");
exit(1);
}
for (unsigned int i=0; i<NoMsg; i++) Msg[i]=NULL;
} }
ErrorMsg::~ErrorMsg() ErrorMsg::~ErrorMsg()
{ {
for (unsigned int i=0;i<NoMsg;i++) {delete[] Msg[i]; Msg[i]=NULL;}; for (unsigned int i=0; i<NoMsg; i++)
delete[] Msg; Msg=NULL; {
delete[] Msg[i];
Msg[i]=NULL;
};
delete[] Msg;
Msg=NULL;
} }
void ErrorMsg::SetMsg(unsigned int nr, const char *Message) void ErrorMsg::SetMsg(unsigned int nr, const char *Message)
{ {
if ((nr<1) || (nr>NoMsg) || (Message==NULL)) ownError(); if ((nr<1) || (nr>NoMsg) || (Message==NULL)) ownError();
Msg[nr-1] = new char[strlen(Message)+1]; if (Msg[nr-1]==NULL) { fprintf(stderr,"Memory allocation failed!! exiting..."); exit(1); } Msg[nr-1] = new char[strlen(Message)+1];
if (Msg[nr-1]==NULL)
{
fprintf(stderr,"Memory allocation failed!! exiting...");
exit(1);
}
Msg[nr-1]=strcpy(Msg[nr-1],Message); Msg[nr-1]=strcpy(Msg[nr-1],Message);
} }

10
tools/ExpenseLog.cpp
View File

@ -104,7 +104,11 @@ ExpenseLog::ExpenseLog(void)
ExpenseLog::~ExpenseLog(void) ExpenseLog::~ExpenseLog(void)
{ {
for (size_t i=0;i<vModules.size();++i) {delete vModules.at(i); vModules.at(i)=NULL;} for (size_t i=0; i<vModules.size(); ++i)
{
delete vModules.at(i);
vModules.at(i)=NULL;
}
vModules.clear(); vModules.clear();
} }
@ -119,7 +123,7 @@ void ExpenseLog::PrintAll(FILE *file)
{ {
double totalAdd=0,totalMul=0,totalBool=0,totalAssign=0; double totalAdd=0,totalMul=0,totalBool=0,totalAssign=0;
fprintf(stderr,"\n ----------------\n Expense Log PrintOut\n Nr of Modules: %d\n",vModules.size()); fprintf(stderr,"\n ----------------\n Expense Log PrintOut\n Nr of Modules: %d\n",vModules.size());
for (size_t i=0;i<vModules.size();++i) for (size_t i=0; i<vModules.size(); ++i)
{ {
totalAdd+=((double)vModules.at(i)->uiIntAdditions+(double)vModules.at(i)->uiDoubleAdditions) + 1e9*((double)vModules.at(i)->uiMrdIA+(double)vModules.at(i)->uiMrdIA); totalAdd+=((double)vModules.at(i)->uiIntAdditions+(double)vModules.at(i)->uiDoubleAdditions) + 1e9*((double)vModules.at(i)->uiMrdIA+(double)vModules.at(i)->uiMrdIA);
totalMul+=((double)vModules.at(i)->uiIntMultiplications+(double)vModules.at(i)->uiDoubleMultiplications) + 1e9*(double)(vModules.at(i)->uiMrdIM+vModules.at(i)->uiMrdIM); totalMul+=((double)vModules.at(i)->uiIntMultiplications+(double)vModules.at(i)->uiDoubleMultiplications) + 1e9*(double)(vModules.at(i)->uiMrdIM+vModules.at(i)->uiMrdIM);
@ -133,7 +137,7 @@ void ExpenseLog::PrintAll(FILE *file)
void ExpenseLog::ClearAll() void ExpenseLog::ClearAll()
{ {
for (size_t i=0;i<vModules.size();++i) for (size_t i=0; i<vModules.size(); ++i)
{ {
vModules.at(i)->Clear(); vModules.at(i)->Clear();
} }

2
tools/ExpenseLog.h
View File

@ -40,7 +40,7 @@ ExpenseModule* AdrOpExpense=EL.AddModule("Adress Operator");
class ExpenseModule class ExpenseModule
{ {
friend class ExpenseLog; friend class ExpenseLog;
public: public:
ExpenseModule(const char* moduleName); ExpenseModule(const char* moduleName);
~ExpenseModule(); ~ExpenseModule();

60
tools/aligned_allocator.h
View File

@ -26,17 +26,18 @@
#ifdef WIN32 #ifdef WIN32
#define __MSVCRT_VERSION__ 0x0700 #define __MSVCRT_VERSION__ 0x0700
#include <malloc.h> #include <malloc.h>
#define MEMALIGN( array, alignment, size ) !(*array = _aligned_malloc( size, alignment )) #define MEMALIGN( array, alignment, size ) !(*array = _aligned_malloc( size, alignment ))
#define FREE( array ) _aligned_free( array ) #define FREE( array ) _aligned_free( array )
#else #else
#define MEMALIGN( array, alignment, size ) posix_memalign( array, alignment, size ) #define MEMALIGN( array, alignment, size ) posix_memalign( array, alignment, size )
#define FREE( array ) free( array ) #define FREE( array ) free( array )
#endif #endif
template <typename T> class aligned_allocator { template <typename T> class aligned_allocator
{
public: public:
// The following will be the same for virtually all allocators. // The following will be the same for virtually all allocators.
typedef T * pointer; typedef T * pointer;
@ -47,30 +48,36 @@ public:
typedef size_t size_type; typedef size_t size_type;
typedef ptrdiff_t difference_type; typedef ptrdiff_t difference_type;
T * address(T& r) const { T * address(T& r) const
{
return &r; return &r;
} }
const T * address(const T& s) const { const T * address(const T& s) const
{
return &s; return &s;
} }
size_t max_size() const { size_t max_size() const
{
// The following has been carefully written to be independent of // The following has been carefully written to be independent of
// the definition of size_t and to avoid signed/unsigned warnings. // the definition of size_t and to avoid signed/unsigned warnings.
return (static_cast<size_t>(0) - static_cast<size_t>(1)) / sizeof(T); return (static_cast<size_t>(0) - static_cast<size_t>(1)) / sizeof(T);
} }
// The following must be the same for all allocators. // The following must be the same for all allocators.
template <typename U> struct rebind { template <typename U> struct rebind
{
typedef aligned_allocator<U> other; typedef aligned_allocator<U> other;
}; };
bool operator!=(const aligned_allocator& other) const { bool operator!=(const aligned_allocator& other) const
{
return !(*this == other); return !(*this == other);
} }
void construct(T * const p, const T& t) const { void construct(T * const p, const T& t) const
{
void * const pv = static_cast<void *>(p); void * const pv = static_cast<void *>(p);
new (pv) T(t); new (pv) T(t);
} }
@ -80,7 +87,8 @@ public:
// Returns true if and only if storage allocated from *this // Returns true if and only if storage allocated from *this
// can be deallocated from other, and vice versa. // can be deallocated from other, and vice versa.
// Always returns true for stateless allocators. // Always returns true for stateless allocators.
bool operator==(const aligned_allocator& other) const { bool operator==(const aligned_allocator& other) const
{
return true; return true;
} }
@ -92,7 +100,8 @@ public:
~aligned_allocator() { } ~aligned_allocator() { }
// The following will be different for each allocator. // The following will be different for each allocator.
T * allocate(const size_t n) const { T * allocate(const size_t n) const
{
// std::cout << "Allocating " << n << (n == 1 ? " object" : "objects") << " of size " << sizeof(T) << "." << std::endl; // std::cout << "Allocating " << n << (n == 1 ? " object" : "objects") << " of size " << sizeof(T) << "." << std::endl;
// The return value of allocate(0) is unspecified. // The return value of allocate(0) is unspecified.
// aligned_allocator returns NULL in order to avoid depending // aligned_allocator returns NULL in order to avoid depending
@ -100,14 +109,16 @@ public:
// (the implementation can define malloc(0) to return NULL, // (the implementation can define malloc(0) to return NULL,
// in which case the bad_alloc check below would fire). // in which case the bad_alloc check below would fire).
// All allocators can return NULL in this case. // All allocators can return NULL in this case.
if (n == 0) { if (n == 0)
{
return NULL; return NULL;
} }
// All allocators should contain an integer overflow check. // All allocators should contain an integer overflow check.
// The Standardization Committee recommends that std::length_error // The Standardization Committee recommends that std::length_error
// be thrown in the case of integer overflow. // be thrown in the case of integer overflow.
if (n > max_size()) { if (n > max_size())
{
throw std::length_error("aligned_allocator<T>::allocate() - Integer overflow."); throw std::length_error("aligned_allocator<T>::allocate() - Integer overflow.");
} }
@ -119,7 +130,8 @@ public:
return static_cast<T *>(pv); return static_cast<T *>(pv);
} }
void deallocate(T * const p, const size_t n) const { void deallocate(T * const p, const size_t n) const
{
// std::cout << "Deallocating " << n << (n == 1 ? " object" : "objects") << " of size " << sizeof(T) << "." << std::endl; // std::cout << "Deallocating " << n << (n == 1 ? " object" : "objects") << " of size " << sizeof(T) << "." << std::endl;
// aligned_allocator wraps free(). // aligned_allocator wraps free().
UNUSED(n); UNUSED(n);
@ -127,7 +139,8 @@ public:
} }
// The following will be the same for all allocators that ignore hints. // The following will be the same for all allocators that ignore hints.
template <typename U> T * allocate(const size_t n, const U * /* const hint */) const { template <typename U> T * allocate(const size_t n, const U * /* const hint */) const
{
return allocate(n); return allocate(n);
} }
@ -145,15 +158,16 @@ private:
// A compiler bug causes it to believe that p->~T() doesn't reference p. // A compiler bug causes it to believe that p->~T() doesn't reference p.
#ifdef _MSC_VER #ifdef _MSC_VER
#pragma warning(push) #pragma warning(push)
#pragma warning(disable: 4100) // unreferenced formal parameter #pragma warning(disable: 4100) // unreferenced formal parameter
#endif #endif
// The definition of destroy() must be the same for all allocators. // The definition of destroy() must be the same for all allocators.
template <typename T> void aligned_allocator<T>::destroy(T * const p) const { template <typename T> void aligned_allocator<T>::destroy(T * const p) const
{
p->~T(); p->~T();
} }
#ifdef _MSC_VER #ifdef _MSC_VER
#pragma warning(pop) #pragma warning(pop)
#endif #endif

43
tools/array_ops.cpp
View File

@ -19,14 +19,14 @@
#include <ostream> #include <ostream>
#ifdef WIN32 #ifdef WIN32
#define __MSVCRT_VERSION__ 0x0700 #define __MSVCRT_VERSION__ 0x0700
#include <malloc.h> #include <malloc.h>
//(void**)&array, 16, sizeof(typeof(f4vector**))*numLines[0] //(void**)&array, 16, sizeof(typeof(f4vector**))*numLines[0]
#define MEMALIGN( array, alignment, size ) !(*array = _aligned_malloc( size, alignment )) #define MEMALIGN( array, alignment, size ) !(*array = _aligned_malloc( size, alignment ))
#define FREE( array ) _aligned_free( array ) #define FREE( array ) _aligned_free( array )
#else #else
#define MEMALIGN( array, alignment, size ) posix_memalign( array, alignment, size ) #define MEMALIGN( array, alignment, size ) posix_memalign( array, alignment, size )
#define FREE( array ) free( array ) #define FREE( array ) free( array )
#endif #endif
void Delete1DArray_v4sf(f4vector* array) void Delete1DArray_v4sf(f4vector* array)
@ -40,9 +40,9 @@ void Delete3DArray_v4sf(f4vector*** array, const unsigned int* numLines)
{ {
if (array==NULL) return; if (array==NULL) return;
unsigned int pos[3]; unsigned int pos[3];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
FREE( array[pos[0]][pos[1]] ); FREE( array[pos[0]][pos[1]] );
//delete[] array[pos[0]][pos[1]]; //delete[] array[pos[0]][pos[1]];
@ -57,7 +57,7 @@ void Delete3DArray_v4sf(f4vector*** array, const unsigned int* numLines)
void Delete_N_3DArray_v4sf(f4vector**** array, const unsigned int* numLines) void Delete_N_3DArray_v4sf(f4vector**** array, const unsigned int* numLines)
{ {
if (array==NULL) return; if (array==NULL) return;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
Delete3DArray_v4sf(array[n],numLines); Delete3DArray_v4sf(array[n],numLines);
} }
@ -68,7 +68,8 @@ void Delete_N_3DArray_v4sf(f4vector**** array, const unsigned int* numLines)
f4vector* Create1DArray_v4sf(const unsigned int numLines) f4vector* Create1DArray_v4sf(const unsigned int numLines)
{ {
f4vector* array=NULL; f4vector* array=NULL;
if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector))*numLines )) { if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector))*numLines ))
{
cerr << "cannot allocate aligned memory" << endl; cerr << "cannot allocate aligned memory" << endl;
exit(3); exit(3);
} }
@ -82,26 +83,29 @@ f4vector*** Create3DArray_v4sf(const unsigned int* numLines)
f4vector*** array=NULL; f4vector*** array=NULL;
unsigned int pos[3]; unsigned int pos[3];
if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector**))*numLines[0] )) { if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector**))*numLines[0] ))
{
cerr << "cannot allocate aligned memory" << endl; cerr << "cannot allocate aligned memory" << endl;
exit(3); exit(3);
} }
//array = new f4vector**[numLines[0]]; //array = new f4vector**[numLines[0]];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{
if (MEMALIGN( (void**)&array[pos[0]], 16, sizeof(typeof(f4vector*))*numLines[1] ))
{ {
if (MEMALIGN( (void**)&array[pos[0]], 16, sizeof(typeof(f4vector*))*numLines[1] )) {
cerr << "cannot allocate aligned memory" << endl; cerr << "cannot allocate aligned memory" << endl;
exit(3); exit(3);
} }
//array[pos[0]] = new f4vector*[numLines[1]]; //array[pos[0]] = new f4vector*[numLines[1]];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{
if (MEMALIGN( (void**)&array[pos[0]][pos[1]], 16, sizeof(typeof(f4vector))*numZ ))
{ {
if (MEMALIGN( (void**)&array[pos[0]][pos[1]], 16, sizeof(typeof(f4vector))*numZ )) {
cerr << "cannot allocate aligned memory" << endl; cerr << "cannot allocate aligned memory" << endl;
exit(3); exit(3);
} }
//array[pos[0]][pos[1]] = new f4vector[numZ]; //array[pos[0]][pos[1]] = new f4vector[numZ];
for (pos[2]=0;pos[2]<numZ;++pos[2]) for (pos[2]=0; pos[2]<numZ; ++pos[2])
{ {
array[pos[0]][pos[1]][pos[2]].f[0] = 0; array[pos[0]][pos[1]][pos[2]].f[0] = 0;
array[pos[0]][pos[1]][pos[2]].f[1] = 0; array[pos[0]][pos[1]][pos[2]].f[1] = 0;
@ -115,12 +119,13 @@ f4vector*** Create3DArray_v4sf(const unsigned int* numLines)
f4vector**** Create_N_3DArray_v4sf(const unsigned int* numLines) f4vector**** Create_N_3DArray_v4sf(const unsigned int* numLines)
{ {
f4vector**** array=NULL; f4vector**** array=NULL;
if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector***))*3 )) { if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector***))*3 ))
{
cerr << "cannot allocate aligned memory" << endl; cerr << "cannot allocate aligned memory" << endl;
exit(3); exit(3);
} }
//array = new f4vector***[3]; //array = new f4vector***[3];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
array[n]=Create3DArray_v4sf(numLines); array[n]=Create3DArray_v4sf(numLines);
} }

30
tools/array_ops.h
View File

@ -20,7 +20,7 @@
#ifdef __SIZEOF_FLOAT__ #ifdef __SIZEOF_FLOAT__
#if __SIZEOF_FLOAT__ != 4 #if __SIZEOF_FLOAT__ != 4
#error wrong size of float #error wrong size of float
#endif #endif
#endif #endif
@ -58,10 +58,10 @@ T** Create2DArray(const unsigned int* numLines)
T** array=NULL; T** array=NULL;
unsigned int pos[3]; unsigned int pos[3];
array = new T*[numLines[0]]; array = new T*[numLines[0]];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
array[pos[0]] = new T[numLines[1]]; array[pos[0]] = new T[numLines[1]];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
array[pos[0]][pos[1]] = 0; array[pos[0]][pos[1]] = 0;
} }
@ -74,7 +74,7 @@ void Delete2DArray(T** array, const unsigned int* numLines)
{ {
if (array==NULL) return; if (array==NULL) return;
unsigned int pos[3]; unsigned int pos[3];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
delete[] array[pos[0]]; delete[] array[pos[0]];
} }
@ -99,13 +99,13 @@ T*** Create3DArray(const unsigned int* numLines)
T*** array=NULL; T*** array=NULL;
unsigned int pos[3]; unsigned int pos[3];
array = new T**[numLines[0]]; array = new T**[numLines[0]];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
array[pos[0]] = new T*[numLines[1]]; array[pos[0]] = new T*[numLines[1]];
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
array[pos[0]][pos[1]] = new T[numLines[2]]; array[pos[0]][pos[1]] = new T[numLines[2]];
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
array[pos[0]][pos[1]][pos[2]] = 0; array[pos[0]][pos[1]][pos[2]] = 0;
} }
@ -119,7 +119,7 @@ T**** Create_N_3DArray(const unsigned int* numLines)
{ {
T**** array=NULL; T**** array=NULL;
array = new T***[3]; array = new T***[3];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
array[n]=Create3DArray<T>( numLines ); array[n]=Create3DArray<T>( numLines );
} }
@ -131,9 +131,9 @@ void Delete3DArray(T*** array, const unsigned int* numLines)
{ {
if (!array) return; if (!array) return;
unsigned int pos[3]; unsigned int pos[3];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
delete[] array[pos[0]][pos[1]]; delete[] array[pos[0]][pos[1]];
} }
@ -146,7 +146,7 @@ template <typename T>
void Delete_N_3DArray(T**** array, const unsigned int* numLines) void Delete_N_3DArray(T**** array, const unsigned int* numLines)
{ {
if (!array) return; if (!array) return;
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
{ {
Delete3DArray<T>(array[n],numLines); Delete3DArray<T>(array[n],numLines);
} }
@ -157,14 +157,14 @@ template <typename T>
void Dump_N_3DArray2File(ostream &file, T**** array, const unsigned int* numLines) void Dump_N_3DArray2File(ostream &file, T**** array, const unsigned int* numLines)
{ {
unsigned int pos[3]; unsigned int pos[3];
for (pos[0]=0;pos[0]<numLines[0];++pos[0]) for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
{ {
for (pos[1]=0;pos[1]<numLines[1];++pos[1]) for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
{ {
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
{ {
file << pos[0] << "\t" << pos[1] << "\t" << pos[2]; file << pos[0] << "\t" << pos[1] << "\t" << pos[2];
for (int n=0;n<3;++n) for (int n=0; n<3; ++n)
file << "\t" << (float)array[n][pos[0]][pos[1]][pos[2]]; file << "\t" << (float)array[n][pos[0]][pos[1]][pos[2]];
file << endl; file << endl;
} }

4
tools/useful.cpp
View File

@ -37,13 +37,13 @@ std::vector<unsigned int> AssignJobs2Threads(unsigned int jobs, unsigned int nrT
std::vector<unsigned int> jpt; //jobs per thread std::vector<unsigned int> jpt; //jobs per thread
unsigned int ui_jpt = jobs/nrThreads; unsigned int ui_jpt = jobs/nrThreads;
for (unsigned int n=0;n<nrThreads;++n) for (unsigned int n=0; n<nrThreads; ++n)
{ {
jpt.push_back(ui_jpt); jpt.push_back(ui_jpt);
jobs-=ui_jpt; jobs-=ui_jpt;
} }
for (unsigned int n=0;n<nrThreads;++n) for (unsigned int n=0; n<nrThreads; ++n)
{ {
if (jobs>0) if (jobs>0)
{ {