diff --git a/FDTD/cartoperator.cpp b/FDTD/cartoperator.cpp index 8502401..3a4efec 100644 --- a/FDTD/cartoperator.cpp +++ b/FDTD/cartoperator.cpp @@ -133,11 +133,15 @@ void CartOperator::ApplyElectricBC(bool* dirs) for (pos[nPP]=0;pos[nPP]numLines[0]-1;++pos[0]) + for (pos[0]=0;pos[0]numLines[0];++pos[0]) { - for (pos[1]=1;pos[1]numLines[1]-1;++pos[1]) + shift[0]=pos[0]; + for (pos[1]=0;pos[1]numLines[1];++pos[1]) { - for (pos[2]=1;pos[2]numLines[2]-1;++pos[2]) + shift[1]=pos[1]; + for (pos[2]=0;pos[2]numLines[2];++pos[2]) { + shift[2]=pos[2]; //do the updates here //for x volt[0][pos[0]][pos[1]][pos[2]] *= Op->vv[0][pos[0]][pos[1]][pos[2]]; - volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-1][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-1]); + volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-shift[1]][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-shift[2]]); //for y volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]]; - volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-1] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-1][pos[1]][pos[2]]); + volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-shift[0]][pos[1]][pos[2]]); //for x volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]]; - volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-1][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-1][pos[2]]); + volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-shift[1]][pos[2]]); } } } diff --git a/FDTD/processfields_td.cpp b/FDTD/processfields_td.cpp index 130cd89..9d5f87c 100644 --- a/FDTD/processfields_td.cpp +++ b/FDTD/processfields_td.cpp @@ -14,6 +14,7 @@ ProcessFieldsTD::~ProcessFieldsTD() void ProcessFieldsTD::Process() { + if (Enabled==false) return; stringstream ss; ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->numTS; diff --git a/FDTD/processing.h b/FDTD/processing.h index 08b7c72..40a2a43 100644 --- a/FDTD/processing.h +++ b/FDTD/processing.h @@ -15,11 +15,18 @@ public: virtual void Process() {}; + //! If Disabled Process() will do nothing... + virtual void SetEnable(bool val) {Enabled=val;} + //! If Disabled Process() will do nothing... + virtual bool GetEnable() {return Enabled;} + protected: Processing(Operator* op, Engine* eng); Operator* Op; Engine* Eng; + bool Enabled; + unsigned int start[3]; unsigned int stop[3]; }; diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp index f9444fb..711568f 100644 --- a/FDTD/processvoltage.cpp +++ b/FDTD/processvoltage.cpp @@ -23,6 +23,7 @@ void ProcessVoltage::OpenFile(string outfile) void ProcessVoltage::Process() { + if (Enabled==false) return; 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; diff --git a/main.cpp b/main.cpp index 3d4e704..53fcf5a 100644 --- a/main.cpp +++ b/main.cpp @@ -14,6 +14,7 @@ using namespace std; void BuildMSL(ContinuousStructure &CSX); void BuildDipol(ContinuousStructure &CSX); +void BuildPlaneWave(ContinuousStructure &CSX); int main(int argc, char *argv[]) { @@ -21,7 +22,7 @@ int main(int argc, char *argv[]) //*************** setup/read geometry ************// ContinuousStructure CSX; - BuildDipol(CSX); + BuildPlaneWave(CSX); //*************** setup operator ************// CartOperator cop; @@ -30,11 +31,13 @@ int main(int argc, char *argv[]) cop.CalcECOperator(); double fmax=1e9; - cop.CalcGaussianPulsExcitation(fmax/2,fmax/2); + cop.CalcGaussianPulsExcitation(fmax/2,fmax/8); time_t OpDoneTime=time(NULL); cop.ShowSize(); + bool bnd[] = {1,1,0,0,0,0}; + cop.ApplyMagneticBC(bnd); cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl; unsigned int NrIter = cop.GetNyquistNum(fmax)/3; @@ -49,14 +52,16 @@ int main(int argc, char *argv[]) //*************** setup processing ************// ProcessVoltage PV(&cop,&eng); PV.OpenFile("tmp/u1"); - double start[]={-0,-75,0}; - double stop[]={-0,75,0}; + double start[]={-100,-250,0}; + double stop[]={-100,250,0}; PV.DefineStartStopCoord(start,stop); - unsigned int maxIter = 5000; + unsigned int maxIter = 2000; ProcessFieldsTD PETD(&cop,&eng); - start[0]=-1000;start[1]=0;start[2]=-1000; - stop[0]=1000;stop[1]=0;stop[2]=1000; + start[0]=-250;start[1]=0;start[2]=-2000; + stop[0]=250;stop[1]=0;stop[2]=4000; +// start[0]=-250;start[1]=-250;start[2]=0; +// stop[0]=250;stop[1]=250;stop[2]=0; PETD.SetDumpType(0); PETD.SetFilePattern("tmp/Et_"); PETD.DefineStartStopCoord(start,stop); @@ -68,9 +73,12 @@ int main(int argc, char *argv[]) PHTD.SetFilePattern("tmp/Ht_"); PHTD.DefineStartStopCoord(start,stop); + PETD.SetEnable(false); + PHTD.SetEnable(false); + PV.Process(); PETD.Process(); - PHTD.Process(); + PHTD.Process(); //*************** simulate ************// for (unsigned int i=0;iSetKappa(0.001); +// CSX.AddProperty(mat); +// +// CSPrimBox* matbox = new CSPrimBox(CSX.GetParameterSet(),mat); +// matbox->SetCoord(0,-1000.0);matbox->SetCoord(1,1000.0); +// matbox->SetCoord(2,-1000.0);matbox->SetCoord(3,1000.0); +// matbox->SetCoord(4,-2000.0);matbox->SetCoord(5,4000.0); +// CSX.AddPrimitive(matbox); + + CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet()); + elec->SetExcitation(1,1); + elec->SetExcitType(0); + elec->SetActiveDir(0,0);//disable x + elec->SetActiveDir(0,2);//disable z +// elec->SetDelay(2.0e-9); + CSX.AddProperty(elec); + + CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec); + box->SetCoord(0,-250.0);box->SetCoord(1,250.0); + box->SetCoord(2,-250.0);box->SetCoord(3,250.0); + box->SetCoord(4,-0000.0);box->SetCoord(5,-0000.0); + CSX.AddPrimitive(box); + + CSRectGrid* grid = CSX.GetGrid(); + + for (int n=-250;n<=250;n+=10) + grid->AddDiscLine(0,(double)n); + for (int n=-250;n<=250;n+=10) + grid->AddDiscLine(1,(double)n); + for (int n=-4000;n<=4000;n+=10) + grid->AddDiscLine(2,(double)n); + + grid->SetDeltaUnit(1e-3); + + CSX.Write2XML("tmp/PlaneWave.xml"); +} + void BuildMSL(ContinuousStructure &CSX) { // CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());