diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 850abe4a0e1861a25c70416a0cd61ab6c609d413..b596ca1f40e00b8d22ca9b6619118303a45262eb 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -221,7 +221,6 @@ EndFor
 If (tag_geom==3) // Split torus weird otherwise
   Mesh.Algorithm = 6;
 EndIf
-Mesh.Optimimze=1;
 // Mesh.SurfaceEdges = 0;
 Mesh.VolumeEdges = 0;
 Mesh.ElementOrder = og;
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 8af300f720f9be44a785ba8175920e04e36b8b2e..5c2d3b9eb419541d7ffce79f30c639899c2db156 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -31,10 +31,14 @@ Group {
 
   
 	SurfDirichlet  = Region[{401,402}];
-  SurfBloch      = Region[{}];
-  SurfExcludeFacets = Region[501];
-  // SurfBloch      = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
+  SurfBloch      = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
 
+  If (FlagLinkFacets==1)
+    SurfExcludeFacets  = Region[{}];
+  Else
+    SurfExcludeFacets  = Region[{SurfBloch}];
+  EndIf
+  
   L_1 = Region[{L_1_temp,SurfIntTop}];
   L_6 = Region[{L_6_temp,SurfIntBot}];
 
@@ -281,17 +285,17 @@ FunctionSpace {
       { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
       { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
       { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
-      If(oi==2)
+      If (FlagLinkFacets==1)
         { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochX; }
         { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint BlochY; }
-        { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
         { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochX; }
         { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint BlochY; }
-        { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
         { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochX; }
         { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint BlochY; }
-        { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
       EndIf
+      { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+      { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+      { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
     }
   }
 }
diff --git a/DiffractionGratings/grating3D_data_2Dlamellar.geo b/DiffractionGratings/grating3D_data_2Dlamellar.geo
index 396f680f689e94b9d94629dadc4f085dccd7ceb4..dd118042f673314cba98289b919d5230eb9e588a 100644
--- a/DiffractionGratings/grating3D_data_2Dlamellar.geo
+++ b/DiffractionGratings/grating3D_data_2Dlamellar.geo
@@ -59,6 +59,7 @@ DefineConstant[
     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 [-]"]},
+    FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0="1",1="2"}},
     
     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_data_bisin.geo b/DiffractionGratings/grating3D_data_bisin.geo
index 46c12218a1b51a34c6302aa0ec6b9e428b5aa87d..b414a938d19706af0ed1afe35906c4f26bfa258e 100644
--- a/DiffractionGratings/grating3D_data_bisin.geo
+++ b/DiffractionGratings/grating3D_data_bisin.geo
@@ -59,6 +59,7 @@ DefineConstant[
     refine_mesh_L_4= {1      , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
     refine_mesh_L_5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0="1",1="2"}},
     
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_checker.geo b/DiffractionGratings/grating3D_data_checker.geo
index 930636e652f49c6caa0e81bb8aebd859a3cfbea6..c670907f05ac40edf23bb0b84bd0b734c69e8892 100644
--- a/DiffractionGratings/grating3D_data_checker.geo
+++ b/DiffractionGratings/grating3D_data_checker.geo
@@ -59,6 +59,7 @@ DefineConstant[
     refine_mesh_L_4= {1      , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
     refine_mesh_L_5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0="1",1="2"}},
     
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_halfellipsoid.geo b/DiffractionGratings/grating3D_data_halfellipsoid.geo
index a2da196699807678544a49f3bad37160610786a4..98307e1792f54d0b739f90ea9c88714936c46e77 100644
--- a/DiffractionGratings/grating3D_data_halfellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_halfellipsoid.geo
@@ -59,6 +59,7 @@ DefineConstant[
     refine_mesh_L_4= {1     , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
     refine_mesh_L_5= {1     , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1     , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0="1",1="2"}},
         
     InterpSampling     = { 10   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_hole.geo b/DiffractionGratings/grating3D_data_hole.geo
index cbc63f305be89e6589b6895fdf24eceb26d5074a..3b42b6a778bba646f7670719c81245f597aebc95 100644
--- a/DiffractionGratings/grating3D_data_hole.geo
+++ b/DiffractionGratings/grating3D_data_hole.geo
@@ -59,6 +59,7 @@ DefineConstant[
     refine_mesh_L_4= {1      , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
     refine_mesh_L_5= {1      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0="1",1="2"}},
     
     InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_pyramid.geo b/DiffractionGratings/grating3D_data_pyramid.geo
index 7a05d4935b1d25e45ef14ec93c2688ced15bf811..bab8d7f7d10ebbe1d72b58edd4a52a7513ecd342 100644
--- a/DiffractionGratings/grating3D_data_pyramid.geo
+++ b/DiffractionGratings/grating3D_data_pyramid.geo
@@ -59,6 +59,7 @@ DefineConstant[
     refine_mesh_L_4= {1     , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
     refine_mesh_L_5= {1     , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1     , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0="1",1="2"}},
 
     InterpSampling     = { 50   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_solarcell.geo b/DiffractionGratings/grating3D_data_solarcell.geo
index 86ec842cbce60c5ed0df24af0dabb371b00b3bb9..6b36721af529ee2f8b87b770f217d383c709fdbc 100644
--- a/DiffractionGratings/grating3D_data_solarcell.geo
+++ b/DiffractionGratings/grating3D_data_solarcell.geo
@@ -61,6 +61,8 @@ DefineConstant[
     refine_mesh_L_4= {2      , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
     refine_mesh_L_5= {2      , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {2      , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    FlagLinkFacets = {0      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], Choices {0="1",1="2"}},
+
     
     InterpSampling     = { 20   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
     Flag_interp_cubic  = { 1    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
diff --git a/DiffractionGratings/grating3D_data_torus.geo b/DiffractionGratings/grating3D_data_torus.geo
index 9f8bca3bd5ec9571411b23b411f3ce2c45582672..509ce89366b2f4b3ccfdd9d3664a0519225610d8 100644
--- a/DiffractionGratings/grating3D_data_torus.geo
+++ b/DiffractionGratings/grating3D_data_torus.geo
@@ -59,6 +59,7 @@ DefineConstant[
     refine_mesh_L_4= {1       , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
     refine_mesh_L_5= {1       , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
     refine_mesh_L_6= {1       , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    FlagLinkFacets = {1      , Name StrCat[pp5,"/8FlagLinkFacets? [-]"], 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} },