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

fix in snapping, E-Dump and new H-Dump

pull/1/head
Thorsten Liebig 2010-03-02 15:37:00 +01:00
parent ca44334da5
commit 6d2e974cc1
5 changed files with 63 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
FDTD/operator.cpp
View File

@ -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]) {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 if (dcoord[n]==discLines[n][numLines[n]-1]) {uicoord[n]=numLines[n]-1; if (lower) uicoord[n]=numLines[n]-2;}
else else
for (unsigned int i=0;i<numLines[n]-1;++i) for (unsigned int i=1;i<numLines[n]-1;++i)
{ {
if (dcoord[n]<=discLines[n][i]) if (dcoord[n]<discLines[n][i])
{ {
if (fabs(dcoord[n]-discLines[n][i])<(fabs(dcoord[n]-discLines[n][i+1]))) if (fabs(dcoord[n]-discLines[n][i])<(fabs(dcoord[n]-discLines[n][i-1])))
uicoord[n]=i; uicoord[n]=i;
else else
uicoord[n]=i+1; uicoord[n]=i-1;
if (lower) uicoord[n]=i; if (lower) uicoord[n]=i-1;
i = numLines[n]; i = numLines[n];
} }
} }

2
FDTD/processfields.cpp
View File

@ -23,7 +23,7 @@ ProcessFields::~ProcessFields()
void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop)
{ {
if (Op->SnapToMesh(dstart,start,true)==false) cerr << "ProcessFields::DefineStartStopCoord: Snapping error in mesh, check start value!!" << endl; if (Op->SnapToMesh(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) for (int n=0;n<3;++n)
{ {

44
FDTD/processfields_td.cpp
View File

@ -41,15 +41,15 @@ void ProcessFieldsTD::Process()
//in x //in x
delta = Op->discLines[0][OpPos[0]+1] - Op->discLines[0][OpPos[0]]; 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]] = 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 //in y
delta = Op->discLines[1][OpPos[1]+1] - Op->discLines[1][OpPos[1]]; 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]] = 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 //in z
delta = Op->discLines[2][OpPos[2]+1] - Op->discLines[2][OpPos[2]]; 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]] = 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); Delete_N_3DArray(E_T,numDLines);
E_T = NULL; 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]<numDLines[0];++pos[0])
{
OpPos[0]=start[0]+pos[0];
for (pos[1]=0;pos[1]<numDLines[1];++pos[1])
{
OpPos[1]=start[1]+pos[1];
for (pos[2]=0;pos[2]<numDLines[2];++pos[2])
{
OpPos[2]=start[2]+pos[2];
//in x
if (OpPos[0]==0) delta = Op->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(); file.close();
} }

2
FDTD/processing.cpp
View File

@ -14,6 +14,6 @@ Processing::~Processing()
void Processing::DefineStartStopCoord(double* dstart, double* dstop) 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(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;
} }

20
main.cpp
View File

@ -57,11 +57,20 @@ int main(int argc, char *argv[])
ProcessFieldsTD PETD(&cop,&eng); ProcessFieldsTD PETD(&cop,&eng);
start[0]=-1000;start[1]=0;start[2]=-1000; start[0]=-1000;start[1]=0;start[2]=-1000;
stop[0]=1000;stop[1]=0;stop[2]=1000; stop[0]=1000;stop[1]=0;stop[2]=1000;
PETD.SetDumpType(0);
PETD.SetFilePattern("tmp/Et_"); PETD.SetFilePattern("tmp/Et_");
PETD.DefineStartStopCoord(start,stop); 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(); PV.Process();
PETD.Process(); PETD.Process();
PHTD.Process();
//*************** simulate ************// //*************** simulate ************//
for (unsigned int i=0;i<maxIter;i+=NrIter) for (unsigned int i=0;i<maxIter;i+=NrIter)
@ -69,6 +78,7 @@ int main(int argc, char *argv[])
eng.IterateTS(NrIter); eng.IterateTS(NrIter);
PV.Process(); PV.Process();
PETD.Process(); PETD.Process();
PHTD.Process();
} }
//*************** postproc ************// //*************** postproc ************//
@ -103,18 +113,18 @@ void BuildDipol(ContinuousStructure &CSX)
CSX.AddProperty(elec); CSX.AddProperty(elec);
CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec); CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec);
box->SetCoord(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(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); CSX.AddPrimitive(box);
CSRectGrid* grid = CSX.GetGrid(); 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); 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); 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->AddDiscLine(1,(double)n);
grid->SetDeltaUnit(1e-3); grid->SetDeltaUnit(1e-3);