diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 96b0f9a..459022a 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -65,15 +65,15 @@ bool Operator::SnapToMesh(double* dcoord, unsigned int* uicoord, bool lower) else if (dcoord[n]>discLines[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; + if (Op->SnapToMesh(dstop,stop,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Snapping error in mesh, check stop value!!" << endl; for (int n=0;n<3;++n) { diff --git a/FDTD/processfields_td.cpp b/FDTD/processfields_td.cpp index 5f5dec7..130cd89 100644 --- a/FDTD/processfields_td.cpp +++ b/FDTD/processfields_td.cpp @@ -41,15 +41,15 @@ void ProcessFieldsTD::Process() //in x delta = Op->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; + E_T[0][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); //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; + E_T[1][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); //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; + E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); } } } @@ -57,5 +57,43 @@ void ProcessFieldsTD::Process() Delete_N_3DArray(E_T,numDLines); E_T = NULL; } + + if (DumpType==1) + { + //create array + FDTD_FLOAT**** H_T = Create_N_3DArray(numDLines); + unsigned int pos[3] = {start[0],start[1],start[2]}; + unsigned int OpPos[3]; + double delta; +// cerr << "processing h-fields... " << endl; + for (pos[0]=0;pos[0]discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; + else delta = 0.5* (Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]-1]); + H_T[0][pos[0]][pos[1]][pos[2]] = Eng->curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[0][OpPos[0]+1][OpPos[1]][OpPos[2]]; + H_T[0][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); + //in y + delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; + H_T[1][pos[0]][pos[1]][pos[2]] = Eng->curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[0][OpPos[0]][OpPos[1]+1][OpPos[2]]; + H_T[1][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); + //in z + delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; + H_T[2][pos[0]][pos[1]][pos[2]] = Eng->curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[0][OpPos[0]][OpPos[1]][OpPos[2]+1]; + H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); + } + } + } + DumpFieldArray2VTK(file,string("H-Field"),H_T,discDLines,numDLines); + Delete_N_3DArray(H_T,numDLines); + H_T = NULL; + } file.close(); } diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp index 4aa4c90..6bf4faf 100644 --- a/FDTD/processing.cpp +++ b/FDTD/processing.cpp @@ -14,6 +14,6 @@ Processing::~Processing() void Processing::DefineStartStopCoord(double* dstart, double* dstop) { 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; + if (Op->SnapToMesh(dstop,stop)==false) cerr << "Processing::DefineStartStopCoord: Snapping error in mesh, check stop value!!" << endl; } diff --git a/main.cpp b/main.cpp index 9221cd9..3d4e704 100644 --- a/main.cpp +++ b/main.cpp @@ -57,11 +57,20 @@ int main(int argc, char *argv[]) ProcessFieldsTD PETD(&cop,&eng); start[0]=-1000;start[1]=0;start[2]=-1000; stop[0]=1000;stop[1]=0;stop[2]=1000; + PETD.SetDumpType(0); PETD.SetFilePattern("tmp/Et_"); PETD.DefineStartStopCoord(start,stop); + ProcessFieldsTD PHTD(&cop,&eng); + start[0]=-1000;start[2]=0;start[1]=-1000; + stop[0]=1000;stop[2]=0;stop[1]=1000; + PHTD.SetDumpType(1); + PHTD.SetFilePattern("tmp/Ht_"); + PHTD.DefineStartStopCoord(start,stop); + PV.Process(); PETD.Process(); + PHTD.Process(); //*************** simulate ************// for (unsigned int i=0;iSetCoord(0,-1.0);box->SetCoord(1,1.0); + box->SetCoord(0,-10.0);box->SetCoord(1,10.0); box->SetCoord(2,-75.0);box->SetCoord(3,75.0); - box->SetCoord(4,-1.0);box->SetCoord(5,1.0); + box->SetCoord(4,-10.0);box->SetCoord(5,10.0); CSX.AddPrimitive(box); CSRectGrid* grid = CSX.GetGrid(); - for (int n=-1000;n<=1000;n+=20) + for (int n=-990;n<=990;n+=20) grid->AddDiscLine(2,(double)n); - for (int n=-1000;n<=1000;n+=20) + for (int n=-990;n<=990;n+=20) grid->AddDiscLine(0,(double)n); - for (int n=-1000;n<=1000;n+=20) + for (int n=-990;n<=990;n+=20) grid->AddDiscLine(1,(double)n); grid->SetDeltaUnit(1e-3);