diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 27e27170362032ee4a0e98ca2f62cafd880c637d..c2b429d36dece9d30532d375f0bb55006e9a0c82 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -151,7 +151,7 @@ If (tag_geom==5)
EndIf
If (tag_geom==7)
- Box(9) = {-rx/2,-period_y/2, hh_L_3, rx, period_y, rz};
+ Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz};
EndIf
@@ -209,11 +209,14 @@ EndIf
// Start with coarsest
Characteristic Length{:} = lc_PML;
+
For k In {0:8}
which_dom = meshing_sequence(k);
Characteristic Length{PointsOf{Physical Volume{which_dom};}} = list_lc(which_dom-1);
EndFor
- If (tag_geom==3) // Split torus weird otherwise
+
+
+If (tag_geom==3) // Split torus weird otherwise
Mesh.Algorithm = 6;
EndIf
// Mesh.Optimimze3D=1;
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index c82f7ed6f2f1f6ecfd74066fc4dc66d8d28beeca..4cc35dca96827ceaf10ed837a7d8a1205069855d 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -181,6 +181,35 @@ Function{
H1d[Omega_super] = Hr[];
H1d[Omega_subs] = Ht[];
+ // ////////////////////////:
+ // k1s[] = (Abs[theta0]>1e-12) ? k1[]/\zhat[] / Norm[ k1[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+ // k2s[] = (Abs[theta0]>1e-12) ? Re[k2[]]/\zhat[] / Norm[ Re[k2[]]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+ // k1rs[] = (Abs[theta0]>1e-12) ? k1r[]/\zhat[] / Norm[ k1r[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+ // k1p[] = (Abs[theta0]>1e-12) ? k1[]/\k1s[] / Norm[ k1[]/\k1s[]] : Vector[ Cos[phi0], Sin[phi0], 0.] ;
+ // k2p[] = (Abs[theta0]>1e-12) ? Re[k2[]]/\k2s[] / Norm[ Re[k2[]]/\k2s[]] : Vector[ Cos[phi0], Sin[phi0], 0.] ;
+ // k1rp[] = (Abs[theta0]>1e-12) ? k1r[]/\k1rs[] / Norm[ k1r[]/\k1rs[]] : Vector[-Cos[phi0],-Sin[phi0], 0.] ;
+ // rsa[] = (CompZ[k1[]]-CompZ[k2[]])/(CompZ[k1[]]+CompZ[k2[]]);
+ // tsa[] = 2.*CompZ[k1[]] /(CompZ[k1[]]+CompZ[k2[]]);
+ // rpa[] = (CompZ[k1[]]*epsr2[]-CompZ[k2[]]*epsr1[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
+ // tpa[] = (2.*CompZ[k1[]]*epsr2[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]) * Sqrt[epsr1[]/epsr2[]];
+ // AmplEi[] = Cos[psi0]* k1p[] + Sin[psi0]* k1s[];
+ // AmplEr[] = rpa[]*Cos[psi0]*k1rp[] + rsa[]*Sin[psi0]*k1rs[];
+ // AmplEt[] = tpa[]*Cos[psi0]* k2p[] + tsa[]*Sin[psi0]* k2s[];
+ // Eia[] = Ae*AmplEi[] * Exp[I[]*k1[]*XYZ[]];
+ // Eta[] = Ae*AmplEt[] * Exp[I[]*k2[]*XYZ[]];
+ // Era[] = Ae*AmplEr[] * Exp[I[]*k1r[]*XYZ[]];
+ // Hia[] = 1/(om0*mu0*mur[]) * k1[] /\ Ei[];
+ // Hra[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
+ // Hta[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
+ // E1a[Omega_super] = Eia[]+Era[];
+ // E1a[Omega_subs] = Eta[];
+ // E1da[Omega_super] = Era[];
+ // E1da[Omega_subs] = Eta[];
+ // H1a[Omega_super] = Hia[]+Hra[];
+ // H1a[Omega_subs] = Hta[];
+ // H1da[Omega_super] = Hra[];
+ // H1da[Omega_subs] = Hta[];
+ // /////////////////////
source[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
// Bloch phase shifts
@@ -326,9 +355,10 @@ PostProcessing {
{ Name u ; Value { Local { [ {u} ]; In Omega; Jacobian JVol; } } }
{ Name Etot ; Value { Local { [ {u}+E1[] ]; In Omega; Jacobian JVol; } } }
{ Name Htot ; Value { Local { [ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]; In Omega; Jacobian JVol; } } }
+ { Name Htotx ; Value { Local { [CompX[ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]]; In Omega; Jacobian JVol; } } }
{ Name Edif ; Value { Local { [ {u}+E1d[] ]; In Omega; Jacobian JVol; } } }
{ Name E1 ; Value { Local { [ E1[] ]; In Omega; Jacobian JVol; } } }
- { Name H1y ; Value { Local { [CompY[H1[]] ]; In Omega; Jacobian JVol; } } }
+ { Name H1x ; Value { Local { [CompX[H1[]] ]; In Omega; Jacobian JVol; } } }
{ Name Em ; Value { Local { [ E1d[] ]; In Omega; Jacobian JVol; } } }
{ Name source ; Value { Local { [ source[] ]; In Omega; Jacobian JVol; } } }
{ Name epsr_xx; Value { Local { [ CompXX[epsr[]] ]; In Omega; Jacobian JVol; } } }
@@ -408,9 +438,10 @@ PostOperation {
Print [ Poy_tot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File StrCat[myDir,"Poy_tot_cut_Y=0.pos"], Name "Poy_tot_cut_Y=0"];
Print [ Poy_tot , OnPlane { {0,-period_y/2,hh_L_6} {0,period_y/2,hh_L_6} {0,-period_y/2,hh_L_1+thick_L_1} } {npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_cut_X=0.pos"], Name "Poy_tot_cut_X=0"];
EndIf
- // // // For DEBUG
- // Print [ H1y , OnElementsOf Omega, File StrCat[myDir,"H1y.pos"]];
+ // // // // For DEBUG
+ // Print [ H1x , OnElementsOf Omega, File StrCat[myDir,"H1x.pos"]];
// Print [ Etot , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
+ // Print [ Htotx , OnElementsOf Omega, File StrCat[myDir,"Htotx.pos"]];
For k In {2:6}
Print[ Abs_L~{k}[L~{k}], OnGlobal, File > StrCat[myDir,Sprintf["temp-Q_L_%g.txt",k]], Format Table ];
EndFor
diff --git a/DiffractionGratings/grating3D_data_2Dlamellar.geo b/DiffractionGratings/grating3D_data_2Dlamellar.geo
index 2572da07b80271d0ff0550f892bf4f8eb075fc45..69ccbe5e93b0cbf3117c8e8eb4ef9bc15c04ad57 100644
--- a/DiffractionGratings/grating3D_data_2Dlamellar.geo
+++ b/DiffractionGratings/grating3D_data_2Dlamellar.geo
@@ -8,20 +8,20 @@ pp6 = "6Output";
DefineConstant[
lambda0 = {1000 , Name StrCat[pp1,"/1lambda0 [nm]"]},
thetadeg = {30 , Name StrCat[pp1,"/2theta0 [deg]"]},
- phideg = {0 , Name StrCat[pp1,"/3phi0 [deg]"]},
+ phideg = {90 , Name StrCat[pp1,"/3phi0 [deg]"]},
psideg = {0 , Name StrCat[pp1,"/4psi0 [deg]"]},
- period_x = {1000 , Name StrCat[pp2,"/1X period [nm]"]},
- period_y = {20 , Name StrCat[pp2,"/2Y period [nm]"]},
+ period_x = {10 , Name StrCat[pp2,"/1X period [nm]"]},
+ period_y = {1000 , Name StrCat[pp2,"/2Y period [nm]"]},
thick_L_1 = {200 , Name StrCat[pp2,"/3thickness layer 1 [nm] (superstrate)"]},
thick_L_2 = {200 , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
thick_L_3 = {1050 , 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 = {100 , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
+ thick_L_4 = {20 , Name StrCat[pp2,"/6thickness layer 4 [nm]"]},
+ thick_L_5 = {20 , Name StrCat[pp2,"/7thickness layer 5 [nm]"]},
+ thick_L_6 = {50 , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
tag_geom = { 7 , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
- rx = {500 , Name StrCat[pp3,"/1rx"]},
- ry = {0 , Name StrCat[pp3,"/2ry"]},
+ rx = { 0 , Name StrCat[pp3,"/1rx"]},
+ ry = {500 , Name StrCat[pp3,"/2ry"]},
rz = {1000 , 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 = {-44.9757 , Name StrCat[pp3,"/7eps_re_Scat"]},
@@ -45,20 +45,26 @@ DefineConstant[
eps_im_L_5 = {2.9524 , Name StrCat[pp4,"/layer 5: imag part of relative permittivity"]},
eps_re_L_6 = {-44.9757 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
eps_im_L_6 = {2.9524 , Name StrCat[pp4,"/layer 6: 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 = {-44.9757 , Name StrCat[pp4,"/layer 6: real part of relative permittivity"]},
+ // eps_im_L_6 = {2.9524 , Name StrCat[pp4,"/layer 6: imag part of relative permittivity"]},
og = {0 , Name StrCat[pp5,"/0geometrical order [-]"] , Choices {0="1",1="2"} },
oi = {1 , Name StrCat[pp5,"/0interpolation order [-]"], Choices {0="1",1="2"} },
- paramaille = {10 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
- lc_scat = {30 , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
- PML_top = {lambda0 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
+ paramaille = {20 , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+ lc_scat = {15 , Name StrCat[pp5,"/2Scatterer absolute mesh size [nm]"]},
+ PML_top = {lambda0*1 , Name StrCat[pp5,"/4PML top thickness [nm]"]},
PML_bot = {lambda0/5 , Name StrCat[pp5,"/5PML bot thickness [nm]"]},
Nmax = {1 , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
refine_mesh_L_1= {1 , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
refine_mesh_L_2= {1 , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
refine_mesh_L_3= {1 , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
- refine_mesh_L_4= {0.5 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
- refine_mesh_L_5= {0.5 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
- refine_mesh_L_6= {0.5 , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+ refine_mesh_L_4= {0.3 , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+ refine_mesh_L_5= {0.3 , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+ refine_mesh_L_6= {0.3 , 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} },
diff --git a/DiffractionGratings/grating3D_runall.sh b/DiffractionGratings/grating3D_runall.sh
index d0f5ed8346e3eb016b49dc492d824a20b70b99c1..6aa7ab871474bc4d059ec59b05560ec28bcd0041 100644
--- a/DiffractionGratings/grating3D_runall.sh
+++ b/DiffractionGratings/grating3D_runall.sh
@@ -1,9 +1,9 @@
GMSHDIR="/home/demesy/programs/gmsh-4.4.1-Linux64/bin"
-# for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar
-# do
-# $GMSHDIR/gmsh grating3D.pro -setstring test_case $t -
-# mv res3D res3D_$t
-# done
+for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar
+do
+ $GMSHDIR/gmsh grating3D.pro -setstring test_case $t -
+ mv res3D res3D_$t
+done
for t in bisin checker halfellipsoid hole pyramid torus 2Dlamellar
do