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

Merge remote branch 'origin/master' into multithreading

pull/1/head
Sebastian Held 2010-03-29 10:18:33 +02:00
parent b8180287e9 5a47a6443a
commit 09364107cf
10 changed files with 195 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
FDTD/engine.h
View File

@ -22,11 +22,6 @@
class Engine class Engine
{ {
friend class Processing;
friend class ProcessVoltage;
friend class ProcessCurrent;
friend class ProcessFields;
friend class ProcessFieldsTD;
public: public:
static Engine* createEngine(Operator* op); static Engine* createEngine(Operator* op);
virtual ~Engine(); virtual ~Engine();
@ -39,6 +34,9 @@ public:
virtual unsigned int GetNumberOfTimesteps() {return numTS;}; virtual unsigned int GetNumberOfTimesteps() {return numTS;};
virtual FDTD_FLOAT**** GetVoltages() {return volt;};
virtual FDTD_FLOAT**** GetCurrents() {return curr;};
protected: protected:
Engine(Operator* op); Engine(Operator* op);
Operator* Op; Operator* Op;

155
FDTD/operator.cpp
View File

@ -106,7 +106,7 @@ bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower)
else if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;} else if (dcoord[n]>discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
else if (dcoord[n]==discLines[n][numLines[n]-1]) {uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;} else if (dcoord[n]==discLines[n][numLines[n]-1]) {uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
else else
for (unsigned int i=1;i<numLines[n]-1;++i) for (unsigned int i=1;i<numLines[n];++i)
{ {
if (dcoord[n]<discLines[n][i]) if (dcoord[n]<discLines[n][i])
{ {
@ -139,19 +139,17 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
double meshStart[] = {discLines[0][uiStart[0]], discLines[1][uiStart[1]], discLines[2][uiStart[2]]}; double meshStart[] = {discLines[0][uiStart[0]], discLines[1][uiStart[1]], discLines[2][uiStart[2]]};
double meshStop[] = {discLines[0][uiStop[0]], discLines[1][uiStop[1]], discLines[2][uiStop[2]]}; double meshStop[] = {discLines[0][uiStop[0]], discLines[1][uiStop[1]], discLines[2][uiStop[2]]};
double foot,dist,minFoot,minDist,minDir; bool UpDir;
double foot=0,dist=0,minFoot=0,minDist=0;
int minDir=0;
unsigned int minPos[3]; unsigned int minPos[3];
double startFoot,stopFoot,currFoot; double startFoot,stopFoot,currFoot;
Point_Line_Distance(meshStart,start,stop,startFoot,dist); Point_Line_Distance(meshStart,start,stop,startFoot,dist);
Point_Line_Distance(meshStop,start,stop,stopFoot,dist); Point_Line_Distance(meshStop,start,stop,stopFoot,dist);
currFoot=startFoot; currFoot=startFoot;
minFoot=startFoot;
double P[3]; double P[3];
// cerr << "start pos " << discLines[0][currPos[0]] << " " << discLines[1][currPos[1]] << " " << discLines[2][currPos[2]] << endl;
//
// FDTD_FLOAT**** array = Create_N_3DArray(numLines);
while (minFoot<stopFoot) while (minFoot<stopFoot)
{ {
minDist=1e300; minDist=1e300;
@ -160,7 +158,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
P[0] = discLines[0][currPos[0]]; P[0] = discLines[0][currPos[0]];
P[1] = discLines[1][currPos[1]]; P[1] = discLines[1][currPos[1]];
P[2] = discLines[2][currPos[2]]; P[2] = discLines[2][currPos[2]];
if ((currPos[n]-1)>=0) if (((int)currPos[n]-1)>=0)
{ {
P[n] = discLines[n][currPos[n]-1]; P[n] = discLines[n][currPos[n]-1];
Point_Line_Distance(P,start,stop,foot,dist); Point_Line_Distance(P,start,stop,foot,dist);
@ -169,11 +167,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
minFoot=foot; minFoot=foot;
minDist=dist; minDist=dist;
minDir = n; minDir = n;
minPos[0]=currPos[0]; UpDir = false;
minPos[1]=currPos[1];
minPos[2]=currPos[2];
minPos[n]=currPos[n]-1;
// array[n][minPos[0]][minPos[1]][minPos[2]] = 1;
} }
} }
if ((currPos[n]+1)<numLines[n]) if ((currPos[n]+1)<numLines[n])
@ -185,30 +179,28 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
minFoot=foot; minFoot=foot;
minDist=dist; minDist=dist;
minDir = n; minDir = n;
minPos[0]=currPos[0]; UpDir = true;
minPos[1]=currPos[1];
minPos[2]=currPos[2];
minPos[n]=currPos[n]+1;
// array[n][minPos[0]][minPos[1]][minPos[2]] = 1;
} }
} }
} }
// cerr << "next best pos " << minDir << " " << " " << discLines[0][minPos[0]] << " " << discLines[1][minPos[1]] << " " << discLines[2][minPos[2]] << endl; minPos[0]=currPos[0];
currPos[0]=minPos[0]; minPos[1]=currPos[1];
currPos[1]=minPos[1]; minPos[2]=currPos[2];
currPos[2]=minPos[2]; if (UpDir)
currFoot=minFoot; {
path.dir.push_back(minDir); currPos[minDir]+=1;
}
else
{
currPos[minDir]+=-1;
minPos[minDir]-=1;
}
path.posPath[0].push_back(minPos[0]); path.posPath[0].push_back(minPos[0]);
path.posPath[1].push_back(minPos[1]); path.posPath[1].push_back(minPos[1]);
path.posPath[2].push_back(minPos[2]); path.posPath[2].push_back(minPos[2]);
currFoot=minFoot;
path.dir.push_back(minDir);
} }
// ofstream file("test.vtk",ios_base::out);
//
// ProcessFields::DumpVectorArray2VTK(file,"path",array,discLines,numLines);
// file.close();
return path; return path;
} }
@ -232,9 +224,9 @@ void Operator::ShowSize()
cout << "-----------------------------" << endl; cout << "-----------------------------" << endl;
} }
void Operator::CalcGaussianPulsExcitation(double f0, double fc) bool Operator::CalcGaussianPulsExcitation(double f0, double fc)
{ {
if (dT==0) return; if (dT==0) return false;
ExciteLength = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT); ExciteLength = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl; cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
@ -246,14 +238,16 @@ void Operator::CalcGaussianPulsExcitation(double f0, double fc)
ExciteSignal[n] = cos(2.0*PI*f0*(n*dT-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*n*dT/3.0-3,2)); ExciteSignal[n] = cos(2.0*PI*f0*(n*dT-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*n*dT/3.0-3,2));
// cerr << ExciteSignal[n] << endl; // cerr << ExciteSignal[n] << endl;
} }
return true;
} }
void Operator::CalcSinusExcitation(double f0, int nTS) bool Operator::CalcSinusExcitation(double f0, int nTS)
{ {
if (dT==0) return; if (dT==0) return false;
if (nTS<=0) return; if (nTS<=0) return false;
ExciteLength = (unsigned int)(nTS); ExciteLength = (unsigned int)(nTS);
cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] ExciteSignal; delete[] ExciteSignal;
ExciteSignal = new FDTD_FLOAT[ExciteLength+1]; ExciteSignal = new FDTD_FLOAT[ExciteLength+1];
ExciteSignal[0]=0.0; ExciteSignal[0]=0.0;
@ -262,6 +256,7 @@ void Operator::CalcSinusExcitation(double f0, int nTS)
ExciteSignal[n] = sin(2.0*PI*f0*n*dT); ExciteSignal[n] = sin(2.0*PI*f0*n*dT);
// cerr << ExciteSignal[n] << endl; // cerr << ExciteSignal[n] << endl;
} }
return true;
} }
void Operator::DumpOperator2File(string filename) void Operator::DumpOperator2File(string filename)
@ -274,10 +269,18 @@ void Operator::DumpOperator2File(string filename)
return; return;
} }
string names[] = {"vv", "vi", "iv" , "ii"}; FDTD_FLOAT**** exc = Create_N_3DArray(numLines);
FDTD_FLOAT**** array[] = {vv,vi,iv,ii}; for (unsigned int n=0;n<E_Exc_Count;++n)
{
exc[E_Exc_dir[n]][E_Exc_index[0][n]][E_Exc_index[1][n]][E_Exc_index[2][n]] = E_Exc_amp[n];
}
ProcessFields::DumpMultiVectorArray2VTK(file, names , array , 4, discLines, numLines); string names[] = {"vv", "vi", "iv" , "ii", "exc"};
FDTD_FLOAT**** array[] = {vv,vi,iv,ii,exc};
ProcessFields::DumpMultiVectorArray2VTK(file, names , array , 5, discLines, numLines);
Delete_N_3DArray(exc,numLines);
file.close(); file.close();
} }
@ -728,7 +731,7 @@ bool Operator::CalcEFieldExcitation()
vector<unsigned int> vDelay; vector<unsigned int> vDelay;
vector<unsigned int> vDir; vector<unsigned int> vDir;
unsigned int ipos; unsigned int ipos;
int pos[3]; unsigned int pos[3];
double coord[3]; double coord[3];
double delta[3]; double delta[3];
double amp=0; double amp=0;
@ -782,6 +785,70 @@ bool Operator::CalcEFieldExcitation()
} }
} }
//special treatment for primitives of type curve (treated as wires) see also Calc_PEC
double p1[3];
double p2[3];
double deltaN=0.0;
int n;
struct Grid_Path path;
CSPropElectrode* elec=NULL;
CSProperties* prop=NULL;
vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::ELECTRODE);
for (size_t p=0;p<vec_prop.size();++p)
{
prop = vec_prop.at(p);
elec = prop->ToElectrode();
for (size_t n=0;n<prop->GetQtyPrimitives();++n)
{
CSPrimitives* prim = prop->GetPrimitive(n);
CSPrimCurve* curv = prim->ToCurve();
if (curv)
{
for (size_t i=1;i<curv->GetNumberOfPoints();++i)
{
curv->GetPoint(i-1,p1);
curv->GetPoint(i,p2);
path = FindPath(p1,p2);
for (size_t t=0;t<path.dir.size();++t)
{
n = path.dir.at(t);
pos[0] = path.posPath[0].at(t);
pos[1] = path.posPath[1].at(t);
pos[2] = path.posPath[2].at(t);
MainOp->SetPos(pos[0],pos[1],pos[2]);
deltaN=fabs(MainOp->GetIndexDelta(n,pos[n]));
coord[0] = discLines[0][pos[0]];
coord[1] = discLines[1][pos[1]];
coord[2] = discLines[2][pos[2]];
coord[n] += 0.5*deltaN;
// cerr << n << " " << coord[0] << " " << coord[1] << " " << coord[2] << endl;
if (elec!=NULL)
{
if ((elec->GetActiveDir(n)) && (pos[n]<numLines[n]-1))
{
amp = elec->GetWeightedExcitation(n,coord)*deltaN*gridDelta;
if (amp!=0)
{
vExcit.push_back(amp);
vDelay.push_back((unsigned int)(elec->GetDelay()/dT));
vDir.push_back(n);
vIndex[0].push_back(pos[0]);
vIndex[1].push_back(pos[1]);
vIndex[2].push_back(pos[2]);
}
if (elec->GetExcitType()==1) //hard excite
{
vv[n][pos[0]][pos[1]][pos[2]] = 0;
vi[n][pos[0]][pos[1]][pos[2]] = 0;
}
}
}
}
}
}
}
}
E_Exc_Count = vExcit.size(); E_Exc_Count = vExcit.size();
cerr << "Operator::CalcEFieldExcitation: Found number of excitations points: " << E_Exc_Count << endl; cerr << "Operator::CalcEFieldExcitation: Found number of excitations points: " << E_Exc_Count << endl;
if (E_Exc_Count==0) if (E_Exc_Count==0)
@ -862,24 +929,18 @@ bool Operator::CalcPEC()
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);
// cerr << p1[0] << " " << p1[1] << " " << p1[2] << endl;
// cerr << p2[0] << " " << p2[1] << " " << p2[2] << endl;
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;
vv[path.dir.at(t)][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0; vv[path.dir.at(t)][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
vi[path.dir.at(t)][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0; vi[path.dir.at(t)][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
vv[0][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
vi[0][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
vv[1][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
vi[1][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
vv[2][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
vi[2][path.posPath[0].at(t)][path.posPath[1].at(t)][path.posPath[2].at(t)] = 0;
} }
// cerr << "found path size: " << path.dir.size() << endl; // cerr << "found path size: " << path.dir.size() << endl;
} }
} }
} }
} }
return true;
} }

4
FDTD/operator.h
View File

@ -38,9 +38,9 @@ public:
virtual int CalcECOperator(); virtual int CalcECOperator();
//! Calculate an excitation with center of f0 and the half bandwidth fc //! Calculate an excitation with center of f0 and the half bandwidth fc
virtual void CalcGaussianPulsExcitation(double f0, double fc); virtual bool CalcGaussianPulsExcitation(double f0, double fc);
//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps //! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps
virtual void CalcSinusExcitation(double f0, int nTS); virtual bool CalcSinusExcitation(double f0, int nTS);
virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
virtual void ApplyMagneticBC(bool* dirs); virtual void ApplyMagneticBC(bool* dirs);

27
FDTD/processcurrent.cpp
View File

@ -45,6 +45,7 @@ void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop)
int ProcessCurrent::Process() int ProcessCurrent::Process()
{ {
FDTD_FLOAT**** curr = Eng->GetCurrents();
if (Enabled==false) return -1; if (Enabled==false) return -1;
if (CheckTimestep()==false) return GetNextInterval(); if (CheckTimestep()==false) return GetNextInterval();
FDTD_FLOAT current=0; FDTD_FLOAT current=0;
@ -65,26 +66,26 @@ int ProcessCurrent::Process()
} }
//x-current //x-current
for (int i=start[0];i<stop[0];++i) for (unsigned int i=start[0];i<stop[0];++i)
current+=Eng->curr[0][i][start[1]][start[2]]; current+=curr[0][i][start[1]][start[2]];
//y-current //y-current
for (int i=start[1];i<stop[1];++i) for (unsigned int i=start[1];i<stop[1];++i)
current+=Eng->curr[1][stop[0]][i][start[2]]; current+=curr[1][stop[0]][i][start[2]];
//z-current //z-current
for (int i=start[2];i<stop[2];++i) for (unsigned int i=start[2];i<stop[2];++i)
current+=Eng->curr[2][stop[0]][stop[1]][i]; current+=curr[2][stop[0]][stop[1]][i];
//x-current //x-current
for (int i=start[0];i<stop[0];++i) for (unsigned int i=start[0];i<stop[0];++i)
current-=Eng->curr[0][i][stop[1]][stop[2]]; current-=curr[0][i][stop[1]][stop[2]];
//y-current //y-current
for (int i=start[1];i<stop[1];++i) for (unsigned int i=start[1];i<stop[1];++i)
current-=Eng->curr[1][start[0]][i][stop[2]]; current-=curr[1][start[0]][i][stop[2]];
//z-current //z-current
for (int i=start[2];i<stop[2];++i) for (unsigned int i=start[2];i<stop[2];++i)
current-=Eng->curr[2][start[0]][start[1]][i]; current-=curr[2][start[0]][start[1]][i];
// cerr << "ts: " << Eng->numTS << " i: " << current << endl; // cerr << "ts: " << Eng->numTS << " i: " << current << endl;
v_current.push_back(current); v_current.push_back(current);
file << (double)Eng->numTS*Op->GetTimestep() << "\t" << current << endl; file << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << current << endl;
return GetNextInterval(); return GetNextInterval();
} }

22
FDTD/processfields.cpp
View File

@ -101,6 +101,8 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
double ProcessFields::CalcTotalEnergy() double ProcessFields::CalcTotalEnergy()
{ {
FDTD_FLOAT**** volt = Eng->GetVoltages();
FDTD_FLOAT**** curr = Eng->GetCurrents();
double energy=0; double energy=0;
if (Eng==NULL) return 0.0; if (Eng==NULL) return 0.0;
unsigned int pos[3]; unsigned int pos[3];
@ -110,12 +112,12 @@ double ProcessFields::CalcTotalEnergy()
{ {
for (pos[2]=0;pos[2]<Op->numLines[2];++pos[2]) for (pos[2]=0;pos[2]<Op->numLines[2];++pos[2])
{ {
energy+=fabs(Eng->volt[0][pos[0]][pos[1]][pos[2]] * Eng->curr[1][pos[0]][pos[1]][pos[2]]); energy+=fabs(volt[0][pos[0]][pos[1]][pos[2]] * curr[1][pos[0]][pos[1]][pos[2]]);
energy+=fabs(Eng->volt[0][pos[0]][pos[1]][pos[2]] * Eng->curr[2][pos[0]][pos[1]][pos[2]]); energy+=fabs(volt[0][pos[0]][pos[1]][pos[2]] * curr[2][pos[0]][pos[1]][pos[2]]);
energy+=fabs(Eng->volt[1][pos[0]][pos[1]][pos[2]] * Eng->curr[0][pos[0]][pos[1]][pos[2]]); energy+=fabs(volt[1][pos[0]][pos[1]][pos[2]] * curr[0][pos[0]][pos[1]][pos[2]]);
energy+=fabs(Eng->volt[1][pos[0]][pos[1]][pos[2]] * Eng->curr[2][pos[0]][pos[1]][pos[2]]); energy+=fabs(volt[1][pos[0]][pos[1]][pos[2]] * curr[2][pos[0]][pos[1]][pos[2]]);
energy+=fabs(Eng->volt[2][pos[0]][pos[1]][pos[2]] * Eng->curr[0][pos[0]][pos[1]][pos[2]]); energy+=fabs(volt[2][pos[0]][pos[1]][pos[2]] * curr[0][pos[0]][pos[1]][pos[2]]);
energy+=fabs(Eng->volt[2][pos[0]][pos[1]][pos[2]] * Eng->curr[1][pos[0]][pos[1]][pos[2]]); energy+=fabs(volt[2][pos[0]][pos[1]][pos[2]] * curr[1][pos[0]][pos[1]][pos[2]]);
} }
} }
} }
@ -183,16 +185,18 @@ bool ProcessFields::DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT*
{ {
WriteVTKHeader(file, discLines, numLines); WriteVTKHeader(file, discLines, numLines);
WriteVTKVectorArray(file, name, array, numLines); WriteVTKVectorArray(file, name, array, numLines);
return true;
} }
bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT**** array[], unsigned int numFields, double** discLines, unsigned int* numLines) bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT**** array[], unsigned int numFields, double** discLines, unsigned int* numLines)
{ {
WriteVTKHeader(file, discLines, numLines); WriteVTKHeader(file, discLines, numLines);
for (int n=0;n<numFields;++n) for (unsigned int n=0;n<numFields;++n)
{ {
WriteVTKVectorArray(file, names[n], array[n], numLines); WriteVTKVectorArray(file, names[n], array[n], numLines);
file << endl; file << endl;
} }
return true;
} }
void ProcessFields::WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT*** array, unsigned int* numLines) void ProcessFields::WriteVTKScalarArray(ofstream &file, string name, FDTD_FLOAT*** array, unsigned int* numLines)
@ -220,16 +224,18 @@ bool ProcessFields::DumpScalarArray2VTK(ofstream &file, string name, FDTD_FLOAT*
{ {
WriteVTKHeader(file, discLines, numLines); WriteVTKHeader(file, discLines, numLines);
WriteVTKScalarArray(file, name, array, numLines); WriteVTKScalarArray(file, name, array, numLines);
return true;
} }
bool ProcessFields::DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT*** array[], unsigned int numFields, double** discLines, unsigned int* numLines) bool ProcessFields::DumpMultiScalarArray2VTK(ofstream &file, string names[], FDTD_FLOAT*** array[], unsigned int numFields, double** discLines, unsigned int* numLines)
{ {
WriteVTKHeader(file, discLines, numLines); WriteVTKHeader(file, discLines, numLines);
for (int n=0;n<numFields;++n) for (unsigned int n=0;n<numFields;++n)
{ {
WriteVTKScalarArray(file, names[n], array[n], numLines); WriteVTKScalarArray(file, names[n], array[n], numLines);
file << endl; file << endl;
} }
return true;
} }

32
FDTD/processfields_td.cpp
View File

@ -31,6 +31,9 @@ ProcessFieldsTD::~ProcessFieldsTD()
void ProcessFieldsTD::DumpCellInterpol(ofstream &file) void ProcessFieldsTD::DumpCellInterpol(ofstream &file)
{ {
FDTD_FLOAT**** volt = Eng->GetVoltages();
FDTD_FLOAT**** curr = Eng->GetCurrents();
if (DumpType==0) if (DumpType==0)
{ {
//create array //create array
@ -50,15 +53,15 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file)
OpPos[2]=start[2]+pos[2]; OpPos[2]=start[2]+pos[2];
//in x //in x
delta = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; delta = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]];
E_T[0][pos[0]][pos[1]][pos[2]] = Eng->volt[0][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]] + Eng->volt[0][OpPos[0]][OpPos[1]][OpPos[2]+1] + Eng->volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]+1]; E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]+1];
E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta);
//in y //in y
delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]];
E_T[1][pos[0]][pos[1]][pos[2]] = Eng->volt[1][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]] + Eng->volt[1][OpPos[0]][OpPos[1]][OpPos[2]+1] + Eng->volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]+1]; E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]+1];
E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta);
//in z //in z
delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]];
E_T[2][pos[0]][pos[1]][pos[2]] = Eng->volt[2][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->volt[2][OpPos[0]][OpPos[1]+1][OpPos[2]] + Eng->volt[2][OpPos[0]+1][OpPos[1]][OpPos[2]] + Eng->volt[2][OpPos[0]+1][OpPos[1]+1][OpPos[2]]; E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]+1][OpPos[2]];
E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta);
} }
} }
@ -88,15 +91,15 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file)
//in x //in x
if (OpPos[0]==0) delta = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; if (OpPos[0]==0) delta = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]];
else delta = 0.5* (Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]-1]); else delta = 0.5* (Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]-1]);
H_T[0][pos[0]][pos[1]][pos[2]] = Eng->curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[0][OpPos[0]+1][OpPos[1]][OpPos[2]]; H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + curr[0][OpPos[0]+1][OpPos[1]][OpPos[2]];
H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta);
//in y //in y
delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]];
H_T[1][pos[0]][pos[1]][pos[2]] = Eng->curr[1][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[1][OpPos[0]][OpPos[1]+1][OpPos[2]]; H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]] + curr[1][OpPos[0]][OpPos[1]+1][OpPos[2]];
H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta);
//in z //in z
delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]];
H_T[2][pos[0]][pos[1]][pos[2]] = Eng->curr[2][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[2][OpPos[0]][OpPos[1]][OpPos[2]+1]; H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]] + curr[2][OpPos[0]][OpPos[1]][OpPos[2]+1];
H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta);
} }
} }
@ -109,6 +112,9 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file)
void ProcessFieldsTD::DumpNoInterpol(ofstream &file) void ProcessFieldsTD::DumpNoInterpol(ofstream &file)
{ {
FDTD_FLOAT**** volt = Eng->GetVoltages();
FDTD_FLOAT**** curr = Eng->GetCurrents();
unsigned int pos[3]; unsigned int pos[3];
double delta[3]; double delta[3];
if (DumpType==0) if (DumpType==0)
@ -124,9 +130,9 @@ void ProcessFieldsTD::DumpNoInterpol(ofstream &file)
for (pos[2]=0;pos[2]<numLines[2];++pos[2]) for (pos[2]=0;pos[2]<numLines[2];++pos[2])
{ {
delta[2]=fabs(Op->MainOp->GetIndexDelta(2,pos[2]+start[2])); delta[2]=fabs(Op->MainOp->GetIndexDelta(2,pos[2]+start[2]));
E_T[0][pos[0]][pos[1]][pos[2]] = Eng->volt[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta; E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta;
E_T[1][pos[0]][pos[1]][pos[2]] = Eng->volt[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta; E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta;
E_T[2][pos[0]][pos[1]][pos[2]] = Eng->volt[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta; E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta;
} }
} }
} }
@ -149,9 +155,9 @@ void ProcessFieldsTD::DumpNoInterpol(ofstream &file)
{ {
delta[2]=fabs(Op->MainOp->GetIndexWidth(2,pos[2]+start[2])); delta[2]=fabs(Op->MainOp->GetIndexWidth(2,pos[2]+start[2]));
//in x //in x
H_T[0][pos[0]][pos[1]][pos[2]] = Eng->curr[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta; H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta;
H_T[1][pos[0]][pos[1]][pos[2]] = Eng->curr[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta; H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta;
H_T[2][pos[0]][pos[1]][pos[2]] = Eng->curr[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta; H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta;
} }
} }
} }
@ -167,7 +173,7 @@ int ProcessFieldsTD::Process()
if (filePattern.empty()) return -1; if (filePattern.empty()) return -1;
if (CheckTimestep()==false) return GetNextInterval(); if (CheckTimestep()==false) return GetNextInterval();
stringstream ss; stringstream ss;
ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->numTS; ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps();
string filename = filePattern + ss.str() + ".vtk"; string filename = filePattern + ss.str() + ".vtk";
ofstream file(filename.c_str()); ofstream file(filename.c_str());

4
FDTD/processing.cpp
View File

@ -51,9 +51,9 @@ double Processing::CalcLineIntegral(unsigned int* start, unsigned int* stop, int
double result=0; double result=0;
FDTD_FLOAT**** array; FDTD_FLOAT**** array;
if (field==0) if (field==0)
array=Eng->volt; array=Eng->GetVoltages();
else if (field==1) else if (field==1)
array=Eng->curr; array=Eng->GetCurrents();
else return 0.0; else return 0.0;
for (int n=0;n<3;++n) for (int n=0;n<3;++n)

2
FDTD/processvoltage.cpp
View File

@ -45,6 +45,6 @@ int ProcessVoltage::Process()
FDTD_FLOAT voltage=CalcLineIntegral(start,stop,0); FDTD_FLOAT voltage=CalcLineIntegral(start,stop,0);
// cerr << voltage << endl; // cerr << voltage << endl;
voltages.push_back(voltage); voltages.push_back(voltage);
file << (double)Eng->numTS*Op->GetTimestep() << "\t" << voltage << endl; file << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl;
return GetNextInterval(); return GetNextInterval();
} }

44
openems.cpp
View File

@ -138,7 +138,12 @@ int openEMS::SetupFDTD(const char* file)
cerr << "Can't read openEMS FDTD Settings... " << endl; cerr << "Can't read openEMS FDTD Settings... " << endl;
exit(-1); exit(-1);
} }
FDTD_Opts->QueryIntAttribute("NumberOfTimesteps",&NrTS); int help=0;
FDTD_Opts->QueryIntAttribute("NumberOfTimesteps",&help);
if (help<0)
NrTS=0;
else
NrTS = help;
FDTD_Opts->QueryDoubleAttribute("endCriteria",&endCrit); FDTD_Opts->QueryDoubleAttribute("endCriteria",&endCrit);
if (endCrit==0) if (endCrit==0)
endCrit=1e-6; endCrit=1e-6;
@ -192,26 +197,39 @@ int openEMS::SetupFDTD(const char* file)
FDTD_Op = new Operator(); FDTD_Op = new Operator();
if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1); if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1);
if (DebugMat)
{
FDTD_Op->DumpMaterial2File("material_dump.vtk");
}
FDTD_Op->CalcECOperator(); FDTD_Op->CalcECOperator();
if (DebugOp)
{
FDTD_Op->DumpOperator2File("operator_dump.vtk");
}
if (Excit_Type==0) if (Excit_Type==0)
FDTD_Op->CalcGaussianPulsExcitation(f0,fc); {
if (!FDTD_Op->CalcGaussianPulsExcitation(f0,fc))
{
cerr << "openEMS: excitation setup failed!!" << endl;
exit(2);
}
}
else if (Excit_Type==1) else if (Excit_Type==1)
FDTD_Op->CalcSinusExcitation(f0,NrTS); {
if (!FDTD_Op->CalcSinusExcitation(f0,NrTS))
{
cerr << "openEMS: excitation setup failed!!" << endl;
exit(2);
}
}
else else
{ {
cerr << "openEMS: Excitation type is unknown" << endl; cerr << "openEMS: Excitation type is unknown" << endl;
exit(-1); exit(-1);
} }
if (DebugMat)
{
FDTD_Op->DumpMaterial2File("material_dump.vtk");
}
if (DebugOp)
{
FDTD_Op->DumpOperator2File("operator_dump.vtk");
}
time_t OpDoneTime=time(NULL); time_t OpDoneTime=time(NULL);
FDTD_Op->ShowSize(); FDTD_Op->ShowSize();
@ -330,14 +348,14 @@ void openEMS::RunFDTD()
//*************** simulate ************// //*************** simulate ************//
int step=PA->Process(); int step=PA->Process();
if ((step<0) || (step>NrTS)) step=NrTS; if ((step<0) || (step>(int)NrTS)) step=NrTS;
while ((FDTD_Eng->GetNumberOfTimesteps()<NrTS) && (change>endCrit)) while ((FDTD_Eng->GetNumberOfTimesteps()<NrTS) && (change>endCrit))
{ {
FDTD_Eng->IterateTS(step); FDTD_Eng->IterateTS(step);
step=PA->Process(); step=PA->Process();
// cout << " do " << step << " steps; current: " << eng.GetNumberOfTimesteps() << endl; // cout << " do " << step << " steps; current: " << eng.GetNumberOfTimesteps() << endl;
currTS = FDTD_Eng->GetNumberOfTimesteps(); currTS = FDTD_Eng->GetNumberOfTimesteps();
if ((step<0) || (step>NrTS - currTS)) step=NrTS - currTS; if ((step<0) || (step>(int)(NrTS - currTS))) step=NrTS - currTS;
gettimeofday(&currTime,NULL); gettimeofday(&currTime,NULL);

2
openems.h
View File

@ -44,7 +44,7 @@ public:
protected: protected:
//! Number of Timesteps //! Number of Timesteps
int NrTS; unsigned int NrTS;
bool Enable_Dumps; bool Enable_Dumps;
bool DebugMat; bool DebugMat;
bool DebugOp; bool DebugOp;