diff --git a/Analyse/PlotVoltage.m b/Analyse/PlotVoltage.m index 34e037b..d2ad72f 100644 --- a/Analyse/PlotVoltage.m +++ b/Analyse/PlotVoltage.m @@ -11,6 +11,7 @@ L=numel(t); subplot(2,1,1); plot(t,u); +grid on; dt=t(2)-t(1); @@ -18,4 +19,5 @@ f = (1:L)/L/dt; fu = fft(u)/L; subplot(2,1,2); plot(f(1:L/2),abs(fu(1:L/2))); +grid on; diff --git a/FDTD/cartoperator.cpp b/FDTD/cartoperator.cpp index bcea047..8502401 100644 --- a/FDTD/cartoperator.cpp +++ b/FDTD/cartoperator.cpp @@ -432,14 +432,22 @@ bool CartOperator::CalcEFieldExcitation() if (prop) { CSPropElectrode* elec = prop->ToElectrode(); - if (elec->GetExcitType()==0) + if ((elec->GetExcitType()==0) || (elec->GetExcitType()==1)) //soft or hard E-Field excite! { vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); for (int n=0;n<3;++n) { vIndex[n].push_back(pos[n]); double delta=MainOp->GetIndexDelta(n,pos[n])*gridDelta; - vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta); + if (elec->GetActiveDir(n)) + vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta); + else + vExcit[n].push_back(0); + if ((elec->GetExcitType()==1) && (elec->GetActiveDir(n))) //hard excite + { + vv[n][pos[0]][pos[1]][pos[2]] = 0; + vi[n][pos[0]][pos[1]][pos[2]] = 0; + } } } } diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index e0fd0f4..4090394 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -62,7 +62,8 @@ bool Engine::IterateTS(unsigned int iterTS) for (unsigned int n=0;nE_Ex_Count;++n) { exc_pos = (int)numTS - (int)Op->E_Ex_delay[n]; - exc_pos*= (exc_pos>0 && exc_posExciteLength); + exc_pos*= (exc_pos>0 && exc_pos<(int)Op->ExciteLength); +// if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl; volt[0][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[0][n]*Op->ExciteSignal[exc_pos]; volt[1][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[1][n]*Op->ExciteSignal[exc_pos]; volt[2][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[2][n]*Op->ExciteSignal[exc_pos]; diff --git a/FDTD/engine.h b/FDTD/engine.h index 5fd4dcb..44c92eb 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -6,6 +6,7 @@ class Engine { friend class ProcessVoltage; + friend class ProcessFieldsTD; public: Engine(Operator* op); virtual ~Engine(); diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 30d15e0..96b0f9a 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -54,27 +54,33 @@ unsigned int Operator::GetNyquistNum(double fmax) return floor(T0/2/dT); } -bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord) +bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower) { bool ok=true; for (int n=0;n<3;++n) { - if (dcoord[n]discLines[n][numLines[n]-1]) {ok=false;uicoord[n]=numLines[n]-1;}; - for (unsigned int i=0;idiscLines[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 + for (unsigned int i=0;iSnapToMesh(dstart,start,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl; + if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl; + + for (int n=0;n<3;++n) + { +// cerr << "start " << start[n] << "stop " << stop[n]; + if (start[n]>stop[n]) + { + unsigned int help = start[n]; + start[n]=stop[n]; + stop[n]=help; + } + ++stop[n]; + numDLines[n]=stop[n]-start[n]; +// cerr << " number of lines " << numDLines[n] << endl; + delete[] discDLines[n]; + discDLines[n] = new double[numDLines[n]]; + for (unsigned int i=0;idiscLines[n][start[n]+i+1] + Op->discLines[n][start[n]+i]); +// cerr << n << " : " << discDLines[n][i] << endl; + } + } +} + +bool ProcessFields::DumpFieldArray2VTK(ofstream &file, string name, FDTD_FLOAT**** array, double** discLines, unsigned int* numLines) +{ + file << "# vtk DataFile Version 2.0" << endl; + file << "Rectilinear Grid openEMS_ProcessFields" << endl; + file << "ASCII" << endl; + file << "DATASET RECTILINEAR_GRID " << endl; + file << "DIMENSIONS " << numLines[0] << " " << numLines[1] << " " << numLines[2] << endl; + file << "X_COORDINATES " << numLines[0] << " float" << endl; + for (unsigned int i=0;i "tmp/efield_000045.vtk" for timestep 45 or "tmp/efield_2.40000e9.vtk" for 2.4GHz E-field dump. + void SetFilePattern(string fp) {filePattern=fp;} + + //! This methode will dump all fields in the center of a main cell (dual-node) using 4 E-field and 2 H-fields per direction. (default) + void SetDumpMode2Cell() {DumpMode=0;} + + //! Set dump type: 0 for E-fields, 1 for H-fields, 2 for D-fields, 3 for B-fields, 4 for J-fields, etc... + void SetDumpType(int type) {DumpType=type;} + +// virtual void Process(); +protected: + ProcessFields(Operator* op, Engine* eng); + + bool DumpFieldArray2VTK(ofstream &file, string name, FDTD_FLOAT**** array, double** discLines, unsigned int* numLines); + + int DumpMode; + int DumpType; + string filePattern; + + unsigned int numDLines[3]; + double* discDLines[3]; +}; + +#endif // PROCESSFIELDS_H diff --git a/FDTD/processfields_td.cpp b/FDTD/processfields_td.cpp new file mode 100644 index 0000000..5f5dec7 --- /dev/null +++ b/FDTD/processfields_td.cpp @@ -0,0 +1,61 @@ +#include "processfields_td.h" +#include +#include +#include + +ProcessFieldsTD::ProcessFieldsTD(Operator* op, Engine* eng) : ProcessFields(op, eng) +{ + pad_length = 8; +} + +ProcessFieldsTD::~ProcessFieldsTD() +{ +} + +void ProcessFieldsTD::Process() +{ + stringstream ss; + ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->numTS; + + string filename = filePattern + ss.str() + ".vtk"; + ofstream file(filename.c_str()); + if (file.is_open()==false) { cerr << "ProcessFieldsTD::Process: can't open file '" << filename << "' for writing... abort! " << endl; return;}; + + if (DumpType==0) + { + //create array + FDTD_FLOAT**** E_T = Create_N_3DArray(numDLines); + unsigned int pos[3] = {start[0],start[1],start[2]}; + unsigned int OpPos[3]; + double delta; +// cerr << "processing e-fields... " << endl; + for (pos[0]=0;pos[0]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]] /= delta; + //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]] /= delta; + //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]] /= delta; + } + } + } + DumpFieldArray2VTK(file,string("E-Field"),E_T,discDLines,numDLines); + Delete_N_3DArray(E_T,numDLines); + E_T = NULL; + } + file.close(); +} diff --git a/FDTD/processfields_td.h b/FDTD/processfields_td.h new file mode 100644 index 0000000..a0db4e3 --- /dev/null +++ b/FDTD/processfields_td.h @@ -0,0 +1,21 @@ +#ifndef PROCESSFIELDS_TD_H +#define PROCESSFIELDS_TD_H + +#include "processfields.h" + +class ProcessFieldsTD : public ProcessFields +{ +public: + ProcessFieldsTD(Operator* op, Engine* eng); + virtual ~ProcessFieldsTD(); + + virtual void Process(); + + //! Set the length of the filename timestep pad filled with zeros (default is 8) + void SetPadLength(int val) {pad_length=val;}; + +protected: + int pad_length; +}; + +#endif // PROCESSFIELDS_TD_H diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp index 34c7e68..4aa4c90 100644 --- a/FDTD/processing.cpp +++ b/FDTD/processing.cpp @@ -8,17 +8,12 @@ Processing::Processing(Operator* op, Engine* eng) Processing::~Processing() { - file.close(); } -void Processing::OpenFile(string outfile) + +void Processing::DefineStartStopCoord(double* dstart, double* dstop) { - if (file.is_open()) file.close(); - - file.open(outfile.c_str()); - if (file.is_open()==false) - { - cerr << "Can't open file: " << outfile << endl; - return; - } + if (Op->SnapToMesh(dstart,start)==false) cerr << "Processing::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl; + if (Op->SnapToMesh(dstop,stop)==false) cerr << "Processing::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl; } + diff --git a/FDTD/processing.h b/FDTD/processing.h index 5eae7e7..08b7c72 100644 --- a/FDTD/processing.h +++ b/FDTD/processing.h @@ -9,17 +9,19 @@ class Processing { public: - Processing(Operator* op, Engine* eng); virtual ~Processing(); - virtual void OpenFile(string outfile); + virtual void DefineStartStopCoord(double* dstart, double* dstop); virtual void Process() {}; protected: + Processing(Operator* op, Engine* eng); Operator* Op; Engine* Eng; - ofstream file; + + unsigned int start[3]; + unsigned int stop[3]; }; #endif // PROCESSING_H diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp index 95ad474..f9444fb 100644 --- a/FDTD/processvoltage.cpp +++ b/FDTD/processvoltage.cpp @@ -6,12 +6,19 @@ ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng) ProcessVoltage::~ProcessVoltage() { + file.close(); } -void ProcessVoltage::DefineStartStopCoord(double* dstart, double* dstop) +void ProcessVoltage::OpenFile(string outfile) { - if (Op->SnapToMesh(dstart,start)==false) cerr << "ProcessVoltage::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl; - if (Op->SnapToMesh(dstop,stop)==false) cerr << "ProcessVoltage::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl; + if (file.is_open()) file.close(); + + file.open(outfile.c_str()); + if (file.is_open()==false) + { + cerr << "Can't open file: " << outfile << endl; + return; + } } void ProcessVoltage::Process() diff --git a/FDTD/processvoltage.h b/FDTD/processvoltage.h index dccedf4..604e827 100644 --- a/FDTD/processvoltage.h +++ b/FDTD/processvoltage.h @@ -3,19 +3,19 @@ #include "processing.h" +//! Process voltage along a line from start to stop coordinates. ATM integration along the axis e.g.: in x, then y then z direction (Future: diagonal integration) class ProcessVoltage : public Processing { public: ProcessVoltage(Operator* op, Engine* eng); virtual ~ProcessVoltage(); - void DefineStartStopCoord(double* dstart, double* dstop); + virtual void OpenFile(string outfile); virtual void Process(); protected: - unsigned int start[3]; - unsigned int stop[3]; + ofstream file; vector voltages; }; diff --git a/main.cpp b/main.cpp index 1ea811a..9221cd9 100644 --- a/main.cpp +++ b/main.cpp @@ -6,6 +6,7 @@ #include "FDTD/cartoperator.h" #include "FDTD/engine.h" #include "FDTD/processvoltage.h" +#include "FDTD/processfields_td.h" #include "ContinuousStructure.h" @@ -35,8 +36,6 @@ int main(int argc, char *argv[]) cop.ShowSize(); -// cop.DumpOperator2File("tmp/Operator"); - cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl; unsigned int NrIter = cop.GetNyquistNum(fmax)/3; @@ -50,18 +49,26 @@ int main(int argc, char *argv[]) //*************** setup processing ************// ProcessVoltage PV(&cop,&eng); PV.OpenFile("tmp/u1"); - double start[]={-500,-150,0}; - double stop[]={-500,150,0}; + double start[]={-0,-75,0}; + double stop[]={-0,75,0}; PV.DefineStartStopCoord(start,stop); unsigned int maxIter = 5000; + ProcessFieldsTD PETD(&cop,&eng); + start[0]=-1000;start[1]=0;start[2]=-1000; + stop[0]=1000;stop[1]=0;stop[2]=1000; + PETD.SetFilePattern("tmp/Et_"); + PETD.DefineStartStopCoord(start,stop); + PV.Process(); -// NrIter=200; + PETD.Process(); + //*************** simulate ************// for (unsigned int i=0;iSetExcitation(1,1); - elec->SetExcitType(0); + elec->SetExcitType(1); + elec->SetActiveDir(0,0);//disable x + elec->SetActiveDir(0,2);//disable z // elec->SetDelay(2.0e-9); CSX.AddProperty(elec); diff --git a/openEMS.pro b/openEMS.pro index e4493f5..d4ac0d6 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -22,7 +22,9 @@ SOURCES += main.cpp \ FDTD/operator.cpp \ tools/array_ops.cpp \ FDTD/processvoltage.cpp \ - FDTD/processing.cpp + FDTD/processing.cpp \ + FDTD/processfields.cpp \ + FDTD/processfields_td.cpp HEADERS += FDTD/cartoperator.h \ tools/ErrorMsg.h \ tools/AdrOp.h \ @@ -31,4 +33,6 @@ HEADERS += FDTD/cartoperator.h \ FDTD/operator.h \ tools/array_ops.h \ FDTD/processvoltage.h \ - FDTD/processing.h + FDTD/processing.h \ + FDTD/processfields.h \ + FDTD/processfields_td.h diff --git a/tools/array_ops.cpp b/tools/array_ops.cpp index 945a8c2..742e1ef 100644 --- a/tools/array_ops.cpp +++ b/tools/array_ops.cpp @@ -3,7 +3,7 @@ FDTD_FLOAT*** Create3DArray(unsigned int* numLines) { - FDTD_FLOAT*** array; + FDTD_FLOAT*** array=NULL; unsigned int pos[3]; array = new FDTD_FLOAT**[numLines[0]]; for (pos[0]=0;pos[0]