Merge branch 'grating3D' into 'master'

Grating3d

See merge request !4

DiffractionGratings/grating3D.geo

+// test_case = "pyramid";
+// test_case = "hole";
+// test_case = "torus";
+// test_case = "halfellipsoid";
+// test_case = "checker";
+// test_case = "bisin";
+Include StrCat["grating3D_data_",test_case,".geo"];
+
+SetFactory("OpenCASCADE");
+
+Macro SetPBCs
+  e = 10*nm;
+  masterX() = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e,-period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+  slaveX()  = Surface In BoundingBox{ period_x/2-e,-period_y/2-e, PML_bot_hh-e, period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+  masterY() = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e, period_x/2+e,-period_y/2+e, PML_top_hh+PML_top+e};
+  slaveY()  = Surface In BoundingBox{-period_x/2-e, period_y/2-e, PML_bot_hh-e, period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+  If (tag_geom==6) // bi-sin : BoundingBox does not catch BSpline Surfaces?
+    masterX()+={48,52};
+    slaveX()+={50,54};
+    masterY()+={51,55};
+    slaveY()+={49,53};
+  EndIf
+  Periodic Surface{slaveX()} = {masterX()} Translate{period_x,        0,        0};
+  Periodic Surface{slaveY()} = {masterY()} Translate{       0, period_y,        0};
+Return
+
+If (tag_geom==6)
+  // Bi-sinusoidal profile by Spline surface filling
+  // not the most general way to implement a freesurface height=f(x,y) and relies on multiple Coherence calls
+  nsin = 20;
+  For i In {0:nsin}
+      xx1~{i} = -period_x/2+i*period_x/2/(nsin-1);
+      xx2~{i} =             i*period_x/2/(nsin-1);
+  EndFor
+  pp=newp;//pp=pp-1;
+  For i In {0:nsin-1}
+    yy=-period_y/2;
+    Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+    Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+    yy=0;
+    Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+    Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+    yy=period_y/2;
+    Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+    Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+    yy=-period_y/2;
+    Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+    Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+    yy=0;
+    Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+    Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+    yy=period_y/2;
+    Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+    Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+  EndFor
+  For k In {1:11}
+    BSpline(newl) = {pp+k:pp+12*nsin:12};
+  EndFor
+  BSpline(newl) = {pp  :pp+12*(nsin-1):12};
+  Coherence;
+  Box(1) = {-period_x/2,-period_y/2,     hh_L_3,period_x,period_y, thick_L_3};
+  Delete{Volume{1}; Surface{1,2,3,4} ; Line{13,15,17,19};}
+  Curve Loop(7) = {21, 20, -23, -16};
+  Curve Loop(8) = {1, 2, -3, -12};
+  Surface(7) = {8};
+  Curve Loop(10) = {6, 5, -8, -3};
+  Surface(8) = {10};
+  Curve Loop(12) = {9, 10, -11, -8};
+  Surface(9) = {12};
+  Curve Loop(14) = {2, 9, -4, -7};
+  Surface(10) = {14};
+  Line(25) = {239, 229};
+  Line(26) = {229, 238};
+  Line(27) = {243, 235};
+  Line(28) = {235, 242};
+  Line(29) = {237, 244};
+  Line(30) = {233, 240};
+  Line(31) = {237, 245};
+  Line(32) = {241, 233};
+  Curve Loop(16) = {25, 12, 6, -27, -21};
+  Surface(11) = {16};
+  Curve Loop(18) = {20, -31, -11, -5, -27};
+  Surface(12) = {18};
+  Curve Loop(20) = {29, -18, -28, 5, 11};
+  Surface(13) = {20};
+  Curve Loop(22) = {10, 29, -24, -30, 4};
+  Surface(14) = {22};
+  Curve Loop(24) = {30, -14, -26, 1, 7};
+  Surface(15) = {24};
+  Curve Loop(26) = {25, 1, 7, -32, -16};
+  Surface(16) = {26};
+  Curve Loop(28) = {26, 22, -28, -6, -12};
+  Surface(17) = {28};
+  Curve Loop(30) = {31, -23, 32, 4, 10};
+  Surface(18) = {30};
+  Surface Loop(3) = {5, 16, 18, 12, 11, 9, 8, 7, 10};
+  Volume(9) = {3};    
+  Surface Loop(2) = {6, 15, 14, 10, 7, 8, 9, 13, 17};
+  Volume(10) = {2};
+EndIf
+
+L1_n     = Sqrt(eps_re_L_1);
+L2_n     = Sqrt(eps_re_L_2);
+L3_n     = Sqrt(eps_re_L_3);
+L4_n     = Sqrt(eps_re_L_4);
+L5_n     = Sqrt(eps_re_L_5);
+L6_n     = Sqrt(eps_re_L_6);
+Scat_n   = 2 ;
+
+PML_lc     = lambda_m/(paramaille*.5);
+L1_lc      = lambda_m/(paramaille*L1_n);
+L2_lc      = lambda_m/(paramaille*L2_n);
+L3_lc      = lambda_m/(paramaille*L3_n);
+L4_lc      = lambda_m/(paramaille*L4_n);
+L5_lc      = lambda_m/(paramaille*L5_n);
+L6_lc      = lambda_m/(paramaille*L6_n);
+
+Box(1) = {-period_x/2,-period_y/2, PML_bot_hh,period_x,period_y, PML_bot  };
+Box(2) = {-period_x/2,-period_y/2,     hh_L_6,period_x,period_y, thick_L_6};
+Box(3) = {-period_x/2,-period_y/2,     hh_L_5,period_x,period_y, thick_L_5};
+Box(4) = {-period_x/2,-period_y/2,     hh_L_4,period_x,period_y, thick_L_4};
+If (tag_geom!=6)
+  Box(5) = {-period_x/2,-period_y/2,     hh_L_3,period_x,period_y, thick_L_3};
+EndIf
+Box(6) = {-period_x/2,-period_y/2,     hh_L_2,period_x,period_y, thick_L_2};
+Box(7) = {-period_x/2,-period_y/2,     hh_L_1,period_x,period_y, thick_L_1};
+Box(8) = {-period_x/2,-period_y/2, PML_top_hh,period_x,period_y, PML_top  };
+
+If (tag_geom==1)
+  p=newp;
+  Point(p) = {0,0,hh_L_3+rz};
+  Line(97) = {25, 65};
+  Line(98) = {65, 27};
+  Line(99) = {65, 29};
+  Line(100) = {31, 65};
+  Curve Loop(49) = {38, -98, -97};
+  Plane Surface(49) = {49};
+  Curve Loop(50) = {98, 48, 100};
+  Plane Surface(50) = {50};
+  Curve Loop(51) = {99, 42, 100};
+  Plane Surface(51) = {51};
+  Curve Loop(52) = {99, -46, 97};
+  Plane Surface(52) = {52};
+  Surface Loop(9) = {52, 51, 50, 49, 24};
+  Volume(9) = {9};
+EndIf
+
+If (tag_geom==2)
+  Cylinder(9) = {0,0,hh_L_3,0,0,thick_L_3,rx};
+EndIf
+
+If (tag_geom==3)
+  Torus(9) = {0,0,hh_L_3+ry-1,rx,ry};
+  Dilate { { 0,0,hh_L_3 }, { 1, 1, rz/ry } } { Volume{9}; }
+  BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+EndIf
+
+If (tag_geom==4)
+  Sphere(9) = {0,0,hh_L_3,rx};
+  Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, rz/rx } } { Volume{9}; }
+  BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+EndIf
+
+If (tag_geom==5)
+  Box(9) = {-rx/2,-ry/2, hh_L_2-rz, rx, ry, rz};
+  Rotate {{0, 0, 1}, {0,0,0}, Pi/4} {Volume{9};}
+EndIf
+
+
+
+
+Coherence;
+// // // // Physical Surface("ENDPML",770) = {69,37};
+
+Call SetPBCs;
+
+Physical Volume("PMLBOT"        ,1) = {1};
+Physical Volume("LAYER_L6_SUBS" ,2) = {2};
+Physical Volume("LAYER_L5_LOW"  ,3) = {3};
+Physical Volume("LAYER_L4_GUI"  ,4) = {4};
+Physical Volume("LAYER_L3_EMB"  ,5) = {10};
+Physical Volume("LAYER_L2_TOP"  ,6) = {6};
+Physical Volume("LAYER_L1_SUPER",7) = {7};
+Physical Volume("PMLTOP"        ,8) = {8};
+Physical Volume("SCAT"          ,9) = {9};
+
+Physical Surface("BXM" ,101 ) = masterX();
+Physical Surface("BXP" ,102 ) = slaveX();
+Physical Surface("BYM" ,201 ) = masterY();
+Physical Surface("BYP" ,202 ) = slaveY();
+Physical Surface("STOP",301 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_top_hh-e,period_x/2+e, period_y/2+e, PML_top_hh+e};
+Physical Surface("SBOT",302 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, hh_L_6-e,period_x/2+e, period_y/2+e, hh_L_6+e};
+Physical Surface("SPMLTOP",401 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_top_hh+PML_top-e,period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+Physical Surface("SPMLBOT",402 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e,period_x/2+e, period_y/2+e, PML_bot_hh+e};
+
+pts_PMLBOT()   = PointsOf{ Physical Volume{1}; };
+pts_LAYER_L6() = PointsOf{ Physical Volume{2}; };
+pts_LAYER_L5() = PointsOf{ Physical Volume{3}; };
+pts_LAYER_L4() = PointsOf{ Physical Volume{4}; };
+pts_LAYER_L3() = PointsOf{ Physical Volume{5}; };
+pts_LAYER_L2() = PointsOf{ Physical Volume{6}; };
+pts_LAYER_L1() = PointsOf{ Physical Volume{7}; };
+pts_PMLTOP()   = PointsOf{ Physical Volume{8}; };
+pts_ROD()      = PointsOf{ Physical Volume{9}; };
+
+Characteristic Length{:}              = PML_lc;
+Characteristic Length{pts_LAYER_L1()} = L1_lc/refine_mesh_L1;
+Characteristic Length{pts_LAYER_L2()} = L2_lc/refine_mesh_L2;
+Characteristic Length{pts_LAYER_L4()} = L4_lc/refine_mesh_L4;
+Characteristic Length{pts_LAYER_L5()} = L5_lc/refine_mesh_L5;
+Characteristic Length{pts_LAYER_L6()} = L6_lc/refine_mesh_L6;
+Characteristic Length{pts_LAYER_L3()} = L3_lc/refine_mesh_L3;
+Characteristic Length{pts_ROD()}      = scat_lc;
+If (tag_geom==3) // Split torus weird otherwise
+  Mesh.Algorithm = 6;
+EndIf
+Mesh.Optimize=1;
+Mesh.SurfaceEdges = 0;
+Mesh.VolumeEdges = 0;

DiffractionGratings/grating3D_data_bisin.geo

+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+    lambda0       = {830   , Name StrCat[pp1,"/1lambda0 [nm]"]},
+    thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
+    phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
+    psideg        = {0     , Name StrCat[pp1,"/4psi0 [deg]"]},
+    period_x      = {1000  , Name StrCat[pp2,"/1X period [nm]"]},
+    period_y      = {1000  , Name StrCat[pp2,"/2Y period [nm]"]},
+    thick_L_1     = {100   , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+    thick_L_2     = {100   , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    thick_L_3     = {500, Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+    thick_L_4     = {100   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+    thick_L_5     = {100   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+    thick_L_6     = {200   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+
+    tag_geom      = {  6   , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+    rx            = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
+    ry            = {1.25*lambda0, Name StrCat[pp3,"/2ry"]},
+    rz            = { 200    , Name StrCat[pp3,"/3rz"]},
+    flag_mat_scat = { 0    , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_Scat   = { 4    , Name StrCat[pp3,"/7eps_re_Scat"]},
+    eps_im_Scat   = { 0    , Name StrCat[pp3,"/8eps_im_Scat"]},
+
+    flag_mat_1    = { 0    , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_2    = { 0    , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_3    = { 0    , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_4    = { 0    , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_5    = { 0    , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_6    = { 0    , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_L_1    = {1     , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+    eps_im_L_1    = {0     , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+    eps_re_L_2    = {1     , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+    eps_im_L_2    = {0     , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+    eps_re_L_3    = {1     , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+    eps_im_L_3    = {0     , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+    eps_re_L_4    = {4     , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+    eps_im_L_4    = {0     , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+    eps_re_L_5    = {4     , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+    eps_im_L_5    = {0     , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+    eps_re_L_6    = {4     , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+    eps_im_L_6    = {0     , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+    
+    paramaille    = {7     , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+    scat_lc       = {lambda0/(2*paramaille)    , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+    PML_top       = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]},
+    PML_bot       = {lambda0*1.2, Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+    Nmax          = {2      , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+    refine_mesh_L1= {1      , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+    refine_mesh_L2= {1      , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+    refine_mesh_L3= {1      , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    refine_mesh_L4= {1      , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+    refine_mesh_L5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+    refine_mesh_L6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    
+    InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+    Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+    FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+    FlagOutEtotFull    = { 0    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutEscaFull    = { 0    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutPoyFull     = { 0    , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+    hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0   = phideg*Pi/180;
+psi0   = psideg*Pi/180;
+
+DomainZsizeSca  = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot  = PML_top_hh-hh_L_6;
+npts_interpX    = period_x/InterpSampling;
+npts_interpY    = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;

DiffractionGratings/grating3D_data_checker.geo

+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+    lambda0       = {500   , Name StrCat[pp1,"/1lambda0 [nm]"]},
+    thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
+    phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
+    psideg        = {45    , Name StrCat[pp1,"/4psi0 [deg]"]},
+    period_x      = {1.25*Sqrt[2]*lambda0, Name StrCat[pp2,"/1X period [nm]"]},
+    period_y      = {1.25*Sqrt[2]*lambda0, Name StrCat[pp2,"/2Y period [nm]"]},
+    thick_L_1     = {50    , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+    thick_L_2     = {50    , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    thick_L_3     = {1.2*lambda0, Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+    thick_L_4     = {50    , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+    thick_L_5     = {50    , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+    thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+
+    tag_geom      = {  5   , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+    rx            = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
+    ry            = {1.25*lambda0, Name StrCat[pp3,"/2ry"]},
+    rz            = {lambda0     , Name StrCat[pp3,"/3rz"]},
+    flag_mat_scat = { 0    , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_Scat   = { 2.25 , Name StrCat[pp3,"/7eps_re_Scat"]},
+    eps_im_Scat   = { 0    , Name StrCat[pp3,"/8eps_im_Scat"]},
+
+    flag_mat_1    = { 0    , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_2    = { 0    , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_3    = { 0    , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_4    = { 0    , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_5    = { 0    , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_6    = { 0    , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_L_1    = {2.25  , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+    eps_im_L_1    = {0     , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+    eps_re_L_2    = {2.25  , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+    eps_im_L_2    = {0     , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+    eps_re_L_3    = {1     , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+    eps_im_L_3    = {0     , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+    eps_re_L_4    = {1     , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+    eps_im_L_4    = {0     , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+    eps_re_L_5    = {1     , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+    eps_im_L_5    = {0     , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+    eps_re_L_6    = {1     , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+    eps_im_L_6    = {0     , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+    
+    paramaille    = {6     , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+    scat_lc       = {lambda0/(1.5*paramaille)    , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+    PML_top       = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]},
+    PML_bot       = {lambda0, Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+    Nmax          = {2      , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+    refine_mesh_L1= {1.5    , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+    refine_mesh_L2= {1.5    , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+    refine_mesh_L3= {1      , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    refine_mesh_L4= {1      , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+    refine_mesh_L5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+    refine_mesh_L6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    
+    InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+    Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+    FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+    FlagOutEtotFull    = { 0    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutEscaFull    = { 0    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutPoyFull     = { 0    , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+    hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0   = phideg*Pi/180;
+psi0   = psideg*Pi/180;
+
+DomainZsizeSca  = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot  = PML_top_hh-hh_L_6;
+npts_interpX    = period_x/InterpSampling;
+npts_interpY    = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;

DiffractionGratings/grating3D_data_halfellipsoid.geo

+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+    lambda0       = {495  , Name StrCat[pp1,"/1lambda0 [nm]"]},
+    thetadeg      = {40   , Name StrCat[pp1,"/2theta0 [deg]"]},
+    phideg        = {36   , Name StrCat[pp1,"/3phi0 [deg]"]},
+    psideg        = {72   , Name StrCat[pp1,"/4psi0 [deg]"]},
+    period_x      = {250  , Name StrCat[pp2,"/1X period [nm]"]},
+    period_y      = {250  , Name StrCat[pp2,"/2Y period [nm]"]},
+    thick_L_1     = {50   , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+    thick_L_2     = {50   , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    thick_L_3     = {100  , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+    thick_L_4     = {50   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+    thick_L_5     = {50   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+    thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+
+    tag_geom      = {  4   , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+    rx            = {107  , Name StrCat[pp3,"/1rx"]},
+    ry            = {47   , Name StrCat[pp3,"/2ry"]},
+    rz            = {40   , Name StrCat[pp3,"/3rz"]},
+    flag_mat_scat = { 0   , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_Scat   = {-2.23, Name StrCat[pp3,"/7eps_re_Scat"]},
+    eps_im_Scat   = { 3.89, Name StrCat[pp3,"/8eps_im_Scat"]},
+
+    flag_mat_1    = { 0    , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_2    = { 0    , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_3    = { 0    , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_4    = { 0    , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_5    = { 0    , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_6    = { 0    , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_L_1    = {1    , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+    eps_im_L_1    = {0    , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+    eps_re_L_2    = {1    , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+    eps_im_L_2    = {0    , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+    eps_re_L_3    = {1    , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+    eps_im_L_3    = {0    , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+    eps_re_L_4    = {1    , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+    eps_im_L_4    = {0    , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+    eps_re_L_5    = {1    , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+    eps_im_L_5    = {0    , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+    eps_re_L_6    = {4    , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+    eps_im_L_6    = {0    , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+    
+    paramaille    = {8    , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+    scat_lc       = {10   , Name StrCat[pp5,"/2metal mesh size [nm]"]},
+    PML_top       = {lambda0 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+    PML_bot       = {lambda0 , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+    Nmax          = {1     , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+    refine_mesh_L1= {1     , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+    refine_mesh_L2= {1     , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+    refine_mesh_L3= {1     , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    refine_mesh_L4= {1     , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+    refine_mesh_L5= {1     , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+    refine_mesh_L6= {1     , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+        
+    InterpSampling     = { 10   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+    Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+    FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+    FlagOutEtotFull    = { 0    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutEscaFull    = { 0    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutPoyFull     = { 0    , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+]; 
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+    hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0   = phideg*Pi/180;
+psi0   = psideg*Pi/180;
+
+DomainZsizeSca  = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot  = PML_top_hh-hh_L_6;
+npts_interpX    = period_x/InterpSampling;
+npts_interpY    = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;

DiffractionGratings/grating3D_data_hole.geo

+nm = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+    lambda0       = {500   , Name StrCat[pp1,"/1lambda0 [nm]"]},
+    thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
+    phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
+    psideg        = {45    , Name StrCat[pp1,"/4psi0 [deg]"]},
+    period_x      = {999.9 , Name StrCat[pp2,"/1X period [nm]"]},
+    period_y      = {999.9 , Name StrCat[pp2,"/2Y period [nm]"]},
+    thick_L_1     = {100   , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+    thick_L_2     = {100   , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    thick_L_3     = {50    , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+    thick_L_4     = {100   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+    thick_L_5     = {100   , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+    thick_L_6     = {100   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+
+    tag_geom      = {  2   , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+    rx            = {250   , Name StrCat[pp3,"/1rx"]},
+    ry            = {47    , Name StrCat[pp3,"/2ry"]},
+    rz            = {500   , Name StrCat[pp3,"/3rz"]},
+    flag_mat_scat = { 0    , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_Scat   = {1     , Name StrCat[pp3,"/7eps_re_Scat"]},
+    eps_im_Scat   = {0     , Name StrCat[pp3,"/8eps_im_Scat"]},
+
+    flag_mat_1    = { 0    , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_2    = { 0    , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_3    = { 0    , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_4    = { 0    , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_5    = { 0    , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_6    = { 0    , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Set Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_L_1    = {1     , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+    eps_im_L_1    = {0     , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+    eps_re_L_2    = {1     , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+    eps_im_L_2    = {0     , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+    eps_re_L_3    = {0.8125, Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+    eps_im_L_3    = {5.25  , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+    eps_re_L_4    = {2.25  , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+    eps_im_L_4    = {0     , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+    eps_re_L_5    = {2.25  , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+    eps_im_L_5    = {0     , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+    eps_re_L_6    = {2.25  , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+    eps_im_L_6    = {0     , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+    
+    paramaille    = {5     , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+    scat_lc       = {lambda0/(3*paramaille)    , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+    PML_top       = {lambda0, Name StrCat[pp5,"/4PML top thickness [nm]"]},
+    PML_bot       = {lambda0, Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+    Nmax          = {2      , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+    refine_mesh_L1= {1      , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+    refine_mesh_L2= {1      , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+    refine_mesh_L3= {5      , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    refine_mesh_L4= {1      , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+    refine_mesh_L5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+    refine_mesh_L6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    
+    InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+    Flag_interp_cubic  = { 0    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+    FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+    FlagOutEtotFull    = { 0    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutEscaFull    = { 0    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutPoyFull     = { 0    , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+    hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0   = phideg*Pi/180;
+psi0   = psideg*Pi/180;
+
+DomainZsizeSca  = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot  = PML_top_hh-hh_L_6;
+npts_interpX    = period_x/InterpSampling;
+npts_interpY    = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;

DiffractionGratings/grating3D_data_pyramid.geo

+nm  = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+    lambda0       = {1533  , Name StrCat[pp1,"/1lambda0 [nm]"]},
+    thetadeg      = {30    , Name StrCat[pp1,"/2theta0 [deg]"]},
+    phideg        = {45    , Name StrCat[pp1,"/3phi0 [deg]"]},
+    psideg        = { 0    , Name StrCat[pp1,"/4psi0 [deg]"]},
+    period_x      = {1500  , Name StrCat[pp2,"/1X period [nm]"]},
+    period_y      = {1000  , Name StrCat[pp2,"/2Y period [nm]"]},
+    thick_L_1     = {50    , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+    thick_L_2     = {50    , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    thick_L_3     = {300   , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+    thick_L_4     = {50    , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+    thick_L_5     = {50    , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+    thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+    
+    tag_geom      = {  1   , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+    rx            = {107   , Name StrCat[pp3,"/1rx"]},
+    ry            = {250   , Name StrCat[pp3,"/2ry"]},
+    rz            = {250   , Name StrCat[pp3,"/3rz"]},
+    flag_mat_scat = { 0    , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_Scat   = {2.25  , Name StrCat[pp3,"/5Custom real part of relative permittivity"]},
+    eps_im_Scat   = { 0    , Name StrCat[pp3,"/6Custom real part of relative permittivity"]},
+
+    flag_mat_1    = { 0    , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_2    = { 0    , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_3    = { 0    , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_4    = { 0    , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_5    = { 0    , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_6    = { 0    , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_L_1    = {1     , Name StrCat[pp4,"/layer 1: real part of relative permittivity"]},
+    eps_im_L_1    = {0     , Name StrCat[pp4,"/layer 1: imag part of relative permittivity"]},
+    eps_re_L_2    = {1     , Name StrCat[pp4,"/layer 2: real part of relative permittivity"]},
+    eps_im_L_2    = {0     , Name StrCat[pp4,"/layer 2: imag part of relative permittivity"]},
+    eps_re_L_3    = {1     , Name StrCat[pp4,"/layer 3: real part of relative permittivity"]},
+    eps_im_L_3    = {0     , Name StrCat[pp4,"/layer 3: imag part of relative permittivity"]},
+    eps_re_L_4    = {2.25  , Name StrCat[pp4,"/layer 4: real part of relative permittivity"]},
+    eps_im_L_4    = {0     , Name StrCat[pp4,"/layer 4: imag part of relative permittivity"]},
+    eps_re_L_5    = {2.25  , Name StrCat[pp4,"/layer 5: real part of relative permittivity"]},
+    eps_im_L_5    = {0     , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
+    eps_re_L_6    = {2.25  , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+    eps_im_L_6    = {0     , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
+
+    paramaille    = {8     , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+    scat_lc       = {lambda0/(1.5*paramaille)    , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+    PML_top       = {lambda0  , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+    PML_bot       = {lambda0  , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+    Nmax          = {1     , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+    refine_mesh_L1= {1     , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+    refine_mesh_L2= {1     , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+    refine_mesh_L3= {1     , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    refine_mesh_L4= {1     , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+    refine_mesh_L5= {1     , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+    refine_mesh_L6= {1     , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+
+    InterpSampling     = { 50   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+    Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+    FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+    FlagOutEtotFull    = { 0    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutEscaFull    = { 0    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutPoyFull     = { 0    , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+]; 
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+    hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0   = phideg*Pi/180;
+psi0   = psideg*Pi/180;
+
+DomainZsizeSca  = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot  = PML_top_hh-hh_L_6;
+npts_interpX    = period_x/InterpSampling;
+npts_interpY    = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;

DiffractionGratings/grating3D_data_torus.geo

+nm  = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+    lambda0       = {1000  , Name StrCat[pp1,"/1lambda0 [nm]"]},
+    thetadeg      = {0     , Name StrCat[pp1,"/2theta0 [deg]"]},
+    phideg        = {0     , Name StrCat[pp1,"/3phi0 [deg]"]},
+    psideg        = {0     , Name StrCat[pp1,"/4psi0 [deg]"]},
+    period_x      = {300   , Name StrCat[pp2,"/1X period [nm]"]},
+    period_y      = {300   , Name StrCat[pp2,"/2Y period [nm]"]},
+    thick_L_1     = {50    , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
+    thick_L_2     = {100   , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    thick_L_3     = {100   , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+    thick_L_4     = {100   , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+    thick_L_5     = {50    , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+    thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+    
+    tag_geom      = {  3   , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal"}},
+    rx            = {150/2 , Name StrCat[pp3,"/1rx"]},
+    ry            = { 50   , Name StrCat[pp3,"/2ry"]},
+    rz            = { 25   , Name StrCat[pp3,"/3rz"]},
+    flag_mat_scat = { 0    , Name StrCat[pp3,"/4Scatterer permittivity model"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_Scat   = {-21   , Name StrCat[pp3,"/5Custom real part of relative permittivity"]},
+    eps_im_Scat   = { 20   , Name StrCat[pp3,"/6Custom real part of relative permittivity"]},
+
+    flag_mat_1    = { 0    , Name StrCat[pp4,"/1Layer 1"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_2    = { 0    , Name StrCat[pp4,"/2Layer 2"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_3    = { 0    , Name StrCat[pp4,"/3Layer 3"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_4    = { 0    , Name StrCat[pp4,"/4Layer 4"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_5    = { 0    , Name StrCat[pp4,"/5Layer 5"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    flag_mat_6    = { 0    , Name StrCat[pp4,"/6Layer 6"], Choices {0="Custom (Value Below)",1="SiO2",2="Ag (palik)",3="Al (palik)",4="Au (johnson)",5="Nb2O5",6="ZnSe",7="MgF2",8="TiO2",9="PMMA",10="Si",11="ITO",12="Cu (palik)"} },
+    eps_re_L_1    = {1     , Name StrCat[pp4,"/7Custom Values/layer 1: real part of relative permittivity"]},
+    eps_im_L_1    = {0     , Name StrCat[pp4,"/7Custom Values/layer 1: imag part of relative permittivity"]},
+    eps_re_L_2    = {1     , Name StrCat[pp4,"/7Custom Values/layer 2: real part of relative permittivity"]},
+    eps_im_L_2    = {0     , Name StrCat[pp4,"/7Custom Values/layer 2: imag part of relative permittivity"]},
+    eps_re_L_3    = {1     , Name StrCat[pp4,"/7Custom Values/layer 3: real part of relative permittivity"]},
+    eps_im_L_3    = {0     , Name StrCat[pp4,"/7Custom Values/layer 3: imag part of relative permittivity"]},
+    eps_re_L_4    = {2.25  , Name StrCat[pp4,"/7Custom Values/layer 4: real part of relative permittivity"]},
+    eps_im_L_4    = {0     , Name StrCat[pp4,"/7Custom Values/layer 4: imag part of relative permittivity"]},
+    eps_re_L_5    = {2.25  , Name StrCat[pp4,"/7Custom Values/layer 5: real part of relative permittivity"]},
+    eps_im_L_5    = {0     , Name StrCat[pp4,"/7Custom Values/layer 5: imag part of relative permittivity"]},
+    eps_re_L_6    = {2.25  , Name StrCat[pp4,"/7Custom Values/layer 6: real part of relative permittivity"]},
+    eps_im_L_6    = {0     , Name StrCat[pp4,"/7Custom Values/layer 6: imag part of relative permittivity"]},
+    
+    paramaille    = {12      , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+    scat_lc       = {10       , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+    PML_top       = {lambda0*1.5 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+    PML_bot       = {lambda0*1.5 , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
+    Nmax          = {1       , Name StrCat[pp5,"/6Max number of non specular order to output [-]"]},
+    refine_mesh_L1= {1       , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+    refine_mesh_L2= {1       , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+    refine_mesh_L3= {1       , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    refine_mesh_L4= {1       , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+    refine_mesh_L5= {1       , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+    refine_mesh_L6= {1       , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    
+    InterpSampling     = {15   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+    Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+    FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+    FlagOutEtotFull    = { 0    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutEscaFull    = { 0    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutPoyFull     = { 0    , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+hh_L_6 = -thick_L_6;
+For k In {1:5}
+    hh_L~{6-k} = hh_L~{7-k}+thick_L~{7-k};
+EndFor
+PML_bot_hh = hh_L_6-PML_bot;
+PML_top_hh = hh_L_1+thick_L_1;
+
+theta0 = thetadeg*Pi/180;
+phi0   = phideg*Pi/180;
+psi0   = psideg*Pi/180;
+
+DomainZsizeSca  = PML_top_hh+PML_bot-(hh_L_6-PML_bot);
+DomainZsizeTot  = PML_top_hh-hh_L_6;
+npts_interpX    = period_x/InterpSampling;
+npts_interpY    = period_y/InterpSampling;
+npts_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;

DiffractionGratings/grating3D_postplot.py

+import numpy as np
+myDir = 'run_results3D/'
+Rnm = np.loadtxt(myDir+'eff_r.txt')[:,1]
+Tnm = np.loadtxt(myDir+'eff_t.txt')[:,1]
+Q = [np.loadtxt(myDir+'temp-Q_L_%g.txt'%k)[1] for k in range(2,7)]
+Q.append(np.loadtxt(myDir+'temp-Q_scat.txt')[1])
+Q=np.array(Q)
+TOT = Rnm.sum()+Tnm.sum()+Q.sum()
+print('Rtot',Rnm.sum())
+print('Ttot',Tnm.sum())
+print('Atot',Q.sum())
+print('TOT ',TOT)
\ No newline at end of file

DiffractionGratings/grating2D_clean_all.sh → DiffractionGratings/gratings_clean_all.sh

 rm -rf ./run_results
+rm -rf ./run_results3D
 rm *.msh
 rm *.pre
 rm *.db