diff --git a/FDTD/engine.h b/FDTD/engine.h index 23c2d24..a917686 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -22,11 +22,6 @@ class Engine { - friend class Processing; - friend class ProcessVoltage; - friend class ProcessCurrent; - friend class ProcessFields; - friend class ProcessFieldsTD; public: static Engine* createEngine(Operator* op); virtual ~Engine(); @@ -39,6 +34,9 @@ public: virtual unsigned int GetNumberOfTimesteps() {return numTS;}; + virtual FDTD_FLOAT**** GetVoltages() {return volt;}; + virtual FDTD_FLOAT**** GetCurrents() {return curr;}; + protected: Engine(Operator* op); Operator* Op; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 71d0ea0..f7b5588 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -106,7 +106,7 @@ 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=1;i=0) + if (((int)currPos[n]-1)>=0) { P[n] = discLines[n][currPos[n]-1]; Point_Line_Distance(P,start,stop,foot,dist); @@ -169,11 +167,7 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[]) minFoot=foot; minDist=dist; minDir = n; - minPos[0]=currPos[0]; - minPos[1]=currPos[1]; - minPos[2]=currPos[2]; - minPos[n]=currPos[n]-1; -// array[n][minPos[0]][minPos[1]][minPos[2]] = 1; + UpDir = false; } } if ((currPos[n]+1) vDelay; vector vDir; unsigned int ipos; - int pos[3]; + unsigned int pos[3]; double coord[3]; double delta[3]; double amp=0; @@ -782,6 +785,70 @@ bool Operator::CalcEFieldExcitation() } } + //special treatment for primitives of type curve (treated as wires) see also Calc_PEC + double p1[3]; + double p2[3]; + double deltaN=0.0; + int n; + struct Grid_Path path; + CSPropElectrode* elec=NULL; + CSProperties* prop=NULL; + vector vec_prop = CSX->GetPropertyByType(CSProperties::ELECTRODE); + for (size_t p=0;pToElectrode(); + for (size_t n=0;nGetQtyPrimitives();++n) + { + CSPrimitives* prim = prop->GetPrimitive(n); + CSPrimCurve* curv = prim->ToCurve(); + if (curv) + { + for (size_t i=1;iGetNumberOfPoints();++i) + { + curv->GetPoint(i-1,p1); + curv->GetPoint(i,p2); + path = FindPath(p1,p2); + for (size_t t=0;tSetPos(pos[0],pos[1],pos[2]); + deltaN=fabs(MainOp->GetIndexDelta(n,pos[n])); + coord[0] = discLines[0][pos[0]]; + coord[1] = discLines[1][pos[1]]; + coord[2] = discLines[2][pos[2]]; + coord[n] += 0.5*deltaN; +// cerr << n << " " << coord[0] << " " << coord[1] << " " << coord[2] << endl; + if (elec!=NULL) + { + if ((elec->GetActiveDir(n)) && (pos[n]GetWeightedExcitation(n,coord)*deltaN*gridDelta; + if (amp!=0) + { + vExcit.push_back(amp); + vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); + vDir.push_back(n); + vIndex[0].push_back(pos[0]); + vIndex[1].push_back(pos[1]); + vIndex[2].push_back(pos[2]); + } + if (elec->GetExcitType()==1) //hard excite + { + vv[n][pos[0]][pos[1]][pos[2]] = 0; + vi[n][pos[0]][pos[1]][pos[2]] = 0; + } + } + } + } + } + } + } + } + E_Exc_Count = vExcit.size(); cerr << "Operator::CalcEFieldExcitation: Found number of excitations points: " << E_Exc_Count << endl; if (E_Exc_Count==0) @@ -862,24 +929,18 @@ bool Operator::CalcPEC() curv->GetPoint(i-1,p1); curv->GetPoint(i,p2); path = FindPath(p1,p2); -// cerr << p1[0] << " " << p1[1] << " " << p1[2] << endl; -// cerr << p2[0] << " " << p2[1] << " " << p2[2] << endl; for (size_t t=0;tGetCurrents(); if (Enabled==false) return -1; if (CheckTimestep()==false) return GetNextInterval(); FDTD_FLOAT current=0; @@ -65,26 +66,26 @@ int ProcessCurrent::Process() } //x-current - for (int i=start[0];icurr[0][i][start[1]][start[2]]; + for (unsigned int i=start[0];icurr[1][stop[0]][i][start[2]]; + for (unsigned int i=start[1];icurr[2][stop[0]][stop[1]][i]; + for (unsigned int i=start[2];icurr[0][i][stop[1]][stop[2]]; + for (unsigned int i=start[0];icurr[1][start[0]][i][stop[2]]; + for (unsigned int i=start[1];icurr[2][start[0]][start[1]][i]; + for (unsigned int i=start[2];inumTS << " i: " << current << endl; v_current.push_back(current); - file << (double)Eng->numTS*Op->GetTimestep() << "\t" << current << endl; + file << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << current << endl; return GetNextInterval(); } diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp index 49214c7..ff70810 100644 --- a/FDTD/processfields.cpp +++ b/FDTD/processfields.cpp @@ -101,6 +101,8 @@ void ProcessFields::DefineStartStopCoord(double* dstart, double* dstop) double ProcessFields::CalcTotalEnergy() { + FDTD_FLOAT**** volt = Eng->GetVoltages(); + FDTD_FLOAT**** curr = Eng->GetCurrents(); double energy=0; if (Eng==NULL) return 0.0; unsigned int pos[3]; @@ -110,12 +112,12 @@ double ProcessFields::CalcTotalEnergy() { for (pos[2]=0;pos[2]numLines[2];++pos[2]) { - energy+=fabs(Eng->volt[0][pos[0]][pos[1]][pos[2]] * Eng->curr[1][pos[0]][pos[1]][pos[2]]); - energy+=fabs(Eng->volt[0][pos[0]][pos[1]][pos[2]] * Eng->curr[2][pos[0]][pos[1]][pos[2]]); - energy+=fabs(Eng->volt[1][pos[0]][pos[1]][pos[2]] * Eng->curr[0][pos[0]][pos[1]][pos[2]]); - energy+=fabs(Eng->volt[1][pos[0]][pos[1]][pos[2]] * Eng->curr[2][pos[0]][pos[1]][pos[2]]); - energy+=fabs(Eng->volt[2][pos[0]][pos[1]][pos[2]] * Eng->curr[0][pos[0]][pos[1]][pos[2]]); - energy+=fabs(Eng->volt[2][pos[0]][pos[1]][pos[2]] * Eng->curr[1][pos[0]][pos[1]][pos[2]]); + energy+=fabs(volt[0][pos[0]][pos[1]][pos[2]] * curr[1][pos[0]][pos[1]][pos[2]]); + energy+=fabs(volt[0][pos[0]][pos[1]][pos[2]] * curr[2][pos[0]][pos[1]][pos[2]]); + energy+=fabs(volt[1][pos[0]][pos[1]][pos[2]] * curr[0][pos[0]][pos[1]][pos[2]]); + energy+=fabs(volt[1][pos[0]][pos[1]][pos[2]] * curr[2][pos[0]][pos[1]][pos[2]]); + energy+=fabs(volt[2][pos[0]][pos[1]][pos[2]] * curr[0][pos[0]][pos[1]][pos[2]]); + energy+=fabs(volt[2][pos[0]][pos[1]][pos[2]] * curr[1][pos[0]][pos[1]][pos[2]]); } } } @@ -183,16 +185,18 @@ bool ProcessFields::DumpVectorArray2VTK(ofstream &file, string name, FDTD_FLOAT* { WriteVTKHeader(file, discLines, numLines); WriteVTKVectorArray(file, name, array, numLines); + return true; } bool ProcessFields::DumpMultiVectorArray2VTK(ofstream &file, string names[], FDTD_FLOAT**** array[], unsigned int numFields, double** discLines, unsigned int* numLines) { WriteVTKHeader(file, discLines, numLines); - for (int n=0;nGetVoltages(); + FDTD_FLOAT**** curr = Eng->GetCurrents(); + if (DumpType==0) { //create array @@ -50,15 +53,15 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file) OpPos[2]=start[2]+pos[2]; //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]] = volt[0][OpPos[0]][OpPos[1]][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[0][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[0][OpPos[0]][OpPos[1]+1][OpPos[2]+1]; 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]] = volt[1][OpPos[0]][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[1][OpPos[0]][OpPos[1]][OpPos[2]+1] + volt[1][OpPos[0]+1][OpPos[1]][OpPos[2]+1]; 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]] = volt[2][OpPos[0]][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]][OpPos[1]+1][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]][OpPos[2]] + volt[2][OpPos[0]+1][OpPos[1]+1][OpPos[2]]; E_T[2][pos[0]][pos[1]][pos[2]] /= (4*delta*Op->gridDelta); } } @@ -88,15 +91,15 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file) //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]] = curr[0][OpPos[0]][OpPos[1]][OpPos[2]] + 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[1][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[1][OpPos[0]][OpPos[1]+1][OpPos[2]]; + H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][OpPos[0]][OpPos[1]][OpPos[2]] + curr[1][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[2][OpPos[0]][OpPos[1]][OpPos[2]] + Eng->curr[2][OpPos[0]][OpPos[1]][OpPos[2]+1]; + H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][OpPos[0]][OpPos[1]][OpPos[2]] + curr[2][OpPos[0]][OpPos[1]][OpPos[2]+1]; H_T[2][pos[0]][pos[1]][pos[2]] /= (2*delta*Op->gridDelta); } } @@ -109,6 +112,9 @@ void ProcessFieldsTD::DumpCellInterpol(ofstream &file) void ProcessFieldsTD::DumpNoInterpol(ofstream &file) { + FDTD_FLOAT**** volt = Eng->GetVoltages(); + FDTD_FLOAT**** curr = Eng->GetCurrents(); + unsigned int pos[3]; double delta[3]; if (DumpType==0) @@ -124,9 +130,9 @@ void ProcessFieldsTD::DumpNoInterpol(ofstream &file) for (pos[2]=0;pos[2]MainOp->GetIndexDelta(2,pos[2]+start[2])); - E_T[0][pos[0]][pos[1]][pos[2]] = Eng->volt[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta; - E_T[1][pos[0]][pos[1]][pos[2]] = Eng->volt[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta; - E_T[2][pos[0]][pos[1]][pos[2]] = Eng->volt[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta; + E_T[0][pos[0]][pos[1]][pos[2]] = volt[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta; + E_T[1][pos[0]][pos[1]][pos[2]] = volt[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta; + E_T[2][pos[0]][pos[1]][pos[2]] = volt[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta; } } } @@ -149,9 +155,9 @@ void ProcessFieldsTD::DumpNoInterpol(ofstream &file) { delta[2]=fabs(Op->MainOp->GetIndexWidth(2,pos[2]+start[2])); //in x - H_T[0][pos[0]][pos[1]][pos[2]] = Eng->curr[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta; - H_T[1][pos[0]][pos[1]][pos[2]] = Eng->curr[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta; - H_T[2][pos[0]][pos[1]][pos[2]] = Eng->curr[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta; + H_T[0][pos[0]][pos[1]][pos[2]] = curr[0][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[0]/Op->gridDelta; + H_T[1][pos[0]][pos[1]][pos[2]] = curr[1][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[1]/Op->gridDelta; + H_T[2][pos[0]][pos[1]][pos[2]] = curr[2][pos[0]+start[0]][pos[1]+start[1]][pos[2]+start[2]]/delta[2]/Op->gridDelta; } } } @@ -167,7 +173,7 @@ int ProcessFieldsTD::Process() if (filePattern.empty()) return -1; if (CheckTimestep()==false) return GetNextInterval(); stringstream ss; - ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->numTS; + ss << std::setw( pad_length ) << std::setfill( '0' ) << Eng->GetNumberOfTimesteps(); string filename = filePattern + ss.str() + ".vtk"; ofstream file(filename.c_str()); diff --git a/FDTD/processing.cpp b/FDTD/processing.cpp index 5be1f68..a54dcca 100644 --- a/FDTD/processing.cpp +++ b/FDTD/processing.cpp @@ -51,9 +51,9 @@ double Processing::CalcLineIntegral(unsigned int* start, unsigned int* stop, int double result=0; FDTD_FLOAT**** array; if (field==0) - array=Eng->volt; + array=Eng->GetVoltages(); else if (field==1) - array=Eng->curr; + array=Eng->GetCurrents(); else return 0.0; for (int n=0;n<3;++n) diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp index 25fe31a..64d41b0 100644 --- a/FDTD/processvoltage.cpp +++ b/FDTD/processvoltage.cpp @@ -45,6 +45,6 @@ int ProcessVoltage::Process() FDTD_FLOAT voltage=CalcLineIntegral(start,stop,0); // cerr << voltage << endl; voltages.push_back(voltage); - file << (double)Eng->numTS*Op->GetTimestep() << "\t" << voltage << endl; + file << (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() << "\t" << voltage << endl; return GetNextInterval(); } diff --git a/openems.cpp b/openems.cpp index 921ae40..d97346a 100644 --- a/openems.cpp +++ b/openems.cpp @@ -138,7 +138,12 @@ int openEMS::SetupFDTD(const char* file) cerr << "Can't read openEMS FDTD Settings... " << endl; exit(-1); } - FDTD_Opts->QueryIntAttribute("NumberOfTimesteps",&NrTS); + int help=0; + FDTD_Opts->QueryIntAttribute("NumberOfTimesteps",&help); + if (help<0) + NrTS=0; + else + NrTS = help; FDTD_Opts->QueryDoubleAttribute("endCriteria",&endCrit); if (endCrit==0) endCrit=1e-6; @@ -192,26 +197,39 @@ int openEMS::SetupFDTD(const char* file) FDTD_Op = new Operator(); if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1); - if (DebugMat) - { - FDTD_Op->DumpMaterial2File("material_dump.vtk"); - } FDTD_Op->CalcECOperator(); - if (DebugOp) - { - FDTD_Op->DumpOperator2File("operator_dump.vtk"); - } if (Excit_Type==0) - FDTD_Op->CalcGaussianPulsExcitation(f0,fc); + { + if (!FDTD_Op->CalcGaussianPulsExcitation(f0,fc)) + { + cerr << "openEMS: excitation setup failed!!" << endl; + exit(2); + } + } else if (Excit_Type==1) - FDTD_Op->CalcSinusExcitation(f0,NrTS); + { + if (!FDTD_Op->CalcSinusExcitation(f0,NrTS)) + { + cerr << "openEMS: excitation setup failed!!" << endl; + exit(2); + } + } else { cerr << "openEMS: Excitation type is unknown" << endl; exit(-1); } + if (DebugMat) + { + FDTD_Op->DumpMaterial2File("material_dump.vtk"); + } + if (DebugOp) + { + FDTD_Op->DumpOperator2File("operator_dump.vtk"); + } + time_t OpDoneTime=time(NULL); FDTD_Op->ShowSize(); @@ -330,14 +348,14 @@ void openEMS::RunFDTD() //*************** simulate ************// int step=PA->Process(); - if ((step<0) || (step>NrTS)) step=NrTS; + if ((step<0) || (step>(int)NrTS)) step=NrTS; while ((FDTD_Eng->GetNumberOfTimesteps()endCrit)) { FDTD_Eng->IterateTS(step); step=PA->Process(); // cout << " do " << step << " steps; current: " << eng.GetNumberOfTimesteps() << endl; currTS = FDTD_Eng->GetNumberOfTimesteps(); - if ((step<0) || (step>NrTS - currTS)) step=NrTS - currTS; + if ((step<0) || (step>(int)(NrTS - currTS))) step=NrTS - currTS; gettimeofday(&currTime,NULL); diff --git a/openems.h b/openems.h index ae2885d..d218c7a 100644 --- a/openems.h +++ b/openems.h @@ -44,7 +44,7 @@ public: protected: //! Number of Timesteps - int NrTS; + unsigned int NrTS; bool Enable_Dumps; bool DebugMat; bool DebugOp;