diff --git a/Analyse/PlotVoltage.m b/Analyse/PlotVoltage.m index d2ad72f..02cec09 100644 --- a/Analyse/PlotVoltage.m +++ b/Analyse/PlotVoltage.m @@ -2,22 +2,62 @@ clear all; clc -tmp = load('../tmp/u1'); +figure(1); +tmpu = load('../tmp/u1'); +tmpi = load('../tmp/i1'); -t = tmp(:,1); -u = tmp(:,2); +t = tmpu(:,1); +u = tmpu(:,2); -L=numel(t); - -subplot(2,1,1); +subplot(2,2,1); +title('u_1 TD'); plot(t,u); +xlabel('t \rightarrow'); +ylabel('ut_1 \rightarrow'); grid on; dt=t(2)-t(1); +u= [u ; zeros(size(u))]; +L=numel(u); +t = (1:L)*dt; f = (1:L)/L/dt; fu = fft(u)/L; -subplot(2,1,2); +subplot(2,2,2); +title('u_1 FD'); plot(f(1:L/2),abs(fu(1:L/2))); +xlabel('f \rightarrow'); +ylabel('|uf_1| \rightarrow'); grid on; + +t = tmpi(:,1); +i = tmpi(:,2); + +subplot(2,2,3); +title('i_1 TD'); +plot(t,i); +xlabel('t \rightarrow'); +ylabel('it_1 \rightarrow'); +grid on; + +dt=t(2)-t(1); +i = [i; zeros(size(t))]; +L=numel(i); +t = (1:L)*dt; +f = (1:L)/L/dt; + +fi = fft(i)/L; +subplot(2,2,4); +title('i_1 FD'); +plot(f(1:L/2),abs(fi(1:L/2))); +xlabel('f \rightarrow'); +ylabel('|if_1| \rightarrow'); +grid on; + +figure(2); +plot(f,real(fu./fi)); +xlim([0 1e9]); +grid on; + + diff --git a/FDTD/cartoperator.cpp b/FDTD/cartoperator.cpp index f3ab6a6..e359227 100644 --- a/FDTD/cartoperator.cpp +++ b/FDTD/cartoperator.cpp @@ -204,8 +204,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC) if (prop) { CSPropMaterial* mat = prop->ToMaterial(); - inEC[0] = mat->GetEpsilon(n)*fabs(deltaP*deltaPP); - inEC[1] = mat->GetKappa(n)*fabs(deltaP*deltaPP); + inEC[0] = mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); + inEC[1] = mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); } else { @@ -220,8 +220,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC) if (prop) { CSPropMaterial* mat = prop->ToMaterial(); - inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP); - inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP); + inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); + inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); } else { @@ -237,8 +237,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC) if (prop) { CSPropMaterial* mat = prop->ToMaterial(); - inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP); - inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP); + inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); + inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); } else { @@ -254,8 +254,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC) if (prop) { CSPropMaterial* mat = prop->ToMaterial(); - inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP); - inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP); + inEC[0] += mat->GetEpsilonWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); + inEC[1] += mat->GetKappaWeighted(n,shiftCoord)*fabs(deltaP*deltaPP); } else { @@ -275,9 +275,9 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC) if (prop) { CSPropMaterial* mat = prop->ToMaterial(); - inEC[2] = fabs(delta_M) / mat->GetMue(n); + inEC[2] = fabs(delta_M) / mat->GetMueWeighted(n,shiftCoord); if (mat->GetSigma(n)) - inEC[3] = fabs(delta_M) / mat->GetSigma(n); + inEC[3] = fabs(delta_M) / mat->GetSigmaWeighted(n,shiftCoord); else inEC[3] = 0; } @@ -295,8 +295,8 @@ bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC) { CSPropMaterial* mat = prop->ToMaterial(); inEC[2] += mat->GetMue(n)*fabs(delta); - if (mat->GetSigma(n)) - inEC[3] += fabs(delta)/mat->GetSigma(n); + if (mat->GetSigmaWeighted(n,shiftCoord)) + inEC[3] += fabs(delta)/mat->GetSigmaWeighted(n,shiftCoord); else inEC[3] = 0; } @@ -454,8 +454,10 @@ bool CartOperator::CalcEFieldExcitation() 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; + vv[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0; + vi[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0; + vv[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0; + vi[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0; } } } diff --git a/FDTD/engine.h b/FDTD/engine.h index 44c92eb..5a8da91 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -6,6 +6,7 @@ class Engine { friend class ProcessVoltage; + friend class ProcessCurrent; friend class ProcessFieldsTD; public: Engine(Operator* op); diff --git a/FDTD/processcurrent.cpp b/FDTD/processcurrent.cpp new file mode 100644 index 0000000..ad74233 --- /dev/null +++ b/FDTD/processcurrent.cpp @@ -0,0 +1,71 @@ +#include "processcurrent.h" + +ProcessCurrent::ProcessCurrent(Operator* op, Engine* eng) : Processing(op, eng) +{ +} + +ProcessCurrent::~ProcessCurrent() +{ +} + +void ProcessCurrent::OpenFile(string outfile) +{ + 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 ProcessCurrent::DefineStartStopCoord(double* dstart, double* dstop) +{ + if (Op->SnapToMesh(dstart,start,true)==false) cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapping problem, check start value!!" << endl; + if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessCurrent::DefineStartStopCoord: Warning: Snapping problem, check stop value!!" << endl; +} + +void ProcessCurrent::Process() +{ + if (Enabled==false) return; + FDTD_FLOAT current=0; +// FDTD_FLOAT help=0; + double sign[3]={1,1,1}; + unsigned int pos[3]={start[0],start[1],start[2]}; + double loc_start[]={start[0],start[1],start[2]}; + double loc_stop[]={stop[0],stop[1],stop[2]}; + + for (int n=0;n<3;++n) + { + if (start[n]>stop[n]) + { + unsigned int help=start[n]; + start[n]=stop[n]; + stop[n]=help; + } + } + + //x-current + for (int i=start[0];icurr[0][i][start[1]][start[2]]; + //y-current + for (int i=start[1];icurr[1][stop[0]][i][start[2]]; + //z-current + for (int i=start[2];icurr[2][stop[0]][stop[1]][i]; + //x-current + for (int i=start[0];icurr[0][i][stop[1]][stop[2]]; + //y-current + for (int i=start[1];icurr[1][start[0]][i][stop[2]]; + //z-current + for (int i=start[2];icurr[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; +} diff --git a/FDTD/processcurrent.h b/FDTD/processcurrent.h new file mode 100644 index 0000000..1910994 --- /dev/null +++ b/FDTD/processcurrent.h @@ -0,0 +1,24 @@ +#ifndef PROCESSCURRENT_H +#define PROCESSCURRENT_H + +#include "processing.h" + +class ProcessCurrent : public Processing +{ +public: + ProcessCurrent(Operator* op, Engine* eng); + virtual ~ProcessCurrent(); + + virtual void OpenFile(string outfile); + + virtual void DefineStartStopCoord(double* dstart, double* dstop); + + virtual void Process(); + +protected: + ofstream file; + + vector v_current; +}; + +#endif // PROCESSCURRENT_H diff --git a/main.cpp b/main.cpp index d372e71..0f5e9a4 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/processcurrent.h" #include "FDTD/processfields_td.h" #include "ContinuousStructure.h" @@ -21,23 +22,26 @@ int main(int argc, char *argv[]) time_t startTime=time(NULL); //*************** setup/read geometry ************// + cerr << "Create Geometry..." << endl; ContinuousStructure CSX; - BuildPlaneWave(CSX); + //BuildPlaneWave(CSX); + BuildMSL(CSX); //*************** setup operator ************// + cerr << "Create Operator..." << endl; CartOperator cop; cop.SetGeometryCSX(&CSX); cop.CalcECOperator(); double fmax=1e9; - cop.CalcGaussianPulsExcitation(fmax/2,fmax/8); + cop.CalcGaussianPulsExcitation(fmax/2,fmax/2); time_t OpDoneTime=time(NULL); cop.ShowSize(); - bool bnd[] = {1,1,0,0,0,0}; + bool bnd[] = {1,1,0,1,0,1}; cop.ApplyMagneticBC(bnd); cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl; @@ -53,16 +57,23 @@ int main(int argc, char *argv[]) //*************** setup processing ************// ProcessVoltage PV(&cop,&eng); PV.OpenFile("tmp/u1"); - double start[]={-100,-250,0}; - double stop[]={-100,250,0}; + double start[]={0,50,0}; + double stop[]={0,0,0}; PV.DefineStartStopCoord(start,stop); - unsigned int maxIter = 2000; + + ProcessCurrent PCurr(&cop,&eng); + PCurr.OpenFile("tmp/i1"); + start[0]=-50;start[1]=40;start[2]=-0; + stop[0]=50;stop[1]=60;stop[2]=-0; + PCurr.DefineStartStopCoord(start,stop); + + unsigned int maxIter = 1800; ProcessFieldsTD PETD(&cop,&eng); - start[0]=-480;start[1]=0;start[2]=-4000; - stop[0]=480;stop[1]=0;stop[2]=4000; -// start[0]=-250;start[1]=-250;start[2]=0; -// stop[0]=250;stop[1]=250;stop[2]=0; + start[0]=-500;start[1]=25;start[2]=-800; + stop[0] = 500;stop[1] =24;stop[2] = 800; +// start[0]=-300;start[1]=0;start[2]=-4000; +// stop[0]=300;stop[1]=300;stop[2]=-4000; PETD.SetDumpType(0); PETD.SetFilePattern("tmp/Et_"); PETD.DefineStartStopCoord(start,stop); @@ -78,6 +89,7 @@ int main(int argc, char *argv[]) PHTD.SetEnable(false); PV.Process(); + PCurr.Process(); PETD.Process(); PHTD.Process(); @@ -86,6 +98,7 @@ int main(int argc, char *argv[]) { eng.IterateTS(NrIter); PV.Process(); + PCurr.Process(); PETD.Process(); PHTD.Process(); } @@ -97,7 +110,7 @@ int main(int argc, char *argv[]) double t_diff = difftime(currTime,prevTime); cerr << "Time for " << eng.GetNumberOfTimesteps() << " iterations with " << cop.GetNumberCells() << " cells : " << t_diff << " sec" << endl; - cerr << "Speed (MCells/s): " << (double)cop.GetNumberCells()*(double)eng.GetNumberOfTimesteps()/t_diff/1e6 << endl; + cerr << "Speed: " << (double)cop.GetNumberCells()*(double)eng.GetNumberOfTimesteps()/t_diff/1e6 << " MCells/s " << endl; return 0; } @@ -192,46 +205,73 @@ void BuildPlaneWave(ContinuousStructure &CSX) void BuildMSL(ContinuousStructure &CSX) { -// 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); + double width = 1000; + double hight = 500; + double length = 2000; + double abs_l = 200; + //substrate.... + CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); +// mat->SetEpsilon(3.6); + double finalKappa = 3/pow(abs_l,4); + mat->SetKappa(finalKappa); + std::ostringstream fct; + fct << "pow(abs(z)-" << length/2.0-abs_l << ",4)"; + mat->SetKappaWeightFunction(fct.str(),0); + mat->SetKappaWeightFunction(fct.str(),1); + mat->SetKappaWeightFunction(fct.str(),2); + mat->SetSigma(finalKappa*__MUE0__/__EPS0__); + mat->SetSigmaWeightFunction(fct.str(),0); + mat->SetSigmaWeightFunction(fct.str(),1); + mat->SetSigmaWeightFunction(fct.str(),2); + CSX.AddProperty(mat); + CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat); + box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); + box->SetCoord(2,0.0);box->SetCoord(3,hight); + box->SetCoord(4,length/2.0-abs_l); box->SetCoord(5,length/2.0); + CSX.AddPrimitive(box); + box = new CSPrimBox(CSX.GetParameterSet(),mat); + box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); + box->SetCoord(2,0.0);box->SetCoord(3,hight); + box->SetCoord(4,length/-2.0+abs_l); box->SetCoord(5,length/-2.0); + CSX.AddPrimitive(box); + +// double coords[]={0,0,length/2}; +// cerr << fct.str() << endl; +// cerr << finalKappa*pow(abs_l,4) << endl; +// cerr << mat->GetKappaWeighted(0,coords) << endl; +// exit(0); + + //MSL + CSPropMetal* MSL = new CSPropMetal(CSX.GetParameterSet()); +// MSL->SetKappa(56e6); + CSX.AddProperty(MSL); + box = new CSPrimBox(CSX.GetParameterSet(),MSL); + box->SetCoord(0,-40.0);box->SetCoord(1,40.0); + box->SetCoord(2,50.0);box->SetCoord(3,50.0); + box->SetCoord(4,length/-2);box->SetCoord(5,length/2.0); + CSX.AddPrimitive(box); + + //MSL excite... CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet()); - elec->SetExcitation(1,1); - elec->SetExcitType(0); + elec->SetExcitation(-1,1); + elec->SetExcitType(1); // elec->SetDelay(2.0e-9); CSX.AddProperty(elec); - - CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec); - box->SetCoord(0,-1.0);box->SetCoord(1,1.0); + box = new CSPrimBox(CSX.GetParameterSet(),elec); + box->SetCoord(0,-40.0);box->SetCoord(1,40.0); box->SetCoord(2,0.0);box->SetCoord(3,50.0); - box->SetCoord(4,-1.0);box->SetCoord(5,1.0); + box->SetCoord(4,length/-2.0);box->SetCoord(5,length/-2.0); CSX.AddPrimitive(box); CSRectGrid* grid = CSX.GetGrid(); - for (int n=-300;n<=300;n+=20) - grid->AddDiscLine(0,(double)n); - for (int n=0;n<=300;n+=20) - grid->AddDiscLine(1,(double)n); - for (int n=-1000;n<=1000;n+=20) - grid->AddDiscLine(2,(double)n); + for (double n=width/-2.0;n<=width/2;n+=20) + grid->AddDiscLine(0,n); + for (double n=0;n<=500;n+=10) + grid->AddDiscLine(1,n); + for (double n=length/-2.0;n<=length/2.0;n+=20) + grid->AddDiscLine(2,n); grid->SetDeltaUnit(1e-3); diff --git a/openEMS.pro b/openEMS.pro index d4ac0d6..35fc409 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -7,7 +7,8 @@ CONFIG += console CONFIG -= app_bundle TEMPLATE = app OBJECTS_DIR = obj -INCLUDEPATH += ../CSXCAD +INCLUDEPATH += ../CSXCAD \ + ../fparser LIBS += -L../CSXCAD \ -lCSXCAD \ -L../fparser \ @@ -24,7 +25,8 @@ SOURCES += main.cpp \ FDTD/processvoltage.cpp \ FDTD/processing.cpp \ FDTD/processfields.cpp \ - FDTD/processfields_td.cpp + FDTD/processfields_td.cpp \ + FDTD/processcurrent.cpp HEADERS += FDTD/cartoperator.h \ tools/ErrorMsg.h \ tools/AdrOp.h \ @@ -35,4 +37,5 @@ HEADERS += FDTD/cartoperator.h \ FDTD/processvoltage.h \ FDTD/processing.h \ FDTD/processfields.h \ - FDTD/processfields_td.h + FDTD/processfields_td.h \ + FDTD/processcurrent.h diff --git a/tools/constants.h b/tools/constants.h index ecc8ec1..8fecc4f 100644 --- a/tools/constants.h +++ b/tools/constants.h @@ -3,6 +3,8 @@ #define __EPS0__ 8.85418781762e-12 #define __MUE0__ 1.256637062e-6 +#define __C0__ 299792458 +#define __Z0__ 376.730313461 #define PI 3.141592653589793238462643383279 #endif // CONSTANTS_H