diff --git a/DiffractionGratings/grating2D.geo b/DiffractionGratings/grating2D.geo
index 4cc2d43d4e9f01ca284906114773393de90b1fdd..f79812265015fe37ef2906aa2167977be565e7ca 100644
--- a/DiffractionGratings/grating2D.geo
+++ b/DiffractionGratings/grating2D.geo
@@ -2,7 +2,11 @@
 // Author : Guillaume Demesy  //
 ////////////////////////////////
 
-Include "grating2D_data.geo";
+// Include "grating2D_data_AnisotropicGrating.geo";
+Include "grating2D_data_LamellarGrating.geo";
+// Include "grating2D_data_PhotonicCrystalSlab.geo";
+// Include "grating2D_data_ResonantGrating.geo";
+// Include "grating2D_data.geo";
 
 paramaille_rods = lambda_min*nm/(paramaille*paramaille_scale_rods);
 lc_rod_out      = lambda_min*nm/(paramaille*paramaille_scale_rod_out);
diff --git a/DiffractionGratings/grating2D.pro b/DiffractionGratings/grating2D.pro
index 34d1292c487384c23e13287ea1449052f5424214..2a2d13317923cce950032b0851045dc6d5d1cf2c 100644
--- a/DiffractionGratings/grating2D.pro
+++ b/DiffractionGratings/grating2D.pro
@@ -2,9 +2,14 @@
 // Author : Guillaume Demesy  //
 ////////////////////////////////
 
-Include "grating2D_data.geo";
+// Include "grating2D_data_AnisotropicGrating.geo";
+Include "grating2D_data_LamellarGrating.geo";
+// Include "grating2D_data_PhotonicCrystalSlab.geo";
+// Include "grating2D_data_ResonantGrating.geo";
+// Include "grating2D_data.geo";
 
 Include "grating2D_materials.pro";
+
 myDir = "run_results/";
 DefineConstant[
   lambda0 = {lambda_min , Min lambda_min, Max lambda_max, Step (lambda_max-lambda_min)/(nb_lambdas-1), Name StrCat[pp2, "0wavelength [nm]"] , Loop 1, Highlight Str[colorpp2],Graph "200000200020", ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}
@@ -371,6 +376,7 @@ PostProcessing {
         { Name Hz_tot       ; Value { Local { [ {u2d}+u1[]          ]; In Omega; Jacobian JVol; } } }
         { Name NormHz_tot   ; Value { Local { [ Norm[{u2d}+u1[]]    ]; In Omega; Jacobian JVol; } } }
         { Name E1           ; Value { Local { [ E1[]                ]; In Omega; Jacobian JVol; } } }
+        { Name u1           ; Value { Local { [ u1[]                ]; In Omega; Jacobian JVol; } } }
         { Name Hz_totp1     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
         { Name Hz_totp2     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
         { Name Hz_totp3     ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
@@ -465,6 +471,9 @@ PostOperation {
         Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"]];
         Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"]];
         Print[ Hz_tot    , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
+        // debug
+        Print[ u1    , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u1%.2fnm.pos", lambda0/nm)];
+        Print[ E1    , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("E1%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("E1%.2fnm.pos", lambda0/nm)];
         Echo[ Str["For i In {0:2}",
                   "  View[i].LineWidth = 4;View[i].ColormapNumber = 15;View[0].AxesFormatX = '%.1f';",
                   "EndFor"], File StrCat[myDir,"tmp0.geo"]] ;
diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index a0538ea679a502a031b818520deb0f876618ae40..6f2153a6a9f51d37a61bee1bc42a52e374e4d4d1 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -99,12 +99,12 @@ If (tag_geom==6)
   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);
+L1_n     = Sqrt(Abs(eps_re_L_1));
+L2_n     = Sqrt(Abs(eps_re_L_2));
+L3_n     = Sqrt(Abs(eps_re_L_3));
+L4_n     = Sqrt(Abs(eps_re_L_4));
+L5_n     = Sqrt(Abs(eps_re_L_5));
+L6_n     = Sqrt(Abs(eps_re_L_6));
 Scat_n   = 2 ;
 
 PML_lc     = lambda_m/(paramaille*.5);
@@ -132,12 +132,12 @@ If (tag_geom==1)
   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};
+  Line(100) = {31, 65};2
+  Curve Loop(49) = {38, -982, -97};
+  Plane Surface(49) = {49};2
+  Curve Loop(50) = {98, 48,2 100};
+  Plane Surface(50) = {50};2
+  Curve Loop(51) = {99, 42,2 100};
   Plane Surface(51) = {51};
   Curve Loop(52) = {99, -46, 97};
   Plane Surface(52) = {52};
@@ -166,7 +166,9 @@ If (tag_geom==5)
   Rotate {{0, 0, 1}, {0,0,0}, Pi/4} {Volume{9};}
 EndIf
 
-
+If (tag_geom==7)
+  Box(9) = {-rx/2,-period_y/2, hh_L_3, rx, period_y, rz};
+EndIf
 
 
 Coherence;
@@ -176,10 +178,10 @@ 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_L5"      ,3) = {3};
+Physical Volume("LAYER_L4"      ,4) = {4};
+Physical Volume("LAYER_L3"      ,5) = {10};
+Physical Volume("LAYER_L2"      ,6) = {6};
 Physical Volume("LAYER_L1_SUPER",7) = {7};
 Physical Volume("PMLTOP"        ,8) = {8};
 Physical Volume("SCAT"          ,9) = {9};
@@ -203,17 +205,18 @@ 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;
+Characteristic Length{:}              = PML_lc/3;
+// 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_L6()} = L6_lc/refine_mesh_L6;
+// Characteristic Length{pts_LAYER_L5()} = L5_lc/refine_mesh_L5;
+// Characteristic Length{pts_LAYER_L4()} = L4_lc/refine_mesh_L4;
+// 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.Optimimze=1;
+// Mesh.SurfaceEdges = 0;
 Mesh.VolumeEdges = 0;
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 30eec74a18c55728ad67456ff3dd8258cc9c910d..1d1d19eae762f70c46aef202ee4dc90a563c88da 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -5,7 +5,7 @@
 // test_case = "checker";
 // test_case = "bisin";
 Include StrCat["grating3D_data_",test_case,".geo"];
-Include "gratings_common_materials.pro"
+Include "grating3D_materials.pro"
 
 myDir = "run_results3D/";
 
@@ -124,6 +124,10 @@ Function{
 	mur[PMLtop]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 	epsr[PMLbot]  = Re[epsr2[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 	mur[PMLbot]   =             TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+
+	// epsr[PMLbot]  = epsr2[];
+	// mur[PMLbot]   = TensorDiag[1,1,1];
+  
 	// epsr[PMLtop]     = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
 	// mur[PMLtop]      = TensorDiag[sz_bermutop[]*sy/sx,sx*sz_bermutop[]/sy,sx*sy/sz_bermutop[]];
 	// epsr[PMLbot]     = TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
@@ -317,6 +321,7 @@ PostProcessing {
       { Name Etot   ; Value { Local { [ {u}+E1[]  ]; 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 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; } } }
@@ -396,6 +401,8 @@ PostOperation {
         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"]];
+      Print [ Etot , OnElementsOf Omega, File StrCat[myDir,"Etot.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"]];