diff --git a/FDTD/engine.h b/FDTD/engine.h index 401fa13..e1660b7 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -22,11 +22,6 @@ class Engine { - friend class Processing; - friend class ProcessVoltage; - friend class ProcessCurrent; - friend class ProcessFields; - friend class ProcessFieldsTD; public: Engine(Operator* op); virtual ~Engine(); @@ -39,6 +34,9 @@ public: unsigned int GetNumberOfTimesteps() {return numTS;}; + virtual FDTD_FLOAT**** GetVoltages() {return volt;}; + virtual FDTD_FLOAT**** GetCurrents() {return curr;}; + protected: Operator* Op; diff --git a/FDTD/processcurrent.cpp b/FDTD/processcurrent.cpp index 66c23a0..a63a997 100644 --- a/FDTD/processcurrent.cpp +++ b/FDTD/processcurrent.cpp @@ -45,6 +45,7 @@ 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; @@ -66,25 +67,25 @@ int ProcessCurrent::Process() //x-current for (unsigned int i=start[0];icurr[0][i][start[1]][start[2]]; + current+=curr[0][i][start[1]][start[2]]; //y-current for (unsigned int i=start[1];icurr[1][stop[0]][i][start[2]]; + current+=curr[1][stop[0]][i][start[2]]; //z-current for (unsigned int i=start[2];icurr[2][stop[0]][stop[1]][i]; + current+=curr[2][stop[0]][stop[1]][i]; //x-current for (unsigned int i=start[0];icurr[0][i][stop[1]][stop[2]]; + current-=curr[0][i][stop[1]][stop[2]]; //y-current for (unsigned int i=start[1];icurr[1][start[0]][i][stop[2]]; + current-=curr[1][start[0]][i][stop[2]]; //z-current for (unsigned int i=start[2];icurr[2][start[0]][start[1]][i]; + current-=curr[2][start[0]][start[1]][i]; // cerr << "ts: " << Eng->numTS << " i: " << current << endl; 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(); } diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp index 3d538c6..4c46261 100644 --- a/FDTD/processfields.cpp +++ b/FDTD/processfields.cpp @@ -101,6 +101,8 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) double ProcessFields::CalcTotalEnergy() { + FDTD_FLOAT**** volt = Eng->GetVoltages(); + FDTD_FLOAT**** curr = Eng->GetCurrents(); double energy=0; if (Eng==NULL) return 0.0; unsigned int pos[3]; @@ -110,12 +112,12 @@ double ProcessFields::CalcTotalEnergy() { for (pos[2]=0;pos[2]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(Eng->volt[0][pos[0]][pos[1]][pos[2]] * Eng->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(Eng->volt[1][pos[0]][pos[1]][pos[2]] * Eng->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(Eng->volt[2][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(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]]); } } } diff --git a/FDTD/processfields_td.cpp b/FDTD/processfields_td.cpp index e95d9b5..bc8a93c 100644 --- a/FDTD/processfields_td.cpp +++ b/FDTD/processfields_td.cpp @@ -31,6 +31,9 @@ ProcessFieldsTD::~ProcessFieldsTD() void ProcessFieldsTD::DumpCellInterpol(ofstream &file) { + FDTD_FLOAT**** volt = Eng->GetVoltages(); + FDTD_FLOAT**** curr = Eng->GetCurrents(); + if (DumpType==0) { //create array @@ -50,15 +53,15 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file) OpPos[2]=start[2]+pos[2]; //in x 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); //in y 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); //in z 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); } } @@ -88,15 +91,15 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file) //in x 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]); - 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); //in y 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); //in z 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); } } @@ -109,6 +112,9 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file) void ProcessFieldsTD::DumpNoInterpol(ofstream &file) { + FDTD_FLOAT**** volt = Eng->GetVoltages(); + FDTD_FLOAT**** curr = Eng->GetCurrents(); + unsigned int pos[3]; double delta[3]; if (DumpType==0) @@ -124,9 +130,9 @@ void ProcessFieldsTD::DumpNoInterpol(ofstream &file) for (pos[2]=0;pos[2]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[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[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[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]] = 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]] = 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])); //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[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[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[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]] = 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]] = 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 (CheckTimestep()==false) return GetNextInterval(); 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"; ofstream file(filename.c_str()); diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp index 5be1f68..a54dcca 100644 --- a/FDTD/processing.cpp +++ b/FDTD/processing.cpp @@ -51,9 +51,9 @@ double Processing::CalcLineIntegral(unsigned int* start, unsigned int* stop, int double result=0; FDTD_FLOAT**** array; if (field==0) - array=Eng->volt; + array=Eng->GetVoltages(); else if (field==1) - array=Eng->curr; + array=Eng->GetCurrents(); else return 0.0; for (int n=0;n<3;++n) diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp index 25fe31a..64d41b0 100644 --- a/FDTD/processvoltage.cpp +++ b/FDTD/processvoltage.cpp @@ -45,6 +45,6 @@ int ProcessVoltage::Process() FDTD_FLOAT voltage=CalcLineIntegral(start,stop,0); // cerr << voltage << endl; voltages.push_back(voltage); - file << (double)Eng->numTS*Op->GetTimestep() << "\t" << voltage << endl; + file << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl; return GetNextInterval(); }