diff --git a/examples/Coax_Cart.xml b/examples/Coax_Cart.xml new file mode 100644 index 0000000..33ecdc4 --- /dev/null +++ b/examples/Coax_Cart.xml @@ -0,0 +1,100 @@ + + + + + + + + -230,-225,-220,-215,-210,-205,-200,-195,-190,-185,-180,-175,-170,-165,-160,-155,-150,-145,-140,-135,-130,-125,-120,-115,-110,-105,-100,-95,-90,-85,-80,-75,-70,-65,-60,-55,-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,105,110,115,120,125,130,135,140,145,150,155,160,165,170,175,180,185,190,195,200,205,210,215,220,225,230 + -230,-225,-220,-215,-210,-205,-200,-195,-190,-185,-180,-175,-170,-165,-160,-155,-150,-145,-140,-135,-130,-125,-120,-115,-110,-105,-100,-95,-90,-85,-80,-75,-70,-65,-60,-55,-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,105,110,115,120,125,130,135,140,145,150,155,160,165,170,175,180,185,190,195,200,205,210,215,220,225,230 + -500,-490,-480,-470,-460,-450,-440,-430,-420,-410,-400,-390,-380,-370,-360,-350,-340,-330,-320,-310,-300,-290,-280,-270,-260,-250,-240,-230,-220,-210,-200,-190,-180,-170,-160,-150,-140,-130,-120,-110,-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,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,510,520,530,540,550,560,570,580,590,600,610,620,630,640,650,660,670,680,690,700,710,720,730,740,750,760,770,780,790,800,810,820,830,840,850,860,870,880,890,900,910,920,930,940,950,960,970,980,990,1000,1010,1020,1030,1040,1050,1060,1070,1080,1090,1100,1110,1120,1130,1140,1150,1160,1170,1180,1190,1200,1210,1220,1230,1240,1250,1260,1270,1280,1290,1300,1310,1320,1330,1340,1350,1360,1370,1380,1390,1400,1410,1420,1430,1440,1450,1460,1470,1480,1490,1500 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/Dipol.xml b/examples/Dipol.xml index 24cdd7c..08c12b7 100644 --- a/examples/Dipol.xml +++ b/examples/Dipol.xml @@ -1,6 +1,6 @@ - + @@ -25,7 +25,7 @@ - + diff --git a/examples/FDTD_examples.cpp b/examples/FDTD_examples.cpp index 12f262a..1bcfd18 100644 --- a/examples/FDTD_examples.cpp +++ b/examples/FDTD_examples.cpp @@ -4,7 +4,8 @@ void BuildDipol(const char* filename) { int maxIter = 1000; - double fmax=1e9; + double f0=0.5e9; + double fc=0.5e9; int Excit_Type=0; int bounds[] = {1,1,0,0,0,0}; @@ -55,7 +56,8 @@ void BuildDipol(const char* filename) TiXmlElement Excite("Excitation"); Excite.SetAttribute("Type",Excit_Type); - Excite.SetAttribute("f0",fmax); + Excite.SetAttribute("f0",f0); + Excite.SetAttribute("fc",fc); FDTD_Opts.InsertEndChild(Excite); TiXmlElement BC("BoundaryCond"); @@ -81,7 +83,8 @@ void BuildDipol(const char* filename) void BuildPlaneWave(const char* filename) { int maxIter = 1000; - double fmax=1e9; + double f0=0.5e9; + double fc=0.5e9; int Excit_Type=0; int bounds[] = {1,1,0,0,0,0}; @@ -194,7 +197,8 @@ void BuildPlaneWave(const char* filename) TiXmlElement Excite("Excitation"); Excite.SetAttribute("Type",Excit_Type); - Excite.SetAttribute("f0",fmax); + Excite.SetAttribute("f0",f0); + Excite.SetAttribute("fc",fc); FDTD_Opts.InsertEndChild(Excite); TiXmlElement BC("BoundaryCond"); @@ -220,7 +224,8 @@ void BuildPlaneWave(const char* filename) void BuildMSL(const char* filename) { int maxIter = 1000; - double fmax=1e9; + double f0=0.5e9; + double fc=0.5e9; int Excit_Type=0; int bounds[] = {1,1,0,0,0,0}; @@ -365,7 +370,8 @@ void BuildMSL(const char* filename) TiXmlElement Excite("Excitation"); Excite.SetAttribute("Type",Excit_Type); - Excite.SetAttribute("f0",fmax); + Excite.SetAttribute("f0",f0); + Excite.SetAttribute("fc",fc); FDTD_Opts.InsertEndChild(Excite); TiXmlElement BC("BoundaryCond"); @@ -388,4 +394,173 @@ void BuildMSL(const char* filename) doc.SaveFile(); } +void BuildCoaxial_Cartesian(const char* filename) +{ + int maxIter = 1000; + double f0=1e9; + double fc=1e9; + int Excit_Type=0; + int bounds[] = {0,0,0,0,0,0}; + + cerr << "Create Geometry..." << endl; + ContinuousStructure CSX; + + double rad[] = {100, 230}; + double length[] = {-500,1500}; + double abs_l = 200; + + double delta[] = {5,5,10}; + + CSPrimBox* box = NULL; + //fake pml.... + CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); +// mat->SetEpsilon(3.6); + double finalKappa = 0.3/pow(abs_l,4); + mat->SetKappa(finalKappa); + std::ostringstream fct; + fct << "pow(abs(z)-" << length[1]-abs_l << ",4)"; + mat->SetKappaWeightFunction(fct.str(),0); + mat->SetKappaWeightFunction(fct.str(),1); + mat->SetKappaWeightFunction(fct.str(),2); + mat->SetSigma(finalKappa*__MUE0__/__EPS0__); + mat->SetSigmaWeightFunction(fct.str(),0); + mat->SetSigmaWeightFunction(fct.str(),1); + mat->SetSigmaWeightFunction(fct.str(),2); + CSX.AddProperty(mat); + + box = new CSPrimBox(CSX.GetParameterSet(),mat); + box->SetCoord(0,-1.0*rad[1]);box->SetCoord(1,1.0*rad[1]); + box->SetCoord(2,-1.0*rad[1]);box->SetCoord(3,1.0*rad[1]); + box->SetCoord(4,length[1]-abs_l); box->SetCoord(5,length[1]); + box->SetPriority(10); + CSX.AddPrimitive(box); + + CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet()); + elec->SetExcitation(1.0,0); + elec->SetExcitation(1.0,1); + elec->SetWeightFunction("x/pow(rho,2)",0); + elec->SetWeightFunction("y/pow(rho,2)",1); + elec->SetExcitType(0); +// elec->SetActiveDir(0,0);//disable x +// elec->SetActiveDir(0,2);//disable z +// elec->SetDelay(2.0e-9); + CSX.AddProperty(elec); + +// double coords[] = {-100,-100,0}; +// cerr << elec->GetWeightedExcitation(0,coords) << endl; +// cerr << elec->GetWeightedExcitation(1,coords) << endl; +// exit(0); + + box = new CSPrimBox(CSX.GetParameterSet(),elec); + box->SetCoord(0,-1.0*rad[1]);box->SetCoord(1,1.0*rad[1]); + box->SetCoord(2,-1.0*rad[1]);box->SetCoord(3,1.0*rad[1]); + box->SetCoord(4,length[0]);box->SetCoord(5,length[0]); + box->SetPriority(5); + CSX.AddPrimitive(box); + + //E-field dump + CSPropDumpBox* Edump = new CSPropDumpBox(CSX.GetParameterSet()); + Edump->SetDumpType(0); + Edump->SetDumpMode(2); + Edump->SetName("tmp/Et_"); + CSX.AddProperty(Edump); + box = new CSPrimBox(CSX.GetParameterSet(),Edump); + box->SetCoord(0,-1*rad[1]);box->SetCoord(1,rad[1]); + box->SetCoord(2,-0*rad[1]);box->SetCoord(3,0*rad[1]); + box->SetCoord(4,length[0]);box->SetCoord(5,length[1]); + 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,rad[0]);box->SetCoord(1,rad[1]); + box->SetCoord(2,0.0);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,-1.5*rad[0]);box->SetCoord(1,1.5*rad[0]); + box->SetCoord(2,-1.5*rad[0]);box->SetCoord(3,1.5*rad[0]); + box->SetCoord(4,0.0);box->SetCoord(5,0.0); + CSX.AddPrimitive(box); + + CSPropMaterial* metal = new CSPropMaterial(CSX.GetParameterSet()); + metal->SetKappa(56e9); +// CSPropMetal* metal = new CSPropMetal(CSX.GetParameterSet()); + CSX.AddProperty(metal); + CSPrimCylinder* cyl = new CSPrimCylinder(CSX.GetParameterSet(),metal); + cyl->SetRadius(rad[0]); + cyl->SetCoord(0,0.0);cyl->SetCoord(1,0.0); + cyl->SetCoord(2,0.0);cyl->SetCoord(3,0.0); + cyl->SetCoord(4,length[0]);cyl->SetCoord(5,length[1]); + cyl->SetPriority(100); + CSX.AddPrimitive(cyl); + box = new CSPrimBox(CSX.GetParameterSet(),metal); + box->SetCoord(0,-1.0*rad[1]);box->SetCoord(1,1.0*rad[1]); + box->SetCoord(2,-1.0*rad[1]);box->SetCoord(3,1.0*rad[1]); + box->SetCoord(4,length[0]);box->SetCoord(5,length[1]); + box->SetPriority(1); + CSX.AddPrimitive(box); + + CSPropMaterial* air = new CSPropMaterial(CSX.GetParameterSet()); + CSX.AddProperty(air); + cyl = new CSPrimCylinder(CSX.GetParameterSet(),air); + cyl->SetRadius(rad[1]); + cyl->SetCoord(0,0.0);cyl->SetCoord(1,0.0); + cyl->SetCoord(2,0.0);cyl->SetCoord(3,0.0); + cyl->SetCoord(4,length[0]);cyl->SetCoord(5,length[1]); + cyl->SetPriority(9); + CSX.AddPrimitive(cyl); + + CSRectGrid* grid = CSX.GetGrid(); + + for (int n=-1.0*rad[1];n<=rad[1];n+=delta[0]) + grid->AddDiscLine(0,(double)n); + for (int n=-1.0*rad[1];n<=rad[1];n+=delta[1]) + grid->AddDiscLine(1,(double)n); + for (int n=length[0];n<=length[1];n+=delta[2]) + grid->AddDiscLine(2,(double)n); + + grid->SetDeltaUnit(1e-3); + + //*************** 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",f0); + Excite.SetAttribute("fc",fc); + 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 69b0cc6..4c1129a 100644 --- a/examples/FDTD_examples.h +++ b/examples/FDTD_examples.h @@ -10,4 +10,7 @@ void BuildPlaneWave(const char* filename); void BuildMSL(const char* filename); +void BuildCoaxial_Cartesian(const char* filename); + + #endif // FDTD_EXAMPLES_H diff --git a/examples/MSL.xml b/examples/MSL.xml index fea7642..4d82b30 100644 --- a/examples/MSL.xml +++ b/examples/MSL.xml @@ -1,6 +1,6 @@ - + @@ -35,7 +35,7 @@ - + diff --git a/examples/PlaneWave.xml b/examples/PlaneWave.xml index 854d8e1..6e2273b 100644 --- a/examples/PlaneWave.xml +++ b/examples/PlaneWave.xml @@ -1,6 +1,6 @@ - + @@ -35,7 +35,7 @@ - + diff --git a/main.cpp b/main.cpp index 2b39372..e58af1c 100644 --- a/main.cpp +++ b/main.cpp @@ -32,7 +32,10 @@ int main(int argc, char *argv[]) const char* fileMSL="examples/MSL.xml"; BuildMSL(fileMSL); - const char* file=fileMSL; + const char* fileCoax="examples/Coax_Cart.xml"; + BuildCoaxial_Cartesian(fileCoax); + + const char* file=fileCoax; // cerr << CSX.ReadFromXML("examples/PlaneWave.xml") << endl; #endif diff --git a/openems.cpp b/openems.cpp index 6223b62..b91841f 100644 --- a/openems.cpp +++ b/openems.cpp @@ -42,7 +42,8 @@ int openEMS::SetupFDTD(const char* file) { if (file==NULL) return -1; Reset(); - double fmax=0; + double f0=0; + double fc=0; int Excit_Type=0; bool EnableDump = true; int bounds[6]; @@ -72,7 +73,11 @@ int openEMS::SetupFDTD(const char* file) exit(-2); } Excite->QueryIntAttribute("Type",&Excit_Type); - Excite->QueryDoubleAttribute("f0",&fmax); + if (Excit_Type==0) + { + Excite->QueryDoubleAttribute("f0",&f0); + Excite->QueryDoubleAttribute("fc",&fc); + } TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond"); if (BC==NULL) @@ -108,7 +113,7 @@ int openEMS::SetupFDTD(const char* file) FDTD_Op->CalcECOperator(); if (Excit_Type==0) - FDTD_Op->CalcGaussianPulsExcitation(fmax/2,fmax/2); + FDTD_Op->CalcGaussianPulsExcitation(f0,fc); else { cerr << "openEMS: Excitation type is unknown" << endl; @@ -121,8 +126,8 @@ int openEMS::SetupFDTD(const char* file) FDTD_Op->ApplyMagneticBC(PMC); - cerr << "Nyquist number of timesteps: " << FDTD_Op->GetNyquistNum(fmax) << endl; - unsigned int Nyquist = FDTD_Op->GetNyquistNum(fmax); + cerr << "Nyquist number of timesteps: " << FDTD_Op->GetNyquistNum(f0+fc) << endl; + unsigned int Nyquist = FDTD_Op->GetNyquistNum(f0+fc); cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl;