From d5df061f4d8ce0053834df6c245ddc33069ebd85 Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Wed, 1 Feb 2023 18:32:38 +0100
Subject: [PATCH] update split torus to U-shape

---
 DiffractionGratings/README.txt                |   2 +-
 DiffractionGratings/grating2D_scalar.pro      |   2 +
 DiffractionGratings/grating3D.geo             |  16 ++-
 DiffractionGratings/grating3D.pro             |  13 +-
 DiffractionGratings/grating3D_data.geo        |   4 +-
 .../grating3D_data_bi_sinusoidal.geo          |   2 +-
 .../grating3D_data_checker.geo                |   2 +-
 .../grating3D_data_convergence.geo            |   2 +-
 .../grating3D_data_half_ellipsoid.geo         |   2 +-
 DiffractionGratings/grating3D_data_hole.geo   |   2 +-
 .../grating3D_data_nanowire_solarcell.geo     |   2 +-
 .../grating3D_data_pyramid.geo                |   2 +-
 .../grating3D_data_retrieve_2D_lamellar.geo   |   2 +-
 .../grating3D_data_skewed_lattice.geo         |   2 +-
 DiffractionGratings/grating3D_data_torus.geo  | 125 ------------------
 .../grating3D_parallel_Mmatrix.sh             |   4 +-
 DiffractionGratings/grating3D_runall.sh       |   4 +-
 .../scattererTmatrix.geo                      |   2 +-
 .../scattererTmatrix_data.geo                 |  10 +-
 ElectromagneticScattering/scattering_data.geo |  10 +-
 20 files changed, 53 insertions(+), 157 deletions(-)
 delete mode 100644 DiffractionGratings/grating3D_data_torus.geo

diff --git a/DiffractionGratings/README.txt b/DiffractionGratings/README.txt
index 51f6270..9c4559e 100644
--- a/DiffractionGratings/README.txt
+++ b/DiffractionGratings/README.txt
@@ -11,7 +11,7 @@ Quick start 3D
 To run the test cases, type in terminal :
 gmsh grating3D.pro -setstring test_case pyramid
 gmsh grating3D.pro -setstring test_case hole
-gmsh grating3D.pro -setstring test_case torus
+gmsh grating3D.pro -setstring test_case U
 gmsh grating3D.pro -setstring test_case half_ellipsoid
 gmsh grating3D.pro -setstring test_case checker
 gmsh grating3D.pro -setstring test_case bisin
diff --git a/DiffractionGratings/grating2D_scalar.pro b/DiffractionGratings/grating2D_scalar.pro
index 740f7d8..44cdaae 100644
--- a/DiffractionGratings/grating2D_scalar.pro
+++ b/DiffractionGratings/grating2D_scalar.pro
@@ -372,6 +372,7 @@ Resolution {
 PostProcessing {
   { Name postpro_energy; NameOfFormulation helmoltz_scalar;
     Quantity {
+      { Name SSize      ; Value { Local { [ $KSPSystemSize  ]; In Omega; Jacobian JVol; } } }
       { Name u           ; Value { Local { [ {u2d}       ]; In Omega; Jacobian JVol; } } }
       { Name u_diff      ; Value { Local { [ {u2d}+u1d[] ]; In Omega; Jacobian JVol; } } }
       { Name u_tot       ; Value { Local { [ {u2d}+u1[]  ]; In Omega; Jacobian JVol; } } }
@@ -422,6 +423,7 @@ PostProcessing {
 PostOperation {
   { Name postop_energy; NameOfPostProcessing postpro_energy ;
     Operation {
+      Print[ SSize, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "SSize.txt"]] ;
       Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
       For i In {0:2*nb_orders}
         Print[ s_r~{i}[SurfCutSuper1], OnGlobal, StoreInVariable $sr~{i}, Format Table , File > StrCat[myDir, Sprintf("temp_s_r_%g.txt", i-nb_orders)]];
diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index a2671b5..36d3e59 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -186,9 +186,17 @@ If (tag_geom==2)
 EndIf
 
 If (tag_geom==3)
-  Torus(9) = {0,0,hh_L_3+ry-1*nm,rx,ry};
-  Dilate { { 0,0,hh_L_3 }, { 1, 1, rz/ry } } { Volume{9}; }
-  BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+  Lx = 210*nm;
+  Ly = 232*nm;
+  Box(9)  = {-Lx/2   , -Ly/2    , hh_L_3 , Lx, ry, rz};
+  Box(10) = {-Lx/2   ,  Ly/2-ry , hh_L_3 , Lx, ry, rz};
+  Box(11) = {-Lx/2   , -Ly/2    , hh_L_3 , rx, Ly, rz};
+  r_fillet = 15*nm;
+  scat() = BooleanUnion{ Volume{9}; Delete; }{ Volume{10,11}; Delete;} ;
+  fscat() = Abs(Boundary{ Volume{scat()}  ; });
+  escat() = Unique(Abs(Boundary{ Surface{fscat()}; }));
+  Fillet{scat()}{escat()}{r_fillet}
+
 EndIf
 
 If (tag_geom==4)
@@ -283,7 +291,7 @@ For k In {0:8}
 EndFor
 
 
-If (tag_geom==3) // Split torus weird otherwise
+If (tag_geom==3) // Split U weird otherwise
   Mesh.Algorithm = 6;
 EndIf
 // Mesh.SurfaceEdges = 0;
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 64f8f70..eb066f2 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -438,9 +438,17 @@ PostProcessing {
         Vector[$int_x_r~{ispecular}~{jspecular},
                $int_y_r~{ispecular}~{jspecular},
                $int_z_r~{ispecular}~{jspecular}] ] ; In SurfIntTop ; } } }
+      { Name et_specular ; Value { Term{Type Global; [ 
+        Vector[$int_x_t~{ispecular}~{jspecular},
+               $int_y_t~{ispecular}~{jspecular},
+               $int_z_t~{ispecular}~{jspecular}] ] ; In SurfIntBot ; } } }
+    
       // Project er_specular on the (s,p) basis
       { Name rp ; Value { Term{ Type Global; [ ppol[] * $er_specular]   ; In SurfIntTop ; } } }
       { Name rs ; Value { Term{ Type Global; [ spol[] * $er_specular]   ; In SurfIntTop ; } } }
+      { Name tp ; Value { Term{ Type Global; [ ppol[] * $et_specular]   ; In SurfIntBot ; } } }
+      { Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular]   ; In SurfIntBot ; } } }
+
 }
   }
 }
@@ -524,8 +532,11 @@ PostOperation {
         EndFor
       EndFor
       Print[ er_specular[SurfIntTop]  , OnRegion SurfIntTop, StoreInVariable $er_specular, Format Table];
+      Print[ et_specular[SurfIntBot]  , OnRegion SurfIntBot, StoreInVariable $et_specular, Format Table];
       Print[ rp[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir,"rp.txt"], Format Table ];
-      Print[ rs[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir,"rs.txt"], Format Table ];      
+      Print[ rs[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir,"rs.txt"], Format Table ];
+      Print[ tp[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"tp.txt"], Format Table ];
+      Print[ ts[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"ts.txt"], Format Table ];
       For i In {0:Nb_ordre-1}
         For j In {0:Nb_ordre-1}
           Print[ numbering_ij~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"numbering_ij.txt"], Format Table ];
diff --git a/DiffractionGratings/grating3D_data.geo b/DiffractionGratings/grating3D_data.geo
index 4862669..a82fe6f 100644
--- a/DiffractionGratings/grating3D_data.geo
+++ b/DiffractionGratings/grating3D_data.geo
@@ -1,8 +1,8 @@
 
 DefineConstant[
-  test_case = {"half_ellipsoid",
+  test_case = {"U",
     Name "1Geometry",
-    Choices {"half_ellipsoid", "hole", "pyramid", "torus", "bi_sinusoidal", "retrieve_2D_lamellar", "nanowire_solarcell", "convergence", "skewed_lattice"},
+    Choices {"half_ellipsoid", "hole", "pyramid", "U", "bi_sinusoidal", "retrieve_2D_lamellar", "nanowire_solarcell", "convergence", "skewed_lattice"},
     GmshOption "Reset", Autocheck 0
   }
 ];
diff --git a/DiffractionGratings/grating3D_data_bi_sinusoidal.geo b/DiffractionGratings/grating3D_data_bi_sinusoidal.geo
index 216ef0b..3fe136d 100644
--- a/DiffractionGratings/grating3D_data_bi_sinusoidal.geo
+++ b/DiffractionGratings/grating3D_data_bi_sinusoidal.geo
@@ -21,7 +21,7 @@ DefineConstant[
     thick_L_6     = {200   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = {  6      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = {  6      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
     ry            = {1.25*lambda0, Name StrCat[pp3,"/2ry"]},
     rz            = { 200    , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_checker.geo b/DiffractionGratings/grating3D_data_checker.geo
index 2e43fc2..2a34d7d 100644
--- a/DiffractionGratings/grating3D_data_checker.geo
+++ b/DiffractionGratings/grating3D_data_checker.geo
@@ -21,7 +21,7 @@ DefineConstant[
     thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = {  5      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = {  5      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {1.25*lambda0, Name StrCat[pp3,"/1rx"]},
     ry            = {1.25*lambda0, Name StrCat[pp3,"/2ry"]},
     rz            = {lambda0     , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_convergence.geo b/DiffractionGratings/grating3D_data_convergence.geo
index fff7de0..ded4ddd 100644
--- a/DiffractionGratings/grating3D_data_convergence.geo
+++ b/DiffractionGratings/grating3D_data_convergence.geo
@@ -21,7 +21,7 @@ DefineConstant[
     thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = {  4  , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = {  4  , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {107  , Name StrCat[pp3,"/1rx"]},
     ry            = {47   , Name StrCat[pp3,"/2ry"]},
     rz            = {40   , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_half_ellipsoid.geo b/DiffractionGratings/grating3D_data_half_ellipsoid.geo
index 6ae275a..07906f9 100644
--- a/DiffractionGratings/grating3D_data_half_ellipsoid.geo
+++ b/DiffractionGratings/grating3D_data_half_ellipsoid.geo
@@ -21,7 +21,7 @@ DefineConstant[
     thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = {  4      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = {  4      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {107  , Name StrCat[pp3,"/1rx"]},
     ry            = {47   , Name StrCat[pp3,"/2ry"]},
     rz            = {40   , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_hole.geo b/DiffractionGratings/grating3D_data_hole.geo
index 32bf712..1298f1a 100644
--- a/DiffractionGratings/grating3D_data_hole.geo
+++ b/DiffractionGratings/grating3D_data_hole.geo
@@ -27,7 +27,7 @@ DefineConstant[
     thick_L_6     = {100   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = { 2      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = { 2      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {250   , Name StrCat[pp3,"/1rx"]},
     ry            = {47    , Name StrCat[pp3,"/2ry"]},
     rz            = {500   , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_nanowire_solarcell.geo b/DiffractionGratings/grating3D_data_nanowire_solarcell.geo
index ccfd3c9..3a218f2 100644
--- a/DiffractionGratings/grating3D_data_nanowire_solarcell.geo
+++ b/DiffractionGratings/grating3D_data_nanowire_solarcell.geo
@@ -22,7 +22,7 @@ DefineConstant[
     thick_L_6     = {25   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = { 2      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = { 2      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = { 50   , Name StrCat[pp3,"/1rx"]},
     ry            = {47    , Name StrCat[pp3,"/2ry"]},
     rz            = {400   , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_pyramid.geo b/DiffractionGratings/grating3D_data_pyramid.geo
index 579361b..342c2bf 100644
--- a/DiffractionGratings/grating3D_data_pyramid.geo
+++ b/DiffractionGratings/grating3D_data_pyramid.geo
@@ -21,7 +21,7 @@ DefineConstant[
     thick_L_6     = {50    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = {  1      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = {  1      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {107   , Name StrCat[pp3,"/1rx"]},
     ry            = {250   , Name StrCat[pp3,"/2ry"]},
     rz            = {250   , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_retrieve_2D_lamellar.geo b/DiffractionGratings/grating3D_data_retrieve_2D_lamellar.geo
index 83f95d7..9a915e1 100644
--- a/DiffractionGratings/grating3D_data_retrieve_2D_lamellar.geo
+++ b/DiffractionGratings/grating3D_data_retrieve_2D_lamellar.geo
@@ -21,7 +21,7 @@ DefineConstant[
     thick_L_6     = {100    , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    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"}},
+    tag_geom      = {  7      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = { 0     , Name StrCat[pp3,"/1rx"]},
     ry            = {500        , Name StrCat[pp3,"/2ry"]},
     rz            = {1000     , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_skewed_lattice.geo b/DiffractionGratings/grating3D_data_skewed_lattice.geo
index 2bf61c7..47cfae8 100644
--- a/DiffractionGratings/grating3D_data_skewed_lattice.geo
+++ b/DiffractionGratings/grating3D_data_skewed_lattice.geo
@@ -22,7 +22,7 @@ DefineConstant[
     thick_L_6     = {50   , Name StrCat[pp2,"/8thickness layer 6 [nm] (substrate)"]},
     xsideg        = {30   , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
 
-    tag_geom      = {  4      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
+    tag_geom      = {  4      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="U",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
     rx            = {100 , Name StrCat[pp3,"/1rx"]},
     ry            = {100   , Name StrCat[pp3,"/2ry"]},
     rz            = {50    , Name StrCat[pp3,"/3rz"]},
diff --git a/DiffractionGratings/grating3D_data_torus.geo b/DiffractionGratings/grating3D_data_torus.geo
deleted file mode 100644
index bc27227..0000000
--- a/DiffractionGratings/grating3D_data_torus.geo
+++ /dev/null
@@ -1,125 +0,0 @@
-nm  = 1000;
-pp1 = "1Incident Plane Wave";
-pp2 = "2Layers Thicknesses";
-pp3 = "3Scatterer Properties";
-pp4 = "4Layer Materials";
-pp5 = "5Computational Parameters";
-pp6 = "6Output";
-DefineConstant[
-    FLAG_TOTAL    = {0     , Name StrCat[pp1,"/0Formulation"],Choices {0="scattered field",1="total field"}},
-    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)"]},
-    xsideg        = {0    , Name StrCat[pp2,"/9skew angle [deg]"],Visible 1},
-
-    tag_geom      = {  3      , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid",5="Checkerboard",6="bi-sinusoidal",7="2D lamellar"}},
-    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"]},
-    
-    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    = {12      , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
-    lc_scat       = {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_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= {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}, Visible 0},
-    PML_TYPE       = {0      , Name StrCat[pp5,"/9PML damping profile [-]"], Choices {0="constant profile",1="Bermudez profile"}, Visible 0},
-
-    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} },
-    FlagOutHtotCuts    = { 0    , Name StrCat[pp6,"/3Output Total Magnetic Field cuts?"] , Choices {0,1} },
-    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/4Output Scattered Electric Field cuts?"] , Choices {0,1} },
-    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/5Output Poynting cuts?"] , Choices {0,1} },
-    FlagOutEtotFull    = { 0    , Name StrCat[pp6,"/6Total Electric Field Full Output?"] , Choices {0,1} },
-    FlagOutEscaFull    = { 0    , Name StrCat[pp6,"/7Scattered Electric Field Full Output?"] , Choices {0,1} },
-    FlagOutPoyFull     = { 0    , Name StrCat[pp6,"/8Poynting Full Output?"] , Choices {0,1} }
-];
-
-lambda0       = nm * lambda0;
-period_x      = nm * period_x;
-period_y      = nm * period_y;
-thick_L_1     = nm * thick_L_1;
-thick_L_2     = nm * thick_L_2;
-thick_L_3     = nm * thick_L_3;
-thick_L_4     = nm * thick_L_4;
-thick_L_5     = nm * thick_L_5;
-thick_L_6     = nm * thick_L_6;
-rx            = nm * rx;
-ry            = nm * ry;
-rz            = nm * rz;
-lc_scat       = nm * lc_scat;
-PML_top       = nm * PML_top;
-PML_bot       = nm * PML_bot;
-InterpSampling= nm * InterpSampling;
-
-lambda_m = lambda0;
-og+=1;
-oi+=1;
-
-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;
-hh_L_7     = PML_bot_hh;
-hh_L_0     = PML_top_hh;
-thick_L_7  = PML_bot;
-thick_L_0  = PML_top;
-
-theta0 = thetadeg*Pi/180;
-phi0   = phideg*Pi/180;
-psi0   = psideg*Pi/180;
-xsi    = xsideg*Pi/180;
-
-dyc = period_y*Cos[xsi];
-dys = period_y*Sin[xsi];
-
-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_checkpoyX  = 50;
-npts_checkpoyY  = 50;
-npts_interpZSca = DomainZsizeSca/InterpSampling;
-npts_interpZTot = DomainZsizeTot/InterpSampling;
diff --git a/DiffractionGratings/grating3D_parallel_Mmatrix.sh b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
index 91287d8..9cbf398 100644
--- a/DiffractionGratings/grating3D_parallel_Mmatrix.sh
+++ b/DiffractionGratings/grating3D_parallel_Mmatrix.sh
@@ -8,7 +8,7 @@ export nb_phi=59
 export lambda_min=350
 export lambda_max=800
 export FLAG_TOTAL=0
-export GRATING_CASE="half_ellipsoid"
+export GRATING_CASE="U"
 export myDir="res_Matrix_nb_lam"$nb_lam"_nb_phi"$nb_phi"_total"$FLAG_TOTAL
 
 rm -r $myDir
@@ -23,7 +23,7 @@ myfunc() {
     local phi=$(echo "scale=10;360/($nb_phi)*$2" | bc )
     local psi=$(echo "scale=10;90*$3" | bc )
     cd $mysubDir
-    getdp grating3D.pro -pre helmholtz_vector -msh grating3D.msh -cal -pos postop_helmholtz_vector -petsc_prealloc 200 -setstring test_case $GRATING_CASE -setnumber lambda0 $lam -setnumber thetadeg 50 -setnumber phideg $phi -setnumber psideg $psi -setnumber FLAG_TOTAL $FLAG_TOTAL
+    getdp grating3D.pro -pre helmholtz_vector -msh grating3D.msh -cal -pos postop_helmholtz_vector -petsc_prealloc 200 -setstring test_case $GRATING_CASE -setnumber lambda0 $lam -setnumber thetadeg 15 -setnumber phideg $phi -setnumber psideg $psi -setnumber FLAG_TOTAL $FLAG_TOTAL
     if [ $3 -eq 0 ]; then cp res3D/rs.txt ../r_pin_sout_lam$1_phi$2.out ; fi
     if [ $3 -eq 0 ]; then cp res3D/rp.txt ../r_pin_pout_lam$1_phi$2.out ; fi
     if [ $3 -eq 1 ]; then cp res3D/rs.txt ../r_sin_sout_lam$1_phi$2.out ; fi
diff --git a/DiffractionGratings/grating3D_runall.sh b/DiffractionGratings/grating3D_runall.sh
index 9a90a02..b200ef0 100644
--- a/DiffractionGratings/grating3D_runall.sh
+++ b/DiffractionGratings/grating3D_runall.sh
@@ -1,6 +1,6 @@
 rm -r res3D*
 
-# for t in bi_sinusoidal checker half_ellipsoid hole pyramid torus retrieve_2D_lamellar nanowire_solarcell convergence skewed_lattice
+# for t in bi_sinusoidal checker half_ellipsoid hole pyramid U retrieve_2D_lamellar nanowire_solarcell convergence skewed_lattice
 for t in half_ellipsoid 
 do
     for f in 0 1
@@ -10,7 +10,7 @@ do
     done
 done
 
-# for t in bi_sinusoidal checker half_ellipsoid hole pyramid torus retrieve_2D_lamellar nanowire_solarcell convergence skewed_lattice
+# for t in bi_sinusoidal checker half_ellipsoid hole pyramid U retrieve_2D_lamellar nanowire_solarcell convergence skewed_lattice
 for t in half_ellipsoid  
 do
     for f in 0 1
diff --git a/ElectromagneticScattering/scattererTmatrix.geo b/ElectromagneticScattering/scattererTmatrix.geo
index 65b5a77..3b5cbce 100644
--- a/ElectromagneticScattering/scattererTmatrix.geo
+++ b/ElectromagneticScattering/scattererTmatrix.geo
@@ -30,7 +30,7 @@ If (flag_shape==CONE)
   Dilate { { 0,0,0 }, { 1 , cone_ry/cone_rx , 1} } { Volume{1}; }
 EndIf
 If (flag_shape==TOR)
-  Torus (1) = {0,0,0,tor_r1,tor_r2x,tor_angle*Pi/180};
+  U (1) = {0,0,0,tor_r1,tor_r2x,tor_angle*Pi/180};
   Dilate { { 0,0,0 }, { 1 , 1 , tor_r2z/tor_r2x} } { Volume{1}; }
 EndIf
 
diff --git a/ElectromagneticScattering/scattererTmatrix_data.geo b/ElectromagneticScattering/scattererTmatrix_data.geo
index 8699e89..764c01a 100644
--- a/ElectromagneticScattering/scattererTmatrix_data.geo
+++ b/ElectromagneticScattering/scattererTmatrix_data.geo
@@ -34,7 +34,7 @@ RES_QNM   = 3;
 
 DefineConstant[
   flag_shape = {ELL , Name StrCat[pp0, "0Scatterer shape"], 
-    Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split torus"},Closed 0}];
+    Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split U"},Closed 0}];
 
 If (flag_shape==ELL)
   DefineConstant[
@@ -88,10 +88,10 @@ If (flag_shape==CONE)
 EndIf
 If (flag_shape==TOR)
   DefineConstant[
-    tor_r1    = { 300 , Name StrCat[pp0 , "1torus radius 1  [nm]"], Highlight Str[colorppOK] , Closed 0},
-    tor_r2x   = { 100 , Name StrCat[pp0 , "2torus radius 2x [nm]"], Highlight Str[colorppOK] , Closed 0},
-    tor_r2z   = { 50  , Name StrCat[pp0 , "3torus radius 2z [nm]"], Highlight Str[colorppOK] , Closed 0},
-    tor_angle = { 340 , Name StrCat[pp0 , "4torus angle [deg]"]   , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
+    tor_r1    = { 300 , Name StrCat[pp0 , "1U radius 1  [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_r2x   = { 100 , Name StrCat[pp0 , "2U radius 2x [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_r2z   = { 50  , Name StrCat[pp0 , "3U radius 2z [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_angle = { 340 , Name StrCat[pp0 , "4U angle [deg]"]   , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
   rbb = tor_r1+tor_r2x;
 EndIf
 DefineConstant[
diff --git a/ElectromagneticScattering/scattering_data.geo b/ElectromagneticScattering/scattering_data.geo
index 8cea03e..ab086e3 100644
--- a/ElectromagneticScattering/scattering_data.geo
+++ b/ElectromagneticScattering/scattering_data.geo
@@ -34,7 +34,7 @@ RES_QNM   = 3;
 
 DefineConstant[
   flag_shape = {ELL , Name StrCat[pp0, "0Scatterer shape"], 
-    Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split torus"},Closed 0}];
+    Choices {ELL="ellispoid",PARALL="parallelepiped",CYL="cylinder",CONE="cone",TOR="split U"},Closed 0}];
 
 If (flag_shape==ELL)
   DefineConstant[
@@ -88,10 +88,10 @@ If (flag_shape==CONE)
 EndIf
 If (flag_shape==TOR)
   DefineConstant[
-    tor_r1    = { 300 , Name StrCat[pp0 , "1torus radius 1  [nm]"], Highlight Str[colorppOK] , Closed 0},
-    tor_r2x   = { 100 , Name StrCat[pp0 , "2torus radius 2x [nm]"], Highlight Str[colorppOK] , Closed 0},
-    tor_r2z   = { 50  , Name StrCat[pp0 , "3torus radius 2z [nm]"], Highlight Str[colorppOK] , Closed 0},
-    tor_angle = { 340 , Name StrCat[pp0 , "4torus angle [deg]"]   , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
+    tor_r1    = { 300 , Name StrCat[pp0 , "1U radius 1  [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_r2x   = { 100 , Name StrCat[pp0 , "2U radius 2x [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_r2z   = { 50  , Name StrCat[pp0 , "3U radius 2z [nm]"], Highlight Str[colorppOK] , Closed 0},
+    tor_angle = { 340 , Name StrCat[pp0 , "4U angle [deg]"]   , Highlight Str[colorppOK] , Closed 0, Min 5, Max 355}];
   rbb = tor_r1+tor_r2x;
 EndIf
 DefineConstant[
-- 
GitLab