zcy
/
openEMS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Engine + BC bug fix

pull/1/head
Thorsten Liebig 2010-03-02 19:01:03 +01:00
parent 6d2e974cc1
commit 1c92ab2920
6 changed files with 91 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
FDTD/cartoperator.cpp
View File

@ -133,11 +133,15 @@ void CartOperator::ApplyElectricBC(bool* dirs)
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP]) for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
{ {
pos[n]=0; pos[n]=0;
vv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; vv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
vi[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; vi[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
vv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
vi[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
pos[n]=numLines[n]-1; pos[n]=numLines[n]-1;
vv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; vv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
vi[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; vi[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
vv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
vi[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
} }
} }
} }
@ -157,12 +161,16 @@ void CartOperator::ApplyMagneticBC(bool* dirs)
for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP]) for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
{ {
pos[n]=0; pos[n]=0;
ii[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; ii[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
iv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; iv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
ii[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
iv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n];
pos[n]=numLines[n]-2; pos[n]=numLines[n]-2;
ii[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; ii[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
iv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; iv[nP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
ii[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
iv[nPP][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1];
} }
} }
} }

16
FDTD/engine.cpp
View File

@ -32,28 +32,32 @@ bool Engine::IterateTS(unsigned int iterTS)
{ {
unsigned int pos[3]; unsigned int pos[3];
int exc_pos; int exc_pos;
bool shift[3];
for (unsigned int iter=0;iter<iterTS;++iter) for (unsigned int iter=0;iter<iterTS;++iter)
{ {
//voltage updates //voltage updates
for (pos[0]=1;pos[0]<Op->numLines[0]-1;++pos[0]) for (pos[0]=0;pos[0]<Op->numLines[0];++pos[0])
{ {
for (pos[1]=1;pos[1]<Op->numLines[1]-1;++pos[1]) shift[0]=pos[0];
for (pos[1]=0;pos[1]<Op->numLines[1];++pos[1])
{ {
for (pos[2]=1;pos[2]<Op->numLines[2]-1;++pos[2]) shift[1]=pos[1];
for (pos[2]=0;pos[2]<Op->numLines[2];++pos[2])
{ {
shift[2]=pos[2];
//do the updates here //do the updates here
//for x //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->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 //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->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 //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->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]]);
} }
} }
} }

1
FDTD/processfields_td.cpp
View File

@ -14,6 +14,7 @@ ProcessFieldsTD::~ProcessFieldsTD()
void ProcessFieldsTD::Process() void ProcessFieldsTD::Process()
{ {
if (Enabled==false) return;
stringstream ss; stringstream ss;
ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->numTS; ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->numTS;

7
FDTD/processing.h
View File

@ -15,11 +15,18 @@ public:
virtual void Process() {}; 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: protected:
Processing(Operator* op, Engine* eng); Processing(Operator* op, Engine* eng);
Operator* Op; Operator* Op;
Engine* Eng; Engine* Eng;
bool Enabled;
unsigned int start[3]; unsigned int start[3];
unsigned int stop[3]; unsigned int stop[3];
}; };

1
FDTD/processvoltage.cpp
View File

@ -23,6 +23,7 @@ void ProcessVoltage::OpenFile(string outfile)
void ProcessVoltage::Process() void ProcessVoltage::Process()
{ {
if (Enabled==false) return;
FDTD_FLOAT voltage=0; FDTD_FLOAT voltage=0;
unsigned int pos[3]={start[0],start[1],start[2]}; unsigned int pos[3]={start[0],start[1],start[2]};
// cerr << Eng->volt[1][pos[0]][pos[1]][pos[2]] << endl; // cerr << Eng->volt[1][pos[0]][pos[1]][pos[2]] << endl;

64
main.cpp
View File

@ -14,6 +14,7 @@ using namespace std;
void BuildMSL(ContinuousStructure &CSX); void BuildMSL(ContinuousStructure &CSX);
void BuildDipol(ContinuousStructure &CSX); void BuildDipol(ContinuousStructure &CSX);
void BuildPlaneWave(ContinuousStructure &CSX);
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
@ -21,7 +22,7 @@ int main(int argc, char *argv[])
//*************** setup/read geometry ************// //*************** setup/read geometry ************//
ContinuousStructure CSX; ContinuousStructure CSX;
BuildDipol(CSX); BuildPlaneWave(CSX);
//*************** setup operator ************// //*************** setup operator ************//
CartOperator cop; CartOperator cop;
@ -30,11 +31,13 @@ int main(int argc, char *argv[])
cop.CalcECOperator(); cop.CalcECOperator();
double fmax=1e9; double fmax=1e9;
cop.CalcGaussianPulsExcitation(fmax/2,fmax/2); cop.CalcGaussianPulsExcitation(fmax/2,fmax/8);
time_t OpDoneTime=time(NULL); time_t OpDoneTime=time(NULL);
cop.ShowSize(); cop.ShowSize();
bool bnd[] = {1,1,0,0,0,0};
cop.ApplyMagneticBC(bnd);
cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl; cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl;
unsigned int NrIter = cop.GetNyquistNum(fmax)/3; unsigned int NrIter = cop.GetNyquistNum(fmax)/3;
@ -49,14 +52,16 @@ int main(int argc, char *argv[])
//*************** setup processing ************// //*************** setup processing ************//
ProcessVoltage PV(&cop,&eng); ProcessVoltage PV(&cop,&eng);
PV.OpenFile("tmp/u1"); PV.OpenFile("tmp/u1");
double start[]={-0,-75,0}; double start[]={-100,-250,0};
double stop[]={-0,75,0}; double stop[]={-100,250,0};
PV.DefineStartStopCoord(start,stop); PV.DefineStartStopCoord(start,stop);
unsigned int maxIter = 5000; unsigned int maxIter = 2000;
ProcessFieldsTD PETD(&cop,&eng); ProcessFieldsTD PETD(&cop,&eng);
start[0]=-1000;start[1]=0;start[2]=-1000; start[0]=-250;start[1]=0;start[2]=-2000;
stop[0]=1000;stop[1]=0;stop[2]=1000; 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.SetDumpType(0);
PETD.SetFilePattern("tmp/Et_"); PETD.SetFilePattern("tmp/Et_");
PETD.DefineStartStopCoord(start,stop); PETD.DefineStartStopCoord(start,stop);
@ -68,9 +73,12 @@ int main(int argc, char *argv[])
PHTD.SetFilePattern("tmp/Ht_"); PHTD.SetFilePattern("tmp/Ht_");
PHTD.DefineStartStopCoord(start,stop); PHTD.DefineStartStopCoord(start,stop);
PETD.SetEnable(false);
PHTD.SetEnable(false);
PV.Process(); PV.Process();
PETD.Process(); PETD.Process();
PHTD.Process(); PHTD.Process();
//*************** simulate ************// //*************** simulate ************//
for (unsigned int i=0;i<maxIter;i+=NrIter) for (unsigned int i=0;i<maxIter;i+=NrIter)
@ -132,6 +140,46 @@ void BuildDipol(ContinuousStructure &CSX)
CSX.Write2XML("tmp/Dipol.xml"); CSX.Write2XML("tmp/Dipol.xml");
} }
void BuildPlaneWave(ContinuousStructure &CSX)
{
// CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());
//// mat->SetKappa(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) void BuildMSL(ContinuousStructure &CSX)
{ {
// CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); // CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet());