diff --git a/examples/FDTD_examples.cpp b/examples/FDTD_examples.cpp index 994d401..05eff48 100644 --- a/examples/FDTD_examples.cpp +++ b/examples/FDTD_examples.cpp @@ -43,15 +43,46 @@ void BuildDipol(ContinuousStructure &CSX) void BuildPlaneWave(ContinuousStructure &CSX) { -// CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); -//// mat->SetKappa(0.001); -// CSX.AddProperty(mat); -// -// CSPrimBox* matbox = new CSPrimBox(CSX.GetParameterSet(),mat); -// matbox->SetCoord(0,-1000.0);matbox->SetCoord(1,1000.0); -// matbox->SetCoord(2,-1000.0);matbox->SetCoord(3,1000.0); -// matbox->SetCoord(4,-4000.0);matbox->SetCoord(5,4000.0); -// CSX.AddPrimitive(matbox); + double width = 1000; + double hight = 1000; + double length = 4000; + double abs_l = 200; + + 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/2.0-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,width/-2.0);box->SetCoord(1,width/2.0); + box->SetCoord(2,hight/-2.0);box->SetCoord(3,hight/2.0); + box->SetCoord(4,length/2.0-abs_l); box->SetCoord(5,length/2.0); + CSX.AddPrimitive(box); + box = new CSPrimBox(CSX.GetParameterSet(),mat); + box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); + box->SetCoord(2,hight/-2.0);box->SetCoord(3,hight/2.0); + box->SetCoord(4,length/-2.0+abs_l); box->SetCoord(5,length/-2.0); + CSX.AddPrimitive(box); + + CSPropMaterial* mat2 = new CSPropMaterial(CSX.GetParameterSet()); + mat2->SetEpsilon(2); + CSX.AddProperty(mat2); + box = new CSPrimBox(CSX.GetParameterSet(),mat2); + box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); + box->SetCoord(2,hight/-2.0);box->SetCoord(3,hight/2.0); + box->SetCoord(4,length/-4.0);box->SetCoord(5,length/4.0); + CSX.AddPrimitive(box); CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet()); elec->SetExcitation(1,1); @@ -61,28 +92,38 @@ void BuildPlaneWave(ContinuousStructure &CSX) // elec->SetDelay(2.0e-9); CSX.AddProperty(elec); - CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec); - box->SetCoord(0,-500.0);box->SetCoord(1,500.0); - box->SetCoord(2,-500.0);box->SetCoord(3,500.0); - box->SetCoord(4,-4000.0);box->SetCoord(5,-4000.0); + box = new CSPrimBox(CSX.GetParameterSet(),elec); + box->SetCoord(0,width/-2.0);box->SetCoord(1,width/2.0); + box->SetCoord(2,hight/-2.0);box->SetCoord(3,hight/2.0); + box->SetCoord(4,0.0);box->SetCoord(5,0.0); CSX.AddPrimitive(box); - CSPropMetal* metal = new CSPropMetal(CSX.GetParameterSet()); - CSX.AddProperty(metal); - CSPrimCylinder* cyl = new CSPrimCylinder(CSX.GetParameterSet(),metal); - cyl->SetRadius(100); - cyl->SetCoord(0,0.0);cyl->SetCoord(1,0.0); - cyl->SetCoord(2,-250.0);cyl->SetCoord(3,250.0); - cyl->SetCoord(4,-0000.0);cyl->SetCoord(5,-0000.0); - CSX.AddPrimitive(cyl); + //E-field dump + CSPropDumpBox* Edump = new CSPropDumpBox(CSX.GetParameterSet()); + Edump->SetEFieldDump(true); + CSX.AddProperty(Edump); + box = new CSPrimBox(CSX.GetParameterSet(),Edump); + box->SetCoord(0,width/-3.0);box->SetCoord(1,width/3.0); + box->SetCoord(2,0.0);box->SetCoord(3,0.0); + box->SetCoord(4,length/-2.0+abs_l);box->SetCoord(5,length/2.0-abs_l); + CSX.AddPrimitive(box); + +// CSPropMetal* metal = new CSPropMetal(CSX.GetParameterSet()); +// CSX.AddProperty(metal); +// CSPrimCylinder* cyl = new CSPrimCylinder(CSX.GetParameterSet(),metal); +// cyl->SetRadius(100); +// cyl->SetCoord(0,0.0);cyl->SetCoord(1,0.0); +// cyl->SetCoord(2,-250.0);cyl->SetCoord(3,250.0); +// cyl->SetCoord(4,-0000.0);cyl->SetCoord(5,-0000.0); +// CSX.AddPrimitive(cyl); CSRectGrid* grid = CSX.GetGrid(); - for (int n=-500;n<=500;n+=20) + for (int n=width/-2.0;n<=width/2;n+=20) grid->AddDiscLine(0,(double)n); - for (int n=-500;n<=500;n+=20) + for (int n=hight/-2.0;n<=hight/2.0;n+=20) grid->AddDiscLine(1,(double)n); - for (int n=-4000;n<=4000;n+=20) + for (int n=length/-2.0;n<=length/2.0;n+=20) grid->AddDiscLine(2,(double)n); grid->SetDeltaUnit(1e-3); @@ -100,7 +141,7 @@ void BuildMSL(ContinuousStructure &CSX) //substrate.... CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); // mat->SetEpsilon(3.6); - double finalKappa = 3.0/pow(abs_l,4); + double finalKappa = 0.3/pow(abs_l,4); mat->SetKappa(finalKappa); std::ostringstream fct; fct << "pow(abs(z)-" << length/2.0-abs_l << ",4)"; diff --git a/main.cpp b/main.cpp index 681dfe7..737cc87 100644 --- a/main.cpp +++ b/main.cpp @@ -23,8 +23,8 @@ int main(int argc, char *argv[]) cerr << "Create Geometry..." << endl; ContinuousStructure CSX; - //BuildPlaneWave(CSX); - BuildMSL(CSX); + BuildPlaneWave(CSX); +// BuildMSL(CSX); //*************** setup operator ************// cerr << "Create Operator..." << endl; @@ -39,7 +39,7 @@ int main(int argc, char *argv[]) time_t OpDoneTime=time(NULL); cop.ShowSize(); - bool bnd[] = {1,1,0,1,0,1}; + bool bnd[] = {1,1,0,0,0,0}; cop.ApplyMagneticBC(bnd); cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl;