diff --git a/.gitignore b/.gitignore index c82e591..1793f37 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ release/ debug/ doc/ obj/ +tmp* Makefile* *~ *.so* diff --git a/FDTD/cartoperator.cpp b/FDTD/cartoperator.cpp index ec919ec..bcea047 100644 --- a/FDTD/cartoperator.cpp +++ b/FDTD/cartoperator.cpp @@ -1,7 +1,7 @@ #include "cartoperator.h" #include "tools/array_ops.h" -CartOperator::CartOperator() +CartOperator::CartOperator() : Operator() { Init(); } @@ -325,6 +325,7 @@ bool CartOperator::Calc_EffMatPos(int n, unsigned int* pos, double* inMat) inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta; inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta; + return true; } diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index 2b63f6f..263af8d 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -78,15 +78,15 @@ bool Engine::IterateTS(unsigned int iterTS) //do the updates here //for x curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]]; - curr[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]+1][pos[2]] - volt[2][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]+1] + volt[1][pos[0]][pos[1]][pos[2]]); + curr[0][pos[0]][pos[1]][pos[2]] -= Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]+1][pos[2]] - volt[2][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]+1] + volt[1][pos[0]][pos[1]][pos[2]]); //for y curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]]; - curr[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]+1] - volt[0][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]+1][pos[1]][pos[2]] + volt[2][pos[0]][pos[1]][pos[2]]); + curr[1][pos[0]][pos[1]][pos[2]] -= Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]+1] - volt[0][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]+1][pos[1]][pos[2]] + volt[2][pos[0]][pos[1]][pos[2]]); //for x curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]]; - curr[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]+1][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]+1][pos[2]] + volt[0][pos[0]][pos[1]][pos[2]]); + curr[2][pos[0]][pos[1]][pos[2]] -= Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]+1][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]+1][pos[2]] + volt[0][pos[0]][pos[1]][pos[2]]); } } } diff --git a/FDTD/engine.h b/FDTD/engine.h index 57db067..5fd4dcb 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -5,6 +5,7 @@ class Engine { + friend class ProcessVoltage; public: Engine(Operator* op); virtual ~Engine(); @@ -15,6 +16,8 @@ public: //!Iterate a number of timesteps virtual bool IterateTS(unsigned int iterTS); + unsigned int GetNumberOfTimesteps() {return numTS;}; + protected: Operator* Op; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 5cd6b3e..30d15e0 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -1,3 +1,4 @@ +#include #include "operator.h" #include "tools/array_ops.h" @@ -46,6 +47,37 @@ void Operator::Reset() Operator::Init(); } +unsigned int Operator::GetNyquistNum(double fmax) +{ + if (dT==0) return 1; + double T0 = 1/fmax; + return floor(T0/2/dT); +} + +bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord) +{ + 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;i +#include +#include "operator.h" +#include "engine.h" + +class Processing +{ +public: + Processing(Operator* op, Engine* eng); + virtual ~Processing(); + + virtual void OpenFile(string outfile); + + virtual void Process() {}; + +protected: + Operator* Op; + Engine* Eng; + ofstream file; +}; + +#endif // PROCESSING_H diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp new file mode 100644 index 0000000..5cb5a02 --- /dev/null +++ b/FDTD/processvoltage.cpp @@ -0,0 +1,42 @@ +#include "processvoltage.h" + +ProcessVoltage::ProcessVoltage(Operator* op, Engine* eng) : Processing(op, eng) +{ +} + +ProcessVoltage::~ProcessVoltage() +{ +} + +void ProcessVoltage::DefineStartStopCoord(double* dstart, double* dstop) +{ + 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; +} + +void ProcessVoltage::Process() +{ + FDTD_FLOAT voltage=0; + unsigned int pos[3]={start[0],start[1],start[2]}; +// cerr << Eng->volt[1][pos[0]][pos[1]][pos[2]] << endl; + for (int n=0;n<3;++n) + { + if (start[n]volt[n][pos[0]][pos[1]][pos[2]]; +// cerr << n << " " << pos[0] << " " << pos[1] << " " << pos[2] << " " << Eng->volt[n][pos[0]][pos[1]][pos[2]] << endl; + } + + } + else + { + for (;pos[n]>stop[n];--pos[n]) + voltage-=Eng->volt[n][pos[0]][pos[1]][pos[2]]; + } + } +// cerr << voltage << endl; + voltages.push_back(voltage); + file << Eng->numTS << "\t" << voltage << endl; +} diff --git a/FDTD/processvoltage.h b/FDTD/processvoltage.h new file mode 100644 index 0000000..dccedf4 --- /dev/null +++ b/FDTD/processvoltage.h @@ -0,0 +1,23 @@ +#ifndef PROCESSVOLTAGE_H +#define PROCESSVOLTAGE_H + +#include "processing.h" + +class ProcessVoltage : public Processing +{ +public: + ProcessVoltage(Operator* op, Engine* eng); + virtual ~ProcessVoltage(); + + void DefineStartStopCoord(double* dstart, double* dstop); + + virtual void Process(); + +protected: + unsigned int start[3]; + unsigned int stop[3]; + + vector voltages; +}; + +#endif // PROCESSVOLTAGE_H diff --git a/main.cpp b/main.cpp index 5b03ce2..46c6eaa 100644 --- a/main.cpp +++ b/main.cpp @@ -1,7 +1,11 @@ #include +#include +#include #include +#include "tools/array_ops.h" #include "FDTD/cartoperator.h" #include "FDTD/engine.h" +#include "FDTD/processvoltage.h" #include "ContinuousStructure.h" @@ -11,65 +15,85 @@ void BuildMSL(ContinuousStructure &CSX); int main(int argc, char *argv[]) { - ContinuousStructure CSX; - time_t startTime=time(NULL); + //*************** setup/read geometry ************// + ContinuousStructure CSX; // CSX.ReadFromXML("csx-files/MSL.xml"); BuildMSL(CSX); + //*************** setup operator ************// CartOperator cop; cop.SetGeometryCSX(&CSX); cop.CalcECOperator(); - cop.CalcGaussianPulsExcitation(1e9,1e9); + double fmax=1e9; + cop.CalcGaussianPulsExcitation(fmax/2,fmax/2); time_t OpDoneTime=time(NULL); cop.ShowSize(); +// cop.DumpOperator2File("tmp/Operator"); + + cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl; + unsigned int NrIter = cop.GetNyquistNum(fmax); + cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl; + //create FDTD engine Engine eng(&cop); time_t currTime = time(NULL); - unsigned int NrIter = 500; + //*************** setup processing ************// + ProcessVoltage PV(&cop,&eng); + PV.OpenFile("tmp/u1"); + double start[]={0,0,0}; + double stop[]={0,50,0}; + PV.DefineStartStopCoord(start,stop); + unsigned int maxIter = 5000; - eng.IterateTS(NrIter); + //*************** simulate ************// + for (unsigned int i=0;iSetEpsilon(3.6); - CSX.AddProperty(mat); - - CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat); - box->SetCoord(0,-200.0);box->SetCoord(1,200.0); - box->SetCoord(2,0.0);box->SetCoord(3,50.0); - box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0); - CSX.AddPrimitive(box); - - CSPropMaterial* MSL = new CSPropMaterial(CSX.GetParameterSet()); - MSL->SetKappa(56e6); - CSX.AddProperty(MSL); - - box = new CSPrimBox(CSX.GetParameterSet(),MSL); - box->SetCoord(0,-20.0);box->SetCoord(1,20.0); - box->SetCoord(2,0.0);box->SetCoord(3,50.0); - box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0); - CSX.AddPrimitive(box); +// CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); +// mat->SetEpsilon(3.6); +// CSX.AddProperty(mat); +// +// CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat); +// box->SetCoord(0,-200.0);box->SetCoord(1,200.0); +// box->SetCoord(2,0.0);box->SetCoord(3,50.0); +// box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0); +// CSX.AddPrimitive(box); +// +// CSPropMaterial* MSL = new CSPropMaterial(CSX.GetParameterSet()); +// MSL->SetKappa(56e6); +// CSX.AddProperty(MSL); +// +// box = new CSPrimBox(CSX.GetParameterSet(),MSL); +// box->SetCoord(0,-20.0);box->SetCoord(1,20.0); +// box->SetCoord(2,0.0);box->SetCoord(3,50.0); +// box->SetCoord(4,-1000.0);box->SetCoord(5,1000.0); +// CSX.AddPrimitive(box); CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet()); elec->SetExcitation(1,1); @@ -77,19 +101,19 @@ void BuildMSL(ContinuousStructure &CSX) // elec->SetDelay(2.0e-9); CSX.AddProperty(elec); - box = new CSPrimBox(CSX.GetParameterSet(),elec); - box->SetCoord(0,-20.0);box->SetCoord(1,20.0); + CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec); + box->SetCoord(0,-1.0);box->SetCoord(1,1.0); box->SetCoord(2,0.0);box->SetCoord(3,50.0); - box->SetCoord(4,0.0);box->SetCoord(5,0.0); + box->SetCoord(4,-1.0);box->SetCoord(5,1.0); CSX.AddPrimitive(box); CSRectGrid* grid = CSX.GetGrid(); - for (int n=-1000;n<=1000;n+=20) + for (int n=-100;n<=100;n+=10) grid->AddDiscLine(2,(double)n); - for (int n=-200;n<=200;n+=10) + for (int n=-100;n<=100;n+=10) grid->AddDiscLine(0,(double)n); - for (int n=0;n<=300;n+=10) + for (int n=0;n<=100;n+=10) grid->AddDiscLine(1,(double)n); grid->SetDeltaUnit(1e-3); diff --git a/openEMS.pro b/openEMS.pro index c50de51..e4493f5 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -20,11 +20,15 @@ SOURCES += main.cpp \ tools/AdrOp.cpp \ FDTD/engine.cpp \ FDTD/operator.cpp \ - tools/array_ops.cpp + tools/array_ops.cpp \ + FDTD/processvoltage.cpp \ + FDTD/processing.cpp HEADERS += FDTD/cartoperator.h \ tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ FDTD/engine.h \ FDTD/operator.h \ - tools/array_ops.h + tools/array_ops.h \ + FDTD/processvoltage.h \ + FDTD/processing.h diff --git a/tools/array_ops.cpp b/tools/array_ops.cpp index 637708a..945a8c2 100644 --- a/tools/array_ops.cpp +++ b/tools/array_ops.cpp @@ -1,4 +1,5 @@ #include "array_ops.h" +#include FDTD_FLOAT*** Create3DArray(unsigned int* numLines) { @@ -56,3 +57,21 @@ void Delete_N_3DArray(FDTD_FLOAT**** array, unsigned int* numLines) } delete[] array; } + +void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, unsigned int* numLines) +{ + unsigned int pos[3]; + for (pos[0]=0;pos[0]