From e9cdadb01a2322ca2da3562889b7e26bcac3fd1c Mon Sep 17 00:00:00 2001 From: Sebastian Held Date: Wed, 21 Apr 2010 14:28:16 +0200 Subject: [PATCH] processing functions are now independent of the underlying engine data format --- FDTD/engine.h | 4 ++-- FDTD/processcurrent.cpp | 13 +++++------ FDTD/processfields.cpp | 37 +++++++++++++++--------------- FDTD/processfields.h | 18 +++++++-------- FDTD/processfields_td.cpp | 30 ++++++++++-------------- FDTD/processing.cpp | 48 +++++++++++++++++++++++++++------------ FDTD/processing.h | 4 +++- 7 files changed, 84 insertions(+), 70 deletions(-) diff --git a/FDTD/engine.h b/FDTD/engine.h index a30c682..0e83249 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -35,8 +35,8 @@ public: virtual unsigned int GetNumberOfTimesteps() {return numTS;}; - virtual FDTD_FLOAT**** GetVoltages() {return volt;}; - virtual FDTD_FLOAT**** GetCurrents() {return curr;}; + inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return volt[n][x][y][z]; } + inline virtual FDTD_FLOAT GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return curr[n][x][y][z]; } protected: Engine(const Operator* op); diff --git a/FDTD/processcurrent.cpp b/FDTD/processcurrent.cpp index eec845c..8176c9d 100644 --- a/FDTD/processcurrent.cpp +++ b/FDTD/processcurrent.cpp @@ -34,7 +34,6 @@ void ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop) int ProcessCurrent::Process() { - FDTD_FLOAT**** curr = Eng->GetCurrents(); if (Enabled==false) return -1; if (CheckTimestep()==false) return GetNextInterval(); FDTD_FLOAT current=0; @@ -56,27 +55,27 @@ int ProcessCurrent::Process() //x-current if (m_start_inside[1] && m_start_inside[2]) for (unsigned int i=start[0];iGetCurr(0,i,start[1],start[2]); //y-current if (m_stop_inside[0] && m_start_inside[2]) for (unsigned int i=start[1];iGetCurr(1,stop[0],i,start[2]); //z-current if (m_stop_inside[0] && m_stop_inside[1]) for (unsigned int i=start[2];iGetCurr(2,stop[0],stop[1],i); //x-current if (m_stop_inside[1] && m_stop_inside[2]) for (unsigned int i=start[0];iGetCurr(0,i,stop[1],stop[2]); //y-current if (m_start_inside[0] && m_stop_inside[2]) for (unsigned int i=start[1];iGetCurr(1,start[0],i,stop[2]); //z-current if (m_start_inside[0] && m_start_inside[1]) for (unsigned int i=start[2];iGetCurr(2,start[0],start[1],i); // cerr << "ts: " << Eng->numTS << " i: " << current << endl; v_current.push_back(current); diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp index 0ef6cb6..77b2bb6 100644 --- a/FDTD/processfields.cpp +++ b/FDTD/processfields.cpp @@ -168,13 +168,12 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) } } - -double ProcessFields::CalcTotalEnergy() +double ProcessFields::CalcTotalEnergy() const { - FDTD_FLOAT**** volt = Eng->GetVoltages(); - FDTD_FLOAT**** curr = Eng->GetCurrents(); + if (!Eng) + return 0; + double energy=0; - if (Eng==NULL) return 0.0; unsigned int pos[3]; for (pos[0]=0;pos[0]GetNumberOfLines(0);++pos[0]) { @@ -182,12 +181,12 @@ double ProcessFields::CalcTotalEnergy() { for (pos[2]=0;pos[2]GetNumberOfLines(2);++pos[2]) { - energy+=fabs(volt[0][pos[0]][pos[1]][pos[2]] * curr[1][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(volt[1][pos[0]][pos[1]][pos[2]] * curr[0][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(volt[2][pos[0]][pos[1]][pos[2]] * curr[0][pos[0]][pos[1]][pos[2]]); - energy+=fabs(volt[2][pos[0]][pos[1]][pos[2]] * curr[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(1,pos[0],pos[1],pos[2]) * Eng->GetCurr(0,pos[0],pos[1],pos[2])); + energy+=fabs(Eng->GetVolt(1,pos[0],pos[1],pos[2]) * Eng->GetCurr(2,pos[0],pos[1],pos[2])); + energy+=fabs(Eng->GetVolt(2,pos[0],pos[1],pos[2]) * Eng->GetCurr(0,pos[0],pos[1],pos[2])); + energy+=fabs(Eng->GetVolt(2,pos[0],pos[1],pos[2]) * Eng->GetCurr(1,pos[0],pos[1],pos[2])); } } } @@ -206,7 +205,7 @@ void ProcessFields::SetSubSampling(unsigned int subSampleRate, int dir) else subSample[dir]=subSampleRate; } -void ProcessFields::WriteVTKHeader(ofstream &file, double** discLines, unsigned int* numLines, unsigned int precision) +void ProcessFields::WriteVTKHeader(ofstream &file, double const* const* discLines, unsigned int const* numLines, unsigned int precision) { file << "# vtk DataFile Version 2.0" << endl; file << "Rectilinear Grid openEMS_ProcessFields" << endl; @@ -228,7 +227,7 @@ void ProcessFields::WriteVTKHeader(ofstream &file, double** discLines, unsigned file << "POINT_DATA " << numLines[0]*numLines[1]*numLines[2] << endl; } -void ProcessFields::WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT**** array, unsigned int* numLines, unsigned int precision) +void ProcessFields::WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT const* const* const* const* array, unsigned int const* numLines, unsigned int precision) { file << "VECTORS " << name << " float " << endl; @@ -251,14 +250,14 @@ void ProcessFields::WriteVTKVectorArray(ofstream &file, string name, FDTD_FLOAT* } -bool ProcessFields::DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT**** array, double** discLines, unsigned int* numLines, unsigned int precision) +bool ProcessFields::DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT const* const* const* const* array, double const* const* discLines, unsigned int const* numLines, unsigned int precision) { WriteVTKHeader(file, discLines, numLines, precision); WriteVTKVectorArray(file, name, array, numLines, precision); return true; } -bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT**** array[], unsigned int numFields, double** discLines, unsigned int* numLines, unsigned int precision) +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) { WriteVTKHeader(file, discLines, numLines, precision); for (unsigned int n=0;nGetVoltages(); - FDTD_FLOAT**** curr = Eng->GetCurrents(); - if (m_DumpType==E_FIELD_DUMP) { //create array @@ -55,21 +52,21 @@ void ProcessFieldsTD::DumpCellInterpol(string filename) delta = Op->GetMeshDelta(0,OpPos); //Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; if (delta) { - 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]] = Eng->GetVolt(0,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetVolt(0,OpPos[0],OpPos[1]+1,OpPos[2]) + Eng->GetVolt(0,OpPos[0],OpPos[1],OpPos[2]+1) + Eng->GetVolt(0,OpPos[0],OpPos[1]+1,OpPos[2]+1); E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); } //in y delta = Op->GetMeshDelta(1,OpPos); //Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; if (delta) { - 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]] = Eng->GetVolt(1,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetVolt(1,OpPos[0]+1,OpPos[1],OpPos[2]) + Eng->GetVolt(1,OpPos[0],OpPos[1],OpPos[2]+1) + Eng->GetVolt(1,OpPos[0]+1,OpPos[1],OpPos[2]+1); E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); } //in z delta = Op->GetMeshDelta(2,OpPos); //Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; if (delta) { - 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]] = Eng->GetVolt(2,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetVolt(2,OpPos[0],OpPos[1]+1,OpPos[2]) + Eng->GetVolt(2,OpPos[0]+1,OpPos[1],OpPos[2]) + Eng->GetVolt(2,OpPos[0]+1,OpPos[1]+1,OpPos[2]); E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta);//*Op->gridDelta); } } @@ -118,7 +115,7 @@ void ProcessFieldsTD::DumpCellInterpol(string filename) delta = Op->GetDiscLine(0,OpPos[0],true); if (delta) { - 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]] = Eng->GetCurr(0,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetCurr(0,OpPos[0]+1,OpPos[1],OpPos[2]); H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta); } //in y @@ -126,7 +123,7 @@ void ProcessFieldsTD::DumpCellInterpol(string filename) delta = Op->GetDiscLine(1,OpPos[1],true); if (delta) { - 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]] = Eng->GetCurr(1,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetCurr(1,OpPos[0],OpPos[1]+1,OpPos[2]); H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta); } //in z @@ -134,7 +131,7 @@ void ProcessFieldsTD::DumpCellInterpol(string filename) delta = Op->GetDiscLine(2,OpPos[2],true); if (delta) { - 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]] = Eng->GetCurr(2,OpPos[0],OpPos[1],OpPos[2]) + Eng->GetCurr(2,OpPos[0],OpPos[1],OpPos[2]+1); H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta);//*Op->gridDelta); } } @@ -162,9 +159,6 @@ void ProcessFieldsTD::DumpCellInterpol(string filename) void ProcessFieldsTD::DumpNoInterpol(string filename) { - FDTD_FLOAT**** volt = Eng->GetVoltages(); - FDTD_FLOAT**** curr = Eng->GetCurrents(); - unsigned int pos[3]; unsigned int OpPos[3]; double delta[3]; @@ -185,11 +179,11 @@ void ProcessFieldsTD::DumpNoInterpol(string filename) OpPos[2]=start[2]+pos[2]*subSample[2]; delta[2]=Op->GetMeshDelta(2,OpPos);//fabs(Op->MainOp->GetIndexDelta(2,OpPos[2])); if (delta[0]) - E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0];// /Op->gridDelta; + E_T[0][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(0,OpPos[0],OpPos[1],OpPos[2])/delta[0];// /Op->gridDelta; if (delta[1]) - E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1];// /Op->gridDelta; + E_T[1][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(1,OpPos[0],OpPos[1],OpPos[2])/delta[1];// /Op->gridDelta; if (delta[2]) - E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2];// /Op->gridDelta; + E_T[2][pos[0]][pos[1]][pos[2]] = Eng->GetVolt(2,OpPos[0],OpPos[1],OpPos[2])/delta[2];// /Op->gridDelta; } } } @@ -231,11 +225,11 @@ void ProcessFieldsTD::DumpNoInterpol(string filename) delta[2]=Op->GetMeshDelta(2,OpPos,true);//fabs(Op->MainOp->GetIndexWidth(2,OpPos[2])); //in x if (delta[0]) - H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]]/delta[0];// /Op->gridDelta; + H_T[0][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(0,OpPos[0],OpPos[1],OpPos[2])/delta[0];// /Op->gridDelta; if (delta[1]) - H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]]/delta[1];// /Op->gridDelta; + H_T[1][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(1,OpPos[0],OpPos[1],OpPos[2])/delta[1];// /Op->gridDelta; if (delta[2]) - H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]]/delta[2];// /Op->gridDelta; + H_T[2][pos[0]][pos[1]][pos[2]] = Eng->GetCurr(2,OpPos[0],OpPos[1],OpPos[2])/delta[2];// /Op->gridDelta; } } } diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp index 8467da9..d557fd6 100644 --- a/FDTD/processing.cpp +++ b/FDTD/processing.cpp @@ -91,33 +91,53 @@ void Processing::DefineStartStopCoord(double* dstart, double* dstop) if (Op->SnapToMesh(dstop,stop)==false) cerr << "Processing::DefineStartStopCoord: Warning: Snapped line outside field domain!!" << endl; } -double Processing::CalcLineIntegral(unsigned int* start, unsigned int* stop, int field) +double Processing::CalcLineIntegral(unsigned int* start, unsigned int* stop, int field) const +{ + switch (field) { + case 0: + return CalcLineIntegral_V( start, stop ); + case 1: + return CalcLineIntegral_I( start, stop ); + } + return 0; +} + +double Processing::CalcLineIntegral_I(unsigned int* start, unsigned int* stop) const { double result=0; - FDTD_FLOAT**** array; - if (field==0) - array=Eng->GetVoltages(); - else if (field==1) - array=Eng->GetCurrents(); - else return 0.0; - for (int n=0;n<3;++n) { if (start[n]GetCurr(n,pos[0],pos[1],pos[2]); } else { unsigned int pos[3]={stop[0],stop[1],stop[2]}; for (;pos[n]GetCurr(n,pos[0],pos[1],pos[2]); + } + } + return result; +} +double Processing::CalcLineIntegral_V(unsigned int* start, unsigned int* stop) const +{ + double result=0; + for (int n=0;n<3;++n) + { + if (start[n]GetVolt(n,pos[0],pos[1],pos[2]); + } + else + { + unsigned int pos[3]={stop[0],stop[1],stop[2]}; + for (;pos[n]GetVolt(n,pos[0],pos[1],pos[2]); } } return result; diff --git a/FDTD/processing.h b/FDTD/processing.h index 5224b51..95069a2 100644 --- a/FDTD/processing.h +++ b/FDTD/processing.h @@ -76,7 +76,9 @@ protected: ofstream file; string m_filename; - double CalcLineIntegral(unsigned int* start, unsigned int* stop, int field); + double CalcLineIntegral(unsigned int* start, unsigned int* stop, int field) const; + double CalcLineIntegral_I(unsigned int* start, unsigned int* stop) const; + double CalcLineIntegral_V(unsigned int* start, unsigned int* stop) const; }; class ProcessingArray