From f38d523cce1be5370efdf13fd802e472ad460920 Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Tue, 19 Nov 2019 20:07:20 +0100
Subject: [PATCH] check all

---
 DiffractionGratings/grating3D.geo            |  9 ++++-
 DiffractionGratings/grating3D.pro            | 42 +++++++++++---------
 DiffractionGratings/grating3D_data_torus.geo | 26 ++++++------
 3 files changed, 44 insertions(+), 33 deletions(-)

diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index e028fcb..a0538ea 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -1,3 +1,9 @@
+// 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");
@@ -144,7 +150,7 @@ If (tag_geom==2)
 EndIf
 
 If (tag_geom==3)
-  Torus(9) = {0,0,hh_L_3+ry-2,rx,ry};
+  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
@@ -208,5 +214,6 @@ 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;
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index f6e0ad8..30eec74 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -1,10 +1,14 @@
-myDir = "run_results3D/";
+// 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"];
-// Include "grating3D_data_hole.geo";
-// Include "grating3D_data_pyramid.geo";
-// Include "grating3D_data_halfellipsoid.geo";
-
 Include "gratings_common_materials.pro"
+
+myDir = "run_results3D/";
+
 Group {
 
 	// SubDomains
@@ -358,18 +362,6 @@ PostProcessing {
 PostOperation {
   { Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
     Operation {
-      If (FlagOutEscaCuts==1)
-        Print [ u , OnPlane { {-period_x/2,0,hh_L_6-PML_bot} {period_x/2,0,hh_L_6-PML_bot} {-period_x/2,0,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpZSca} , File StrCat[myDir,"u_cut_Y=0.pos"], Name "u_cut_Y=0"];
-        Print [ u , OnPlane { {0,-period_y/2,hh_L_6-PML_bot} {0,period_y/2,hh_L_6-PML_bot} {0,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_cut_X=0.pos"], Name "u_cut_X=0"];
-      EndIf
-      If (FlagOutEtotCuts==1)
-        Print [ Etot , 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,"Etot_cut_Y=0.pos"], Name "Etot_cut_Y=0"];
-        Print [ Etot , 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,"Etot_cut_X=0.pos"], Name "Etot_cut_X=0"];
-      EndIf
-      If (FlagOutPoyCut==1)
-        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
       If (FlagOutEscaFull==1)
         If (Flag_interp_cubic==1)
           Print [ u , OnBox { {-period_x/2,-period_y/2,hh_L_6-PML_bot} {period_x/2,-period_y/2,hh_L_6-PML_bot} {-period_x/2,period_y/2,hh_L_6-PML_bot} {-period_x/2,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_grid.pos"], Name "u_grid"];
@@ -386,13 +378,25 @@ PostOperation {
       EndIf
       If (FlagOutPoyFull==1)
         If (Flag_interp_cubic==1)
-          Print [ Poy_tot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_grid.pos", Name "Poy_tot_grid"];
+          Print [ Poy_tot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File StrCat[myDir,"Poy_tot_grid.pos"], Name "Poy_tot_grid"];
         Else
           Print [ Poy_tot , OnElementsOf Omega_plot, File StrCat[myDir,"Poy_tot.pos"]];
         EndIf
       EndIf
+      If (FlagOutEscaCuts==1)
+        Print [ u , OnPlane { {-period_x/2,0,hh_L_6-PML_bot} {period_x/2,0,hh_L_6-PML_bot} {-period_x/2,0,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpZSca} , File StrCat[myDir,"u_cut_Y=0.pos"], Name "u_cut_Y=0"];
+        Print [ u , OnPlane { {0,-period_y/2,hh_L_6-PML_bot} {0,period_y/2,hh_L_6-PML_bot} {0,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpY,npts_interpZSca} , File StrCat[myDir,"u_cut_X=0.pos"], Name "u_cut_X=0"];
+      EndIf
+      If (FlagOutEtotCuts==1)
+        Print [ Etot , 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,"Etot_cut_Y=0.pos"], Name "Etot_cut_Y=0"];
+        Print [ Etot , 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,"Etot_cut_X=0.pos"], Name "Etot_cut_X=0"];
+      EndIf
+      If (FlagOutPoyCut==1)
+        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 [ epsr_xx , OnElementsOf Omega, File StrCat[myDir,"epsr_xx.pos"]];
+      // Print [ epsr_xx , OnElementsOf Omega, File StrCat[myDir,"epsr_xx.pos"]];
       // Print [ epsr_0xx , OnElementsOf Scat, File StrCat[myDir,"epsr_xx.pos"]];
       // Print [ Ecm , OnElementsOf Omega_plot, File StrCat[myDir,"Ecm.pos"]];
       // Print [ E1, OnElementsOf Omega_plot, File StrCat[myDir,"E1.pos"]];
diff --git a/DiffractionGratings/grating3D_data_torus.geo b/DiffractionGratings/grating3D_data_torus.geo
index a5216cd..70fb906 100644
--- a/DiffractionGratings/grating3D_data_torus.geo
+++ b/DiffractionGratings/grating3D_data_torus.geo
@@ -10,12 +10,12 @@ DefineConstant[
     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]"]},
+    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     = {200   , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    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     = {200   , Name StrCat[pp2,"/6thickness layer 4 [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)"]},
     
@@ -46,10 +46,10 @@ DefineConstant[
     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    = {10      , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
-    scat_lc       = {20      , 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]"]},
+    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 [-]"]},
@@ -58,14 +58,14 @@ DefineConstant[
     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} },
+    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    = { 1    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
-    FlagOutEscaFull    = { 1    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
-    FlagOutPoyFull     = { 1    , Name StrCat[pp6,"/7Poynting Full Output?"] , 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;
-- 
GitLab