diff --git a/examples/Dipol.xml b/examples/Dipol.xml new file mode 100644 index 0000000..24cdd7c --- /dev/null +++ b/examples/Dipol.xml @@ -0,0 +1,41 @@ + + + + + + + + -990,-970,-950,-930,-910,-890,-870,-850,-830,-810,-790,-770,-750,-730,-710,-690,-670,-650,-630,-610,-590,-570,-550,-530,-510,-490,-470,-450,-430,-410,-390,-370,-350,-330,-310,-290,-270,-250,-230,-210,-190,-170,-150,-130,-110,-90,-70,-50,-30,-10,10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330,350,370,390,410,430,450,470,490,510,530,550,570,590,610,630,650,670,690,710,730,750,770,790,810,830,850,870,890,910,930,950,970,990 + -990,-970,-950,-930,-910,-890,-870,-850,-830,-810,-790,-770,-750,-730,-710,-690,-670,-650,-630,-610,-590,-570,-550,-530,-510,-490,-470,-450,-430,-410,-390,-370,-350,-330,-310,-290,-270,-250,-230,-210,-190,-170,-150,-130,-110,-90,-70,-50,-30,-10,10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330,350,370,390,410,430,450,470,490,510,530,550,570,590,610,630,650,670,690,710,730,750,770,790,810,830,850,870,890,910,930,950,970,990 + -990,-970,-950,-930,-910,-890,-870,-850,-830,-810,-790,-770,-750,-730,-710,-690,-670,-650,-630,-610,-590,-570,-550,-530,-510,-490,-470,-450,-430,-410,-390,-370,-350,-330,-310,-290,-270,-250,-230,-210,-190,-170,-150,-130,-110,-90,-70,-50,-30,-10,10,30,50,70,90,110,130,150,170,190,210,230,250,270,290,310,330,350,370,390,410,430,450,470,490,510,530,550,570,590,610,630,650,670,690,710,730,750,770,790,810,830,850,870,890,910,930,950,970,990 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/FDTD_examples.cpp b/examples/FDTD_examples.cpp index 10fdc93..12f262a 100644 --- a/examples/FDTD_examples.cpp +++ b/examples/FDTD_examples.cpp @@ -1,8 +1,16 @@ #include "FDTD_examples.h" #include "../tools/constants.h" -void BuildDipol(ContinuousStructure &CSX) +void BuildDipol(const char* filename) { + int maxIter = 1000; + double fmax=1e9; + int Excit_Type=0; + int bounds[] = {1,1,0,0,0,0}; + + cerr << "Create Geometry..." << endl; + ContinuousStructure CSX; + CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); mat->SetKappa(0.001); CSX.AddProperty(mat); @@ -38,11 +46,48 @@ void BuildDipol(ContinuousStructure &CSX) grid->SetDeltaUnit(1e-3); - CSX.Write2XML("tmp/Dipol.xml"); + //*************** Create XML file ********************** + TiXmlDocument doc(filename); + doc.InsertEndChild(TiXmlDeclaration("1.0","ISO-8859-1","yes")); + + TiXmlElement FDTD_Opts("openEMS-Parameter"); + FDTD_Opts.SetAttribute("NumberOfTimesteps",maxIter); + + TiXmlElement Excite("Excitation"); + Excite.SetAttribute("Type",Excit_Type); + Excite.SetAttribute("f0",fmax); + FDTD_Opts.InsertEndChild(Excite); + + TiXmlElement BC("BoundaryCond"); + BC.SetAttribute("xmin",bounds[0]); + BC.SetAttribute("xmax",bounds[1]); + BC.SetAttribute("ymin",bounds[2]); + BC.SetAttribute("ymax",bounds[3]); + BC.SetAttribute("zmin",bounds[4]); + BC.SetAttribute("zmax",bounds[5]); + FDTD_Opts.InsertEndChild(BC); + + doc.InsertEndChild(FDTD_Opts); + + if (CSX.Write2XML(&doc,true)==false) + { + cerr << "writing failed" << endl; + exit(-1); + } + + doc.SaveFile(); } -void BuildPlaneWave(ContinuousStructure &CSX) +void BuildPlaneWave(const char* filename) { + int maxIter = 1000; + double fmax=1e9; + int Excit_Type=0; + int bounds[] = {1,1,0,0,0,0}; + + cerr << "Create Geometry..." << endl; + ContinuousStructure CSX; + double width = 1000; double hight = 1000; double length = 4000; @@ -140,15 +185,55 @@ void BuildPlaneWave(ContinuousStructure &CSX) grid->SetDeltaUnit(1e-3); - CSX.Write2XML("tmp/PlaneWave.xml"); + //*************** Create XML file ********************** + TiXmlDocument doc(filename); + doc.InsertEndChild(TiXmlDeclaration("1.0","ISO-8859-1","yes")); + + TiXmlElement FDTD_Opts("openEMS-Parameter"); + FDTD_Opts.SetAttribute("NumberOfTimesteps",maxIter); + + TiXmlElement Excite("Excitation"); + Excite.SetAttribute("Type",Excit_Type); + Excite.SetAttribute("f0",fmax); + FDTD_Opts.InsertEndChild(Excite); + + TiXmlElement BC("BoundaryCond"); + BC.SetAttribute("xmin",bounds[0]); + BC.SetAttribute("xmax",bounds[1]); + BC.SetAttribute("ymin",bounds[2]); + BC.SetAttribute("ymax",bounds[3]); + BC.SetAttribute("zmin",bounds[4]); + BC.SetAttribute("zmax",bounds[5]); + FDTD_Opts.InsertEndChild(BC); + + doc.InsertEndChild(FDTD_Opts); + + if (CSX.Write2XML(&doc,true)==false) + { + cerr << "writing failed" << endl; + exit(-1); + } + + doc.SaveFile(); } -void BuildMSL(ContinuousStructure &CSX) +void BuildMSL(const char* filename) { + int maxIter = 1000; + double fmax=1e9; + int Excit_Type=0; + int bounds[] = {1,1,0,0,0,0}; + + cerr << "Create Geometry..." << endl; + ContinuousStructure CSX; + double width = 1000; - double hight = 500; + double height = 500; double length = 2000; double abs_l = 200; + double MSL_height=50; + double MSL_width=80; + double delta[] = {20,10,20}; //substrate.... CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); @@ -167,7 +252,7 @@ void BuildMSL(ContinuousStructure &CSX) CSX.AddProperty(mat); CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),mat); box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); - box->SetCoord(2,0.0);box->SetCoord(3,hight); + box->SetCoord(2,0.0);box->SetCoord(3,height); box->SetCoord(4,length/2.0-abs_l); box->SetCoord(5,length/2.0); CSX.AddPrimitive(box); // box = new CSPrimBox(CSX.GetParameterSet(),mat); @@ -177,11 +262,16 @@ void BuildMSL(ContinuousStructure &CSX) // CSX.AddPrimitive(box); //MSL - CSPropMetal* MSL = new CSPropMetal(CSX.GetParameterSet()); + CSProperties* MSL = NULL; + CSPropMaterial* MSL_mat = new CSPropMaterial(CSX.GetParameterSet()); + MSL_mat->SetKappa(56e6); + MSL = MSL_mat; +// MSL = new CSPropMetal(CSX.GetParameterSet()); + CSX.AddProperty(MSL); box = new CSPrimBox(CSX.GetParameterSet(),MSL); - box->SetCoord(0,-40.0);box->SetCoord(1,40.0); - box->SetCoord(2,50.0);box->SetCoord(3,50.0); + box->SetCoord(0,MSL_width/-2.0);box->SetCoord(1,MSL_width/2.0); + box->SetCoord(2,MSL_height);box->SetCoord(3,MSL_height+delta[1]); box->SetCoord(4,length/-2);box->SetCoord(5,length/2.0); box->SetPriority(100); CSX.AddPrimitive(box); @@ -198,29 +288,104 @@ void BuildMSL(ContinuousStructure &CSX) box->SetCoord(4,length/-2.0);box->SetCoord(5,length/-2.0); CSX.AddPrimitive(box); - //E-field dump - CSPropDumpBox* Edump = new CSPropDumpBox(CSX.GetParameterSet()); + CSPropDumpBox* Edump = NULL; +// //E-field dump xz +// Edump = new CSPropDumpBox(CSX.GetParameterSet()); +// Edump->SetDumpType(0); +// Edump->SetName("tmp/Et_xz_"); +// CSX.AddProperty(Edump); +// box = new CSPrimBox(CSX.GetParameterSet(),Edump); +// box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); +// box->SetCoord(2,25.0);box->SetCoord(3,25.); +// box->SetCoord(4,length/-2.0);box->SetCoord(5,length/2.0); +// CSX.AddPrimitive(box); +// +// //E-field dump xy +// Edump = new CSPropDumpBox(CSX.GetParameterSet()); +// Edump->SetDumpType(0); +// Edump->SetName("tmp/Et_xy_"); +// CSX.AddProperty(Edump); +// box = new CSPrimBox(CSX.GetParameterSet(),Edump); +// box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); +// box->SetCoord(2,0.0);box->SetCoord(3,height); +// box->SetCoord(4,0.0);box->SetCoord(5,0.0); +// CSX.AddPrimitive(box); + + //E-field dump 3D + Edump = new CSPropDumpBox(CSX.GetParameterSet()); Edump->SetDumpType(0); + Edump->SetDumpMode(2); //cell interpolated dump Edump->SetName("tmp/Et_"); CSX.AddProperty(Edump); box = new CSPrimBox(CSX.GetParameterSet(),Edump); - box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); - box->SetCoord(2,25.0);box->SetCoord(3,25.); + box->SetCoord(0,MSL_width*-1.5);box->SetCoord(1,MSL_width*1.5); + box->SetCoord(2,0.0);box->SetCoord(3,MSL_height*1.5); box->SetCoord(4,length/-2.0);box->SetCoord(5,length/2.0); CSX.AddPrimitive(box); + //voltage calc + CSPropProbeBox* volt = new CSPropProbeBox(CSX.GetParameterSet()); + volt->SetProbeType(0); + volt->SetName("tmp/u1"); + CSX.AddProperty(volt); + box = new CSPrimBox(CSX.GetParameterSet(),volt); + box->SetCoord(0,0.0);box->SetCoord(1,0.0); + box->SetCoord(2,MSL_height);box->SetCoord(3,0.0); + box->SetCoord(4,0.0);box->SetCoord(5,0.0); + CSX.AddPrimitive(box); + + //current calc + CSPropProbeBox* curr = new CSPropProbeBox(CSX.GetParameterSet()); + curr->SetProbeType(1); + curr->SetName("tmp/i1"); + CSX.AddProperty(curr); + box = new CSPrimBox(CSX.GetParameterSet(),curr); + box->SetCoord(0,MSL_width*-1.5);box->SetCoord(1,MSL_width*1.5); + box->SetCoord(2,MSL_height/2.0);box->SetCoord(3,MSL_height*1.5); + box->SetCoord(4,0.0);box->SetCoord(5,0.0); + CSX.AddPrimitive(box); + CSRectGrid* grid = CSX.GetGrid(); - for (double n=width/-2.0;n<=width/2;n+=20) + for (double n=width/-2.0;n<=width/2;n+=delta[0]) grid->AddDiscLine(0,n); - for (double n=0;n<=500;n+=10) + for (double n=0;n<=height;n+=delta[1]) grid->AddDiscLine(1,n); - for (double n=length/-2.0;n<=length/2.0;n+=20) + for (double n=length/-2.0;n<=length/2.0;n+=delta[2]) grid->AddDiscLine(2,n); grid->SetDeltaUnit(1e-3); - CSX.Write2XML("tmp/MSL.xml"); + //*************** Create XML file ********************** + TiXmlDocument doc(filename); + doc.InsertEndChild(TiXmlDeclaration("1.0","ISO-8859-1","yes")); + + TiXmlElement FDTD_Opts("openEMS-Parameter"); + FDTD_Opts.SetAttribute("NumberOfTimesteps",maxIter); + + TiXmlElement Excite("Excitation"); + Excite.SetAttribute("Type",Excit_Type); + Excite.SetAttribute("f0",fmax); + FDTD_Opts.InsertEndChild(Excite); + + TiXmlElement BC("BoundaryCond"); + BC.SetAttribute("xmin",bounds[0]); + BC.SetAttribute("xmax",bounds[1]); + BC.SetAttribute("ymin",bounds[2]); + BC.SetAttribute("ymax",bounds[3]); + BC.SetAttribute("zmin",bounds[4]); + BC.SetAttribute("zmax",bounds[5]); + FDTD_Opts.InsertEndChild(BC); + + doc.InsertEndChild(FDTD_Opts); + + if (CSX.Write2XML(&doc,true)==false) + { + cerr << "writing failed" << endl; + exit(-1); + } + + doc.SaveFile(); } diff --git a/examples/FDTD_examples.h b/examples/FDTD_examples.h index e4737ff..69b0cc6 100644 --- a/examples/FDTD_examples.h +++ b/examples/FDTD_examples.h @@ -2,11 +2,12 @@ #define FDTD_EXAMPLES_H #include "ContinuousStructure.h" +#include "tinyxml.h" -void BuildDipol(ContinuousStructure &CSX); +void BuildDipol(const char* filename); -void BuildPlaneWave(ContinuousStructure &CSX); +void BuildPlaneWave(const char* filename); -void BuildMSL(ContinuousStructure &CSX); +void BuildMSL(const char* filename); #endif // FDTD_EXAMPLES_H diff --git a/examples/MSL.xml b/examples/MSL.xml new file mode 100644 index 0000000..fea7642 --- /dev/null +++ b/examples/MSL.xml @@ -0,0 +1,82 @@ + + + + + + + + -500,-480,-460,-440,-420,-400,-380,-360,-340,-320,-300,-280,-260,-240,-220,-200,-180,-160,-140,-120,-100,-80,-60,-40,-20,0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360,380,400,420,440,460,480,500 + 0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240,250,260,270,280,290,300,310,320,330,340,350,360,370,380,390,400,410,420,430,440,450,460,470,480,490,500 + -1000,-980,-960,-940,-920,-900,-880,-860,-840,-820,-800,-780,-760,-740,-720,-700,-680,-660,-640,-620,-600,-580,-560,-540,-520,-500,-480,-460,-440,-420,-400,-380,-360,-340,-320,-300,-280,-260,-240,-220,-200,-180,-160,-140,-120,-100,-80,-60,-40,-20,0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360,380,400,420,440,460,480,500,520,540,560,580,600,620,640,660,680,700,720,740,760,780,800,820,840,860,880,900,920,940,960,980,1000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/PlaneWave.xml b/examples/PlaneWave.xml new file mode 100644 index 0000000..854d8e1 --- /dev/null +++ b/examples/PlaneWave.xml @@ -0,0 +1,81 @@ + + + + + + + + -500,-480,-460,-440,-420,-400,-380,-360,-340,-320,-300,-280,-260,-240,-220,-200,-180,-160,-140,-120,-100,-80,-60,-40,-20,0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360,380,400,420,440,460,480,500 + -500,-480,-460,-440,-420,-400,-380,-360,-340,-320,-300,-280,-260,-240,-220,-200,-180,-160,-140,-120,-100,-80,-60,-40,-20,0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360,380,400,420,440,460,480,500 + -2000,-1980,-1960,-1940,-1920,-1900,-1880,-1860,-1840,-1820,-1800,-1780,-1760,-1740,-1720,-1700,-1680,-1660,-1640,-1620,-1600,-1580,-1560,-1540,-1520,-1500,-1480,-1460,-1440,-1420,-1400,-1380,-1360,-1340,-1320,-1300,-1280,-1260,-1240,-1220,-1200,-1180,-1160,-1140,-1120,-1100,-1080,-1060,-1040,-1020,-1000,-980,-960,-940,-920,-900,-880,-860,-840,-820,-800,-780,-760,-740,-720,-700,-680,-660,-640,-620,-600,-580,-560,-540,-520,-500,-480,-460,-440,-420,-400,-380,-360,-340,-320,-300,-280,-260,-240,-220,-200,-180,-160,-140,-120,-100,-80,-60,-40,-20,0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360,380,400,420,440,460,480,500,520,540,560,580,600,620,640,660,680,700,720,740,760,780,800,820,840,860,880,900,920,940,960,980,1000,1020,1040,1060,1080,1100,1120,1140,1160,1180,1200,1220,1240,1260,1280,1300,1320,1340,1360,1380,1400,1420,1440,1460,1480,1500,1520,1540,1560,1580,1600,1620,1640,1660,1680,1700,1720,1740,1760,1780,1800,1820,1840,1860,1880,1900,1920,1940,1960,1980,2000 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/main.cpp b/main.cpp index b759e90..2b39372 100644 --- a/main.cpp +++ b/main.cpp @@ -1,133 +1,46 @@ #include #include #include -#include -#include "tools/array_ops.h" -#include "FDTD/operator.h" -#include "FDTD/engine.h" -#include "FDTD/processvoltage.h" -#include "FDTD/processcurrent.h" -#include "FDTD/processfields_td.h" -#include "ContinuousStructure.h" +#include "openems.h" #include "examples/FDTD_examples.h" +//#define STANDALONE + using namespace std; int main(int argc, char *argv[]) { - time_t startTime=time(NULL); +#ifdef STANDALONE + if (argc<=1) + { + cerr << " argc= " << argc << endl; + exit(-1); + } + char* file = argv[1]; +#else //*************** setup/read geometry ************// - cerr << "Create Geometry..." << endl; - ContinuousStructure CSX; - BuildPlaneWave(CSX); -// BuildMSL(CSX); + const char* fileDP="examples/Dipol.xml"; + BuildDipol(fileDP); - //*************** setup operator ************// - cerr << "Create Operator..." << endl; - Operator cop; - cop.SetGeometryCSX(&CSX); + const char* filePW="examples/PlaneWave.xml"; + BuildPlaneWave(filePW); - cop.CalcECOperator(); + const char* fileMSL="examples/MSL.xml"; + BuildMSL(fileMSL); - double fmax=1e9; - cop.CalcGaussianPulsExcitation(fmax/2,fmax/2); + const char* file=fileMSL; - time_t OpDoneTime=time(NULL); +// cerr << CSX.ReadFromXML("examples/PlaneWave.xml") << endl; +#endif + openEMS FDTD; - cop.ShowSize(); - bool bnd[] = {1,1,0,0,0,0}; - cop.ApplyMagneticBC(bnd); + int EC = FDTD.SetupFDTD(file); + if (EC) return EC; + FDTD.RunFDTD(); - cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl; - unsigned int Nyquist = cop.GetNyquistNum(fmax); - - cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl; - - //create FDTD engine - Engine eng(&cop); - - time_t currTime = time(NULL); - - //*************** setup processing ************// - ProcessingArray PA; - - ProcessVoltage* PV = new ProcessVoltage(&cop,&eng); - PA.AddProcessing(PV); - PV->SetProcessInterval(Nyquist/3); //three times as often as nyquist dictates - PV->OpenFile("tmp/u1"); - double start[]={0,50,0}; - double stop[]={0,0,0}; - PV->DefineStartStopCoord(start,stop); - - ProcessCurrent* PCurr = new ProcessCurrent(&cop,&eng); - PA.AddProcessing(PCurr); - PCurr->SetProcessInterval(Nyquist/3); //three times as often as nyquist dictates - PCurr->OpenFile("tmp/i1"); - start[0]=-50;start[1]=40;start[2]=-0; - stop[0]=50;stop[1]=60;stop[2]=-0; - PCurr->DefineStartStopCoord(start,stop); - - unsigned int maxIter = 1800; - - bool EnableDump = true; - - vector DumpProps = CSX.GetPropertyByType(CSProperties::DUMPBOX); - for (size_t i=0;iSetEnable(EnableDump); - ProcTD->SetProcessInterval(Nyquist/2); //twice as often as nyquist dictates - - //only looking for one prim atm - CSPrimitives* prim = DumpProps.at(i)->GetPrimitive(0); - if (prim==NULL) - delete ProcTD; - else - { - bool acc; - double* bnd = prim->GetBoundBox(acc); - start[0]= bnd[0];start[1]=bnd[2];start[2]=bnd[4]; - stop[0] = bnd[1];stop[1] =bnd[3];stop[2] =bnd[5]; - ProcTD->DefineStartStopCoord(start,stop); - CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox(); - if (db) - { - ProcTD->SetDumpType(db->GetDumpType()); - ProcTD->SetDumpMode(db->GetDumpMode()); - ProcTD->SetFilePattern(db->GetName()); - PA.AddProcessing(ProcTD); - } - else - delete ProcTD; - } - } - - int step=PA.Process(); - if ((step<0) || (step>maxIter)) step=maxIter; - //*************** simulate ************// - while (eng.GetNumberOfTimesteps()maxIter - eng.GetNumberOfTimesteps())) step=maxIter - eng.GetNumberOfTimesteps(); - } - - //*************** postproc ************// - time_t prevTime = currTime; - currTime = time(NULL); - - double t_diff = difftime(currTime,prevTime); - - cerr << "Time for " << eng.GetNumberOfTimesteps() << " iterations with " << cop.GetNumberCells() << " cells : " << t_diff << " sec" << endl; - cerr << "Speed: " << (double)cop.GetNumberCells()*(double)eng.GetNumberOfTimesteps()/t_diff/1e6 << " MCells/s " << endl; - - //cleanup - PA.DeleteAll(); - - return 0; + return 0; } diff --git a/openEMS.pro b/openEMS.pro index 5e1cb52..cbc89fa 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -8,7 +8,8 @@ CONFIG -= app_bundle TEMPLATE = app OBJECTS_DIR = obj INCLUDEPATH += ../CSXCAD \ - ../fparser + ../fparser \ + ../tinyxml LIBS += -L../CSXCAD \ -lCSXCAD \ -L../fparser \ @@ -26,7 +27,8 @@ SOURCES += main.cpp \ FDTD/processfields.cpp \ FDTD/processfields_td.cpp \ FDTD/processcurrent.cpp \ - examples/FDTD_examples.cpp + examples/FDTD_examples.cpp \ + openems.cpp HEADERS += tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ @@ -38,4 +40,5 @@ HEADERS += tools/ErrorMsg.h \ FDTD/processfields.h \ FDTD/processfields_td.h \ FDTD/processcurrent.h \ - examples/FDTD_examples.h + examples/FDTD_examples.h \ + openems.h diff --git a/openems.cpp b/openems.cpp new file mode 100644 index 0000000..6223b62 --- /dev/null +++ b/openems.cpp @@ -0,0 +1,230 @@ +#include "openems.h" +#include "tools/array_ops.h" +#include "FDTD/operator.h" +#include "FDTD/engine.h" +#include "FDTD/processvoltage.h" +#include "FDTD/processcurrent.h" +#include "FDTD/processfields_td.h" + +//external libs +#include "tinyxml.h" +#include "ContinuousStructure.h" + +openEMS::openEMS() +{ + FDTD_Op=NULL; + FDTD_Eng=NULL; + PA=NULL; +} + +openEMS::~openEMS() +{ + delete FDTD_Op; + FDTD_Op=NULL; + delete FDTD_Eng; + FDTD_Eng=NULL; + delete PA; + PA=NULL; +} + +void openEMS::Reset() +{ + delete FDTD_Op; + FDTD_Op=NULL; + delete FDTD_Eng; + FDTD_Eng=NULL; + if (PA) PA->DeleteAll(); + delete PA; + PA=NULL; +} + +int openEMS::SetupFDTD(const char* file) +{ + if (file==NULL) return -1; + Reset(); + double fmax=0; + int Excit_Type=0; + bool EnableDump = true; + int bounds[6]; + + time_t startTime=time(NULL); + + TiXmlDocument doc(file); + if (!doc.LoadFile()) + { + cerr << "openEMS: Error File-Loading failed!!! File: " << file << endl; + exit(-1); + } + + cerr << "Read openEMS Settings..." << endl; + TiXmlElement* FDTD_Opts = doc.FirstChildElement("openEMS-Parameter"); + if (FDTD_Opts==NULL) + { + cerr << "Can't read openEMS Settings... " << endl; + exit(-1); + } + FDTD_Opts->QueryIntAttribute("NumberOfTimesteps",&NrTS); + + TiXmlElement* Excite = FDTD_Opts->FirstChildElement("Excitation"); + if (Excite==NULL) + { + cerr << "Can't read openEMS Excitation Settings... " << endl; + exit(-2); + } + Excite->QueryIntAttribute("Type",&Excit_Type); + Excite->QueryDoubleAttribute("f0",&fmax); + + TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond"); + if (BC==NULL) + { + cerr << "Can't read openEMS boundary cond Settings... " << endl; + exit(-3); + } + BC->QueryIntAttribute("xmin",&bounds[0]); + BC->QueryIntAttribute("xmax",&bounds[1]); + BC->QueryIntAttribute("ymin",&bounds[2]); + BC->QueryIntAttribute("ymax",&bounds[3]); + BC->QueryIntAttribute("zmin",&bounds[4]); + BC->QueryIntAttribute("zmax",&bounds[5]); + + cerr << "Read Geometry..." << endl; + ContinuousStructure CSX; + string EC(CSX.ReadFromXML(&doc)); + if (EC.empty()==false) + { + cerr << EC << endl; + return(-2); + } + + bool PMC[6]; + for (int n=0;n<6;++n) + PMC[n]=(bounds[n]==1); + + //*************** setup operator ************// + cerr << "Create Operator..." << endl; + FDTD_Op = new Operator(); + if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1); + + FDTD_Op->CalcECOperator(); + + if (Excit_Type==0) + FDTD_Op->CalcGaussianPulsExcitation(fmax/2,fmax/2); + else + { + cerr << "openEMS: Excitation type is unknown" << endl; + exit(-1); + } + + time_t OpDoneTime=time(NULL); + + FDTD_Op->ShowSize(); + + FDTD_Op->ApplyMagneticBC(PMC); + + cerr << "Nyquist number of timesteps: " << FDTD_Op->GetNyquistNum(fmax) << endl; + unsigned int Nyquist = FDTD_Op->GetNyquistNum(fmax); + + cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl; + + //create FDTD engine + FDTD_Eng = new Engine(FDTD_Op); + + time_t currTime = time(NULL); + + //*************** setup processing ************// + PA = new ProcessingArray(); + + double start[3]; + double stop[3]; + vector Probes = CSX.GetPropertyByType(CSProperties::PROBEBOX); + for (size_t i=0;iGetPrimitive(0); + if (prim!=NULL) + { + bool acc; + double* bnd = prim->GetBoundBox(acc,true); + start[0]= bnd[0];start[1]=bnd[2];start[2]=bnd[4]; + stop[0] = bnd[1];stop[1] =bnd[3];stop[2] =bnd[5]; + CSPropProbeBox* pb = Probes.at(i)->ToProbeBox(); + Processing* proc = NULL; + if (pb) + { + if (pb->GetProbeType()==0) + { + ProcessVoltage* procVolt = new ProcessVoltage(FDTD_Op,FDTD_Eng); + procVolt->OpenFile(pb->GetName()); + proc=procVolt; + } + if (pb->GetProbeType()==1) + { + ProcessCurrent* procCurr = new ProcessCurrent(FDTD_Op,FDTD_Eng); + procCurr->OpenFile(pb->GetName()); + proc=procCurr; + } + proc->SetProcessInterval(Nyquist/3); //three times as often as nyquist dictates + proc->DefineStartStopCoord(start,stop); + PA->AddProcessing(proc); + } + else + delete proc; + } + } + + vector DumpProps = CSX.GetPropertyByType(CSProperties::DUMPBOX); + for (size_t i=0;iSetEnable(EnableDump); + ProcTD->SetProcessInterval(Nyquist/2); //twice as often as nyquist dictates + + //only looking for one prim atm + CSPrimitives* prim = DumpProps.at(i)->GetPrimitive(0); + if (prim==NULL) + delete ProcTD; + else + { + bool acc; + double* bnd = prim->GetBoundBox(acc); + start[0]= bnd[0];start[1]=bnd[2];start[2]=bnd[4]; + stop[0] = bnd[1];stop[1] =bnd[3];stop[2] =bnd[5]; + CSPropDumpBox* db = DumpProps.at(i)->ToDumpBox(); + if (db) + { + ProcTD->SetDumpType(db->GetDumpType()); + ProcTD->SetDumpMode(db->GetDumpMode()); + ProcTD->SetFilePattern(db->GetName()); + ProcTD->DefineStartStopCoord(start,stop); + PA->AddProcessing(ProcTD); + } + else + delete ProcTD; + } + } + return 0; +} + +void openEMS::RunFDTD() +{ + time_t currTime = time(NULL); + //*************** simulate ************// + int step=PA->Process(); + if ((step<0) || (step>NrTS)) step=NrTS; + while (FDTD_Eng->GetNumberOfTimesteps()IterateTS(step); + step=PA->Process(); +// cerr << " do " << step << " steps; current: " << eng.GetNumberOfTimesteps() << endl; + if ((step<0) || (step>NrTS - FDTD_Eng->GetNumberOfTimesteps())) step=NrTS - FDTD_Eng->GetNumberOfTimesteps(); + } + + //*************** postproc ************// + time_t prevTime = currTime; + currTime = time(NULL); + + double t_diff = difftime(currTime,prevTime); + + cerr << "Time for " << FDTD_Eng->GetNumberOfTimesteps() << " iterations with " << FDTD_Op->GetNumberCells() << " cells : " << t_diff << " sec" << endl; + cerr << "Speed: " << (double)FDTD_Op->GetNumberCells()*(double)FDTD_Eng->GetNumberOfTimesteps()/t_diff/1e6 << " MCells/s " << endl; +} diff --git a/openems.h b/openems.h new file mode 100644 index 0000000..f33faf6 --- /dev/null +++ b/openems.h @@ -0,0 +1,28 @@ +#ifndef OPENEMS_H +#define OPENEMS_H + +class Operator; +class Engine; +class ProcessingArray; + +class openEMS +{ +public: + openEMS(); + ~openEMS(); + + int SetupFDTD(const char* file); + + void RunFDTD(); + + void Reset(); + +protected: + //! Number of Timesteps + int NrTS; + Operator* FDTD_Op; + Engine* FDTD_Eng; + ProcessingArray* PA; +}; + +#endif // OPENEMS_H