make grating3D work interactively

DiffractionGratings/grating3D.geo

 // Copyright (C) 2019 Guillaume Demésy
 //
 // This file is part of the model grating3D.pro.
-// 
+//
 // This program is free software: you can redistribute it and/or modify
 // it under the terms of the GNU General Public License as published by
 // the Free Software Foundation, either version 3 of the License, or
 // (at your option) any later version.
-// 
+//
 // This program is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 // GNU General Public License for more details.
-// 
+//
 // You should have received a copy of the GNU General Public License
 // along with This program. If not, see <https://www.gnu.org/licenses/>.
 
-Include StrCat["grating3D_data_",test_case,".geo"];
+Include "grating3D_data.geo";
 
 SetFactory("OpenCASCADE");
 e = 5*nm;
@@ -119,7 +119,7 @@ If (tag_geom==6)
   Curve Loop(30) = {31, -23, 32, 4, 10};
   Surface(18) = {30};
   Surface Loop(3) = {5, 16, 18, 12, 11, 9, 8, 7, 10};
-  Volume(1) = {3};    
+  Volume(1) = {3};
   Surface Loop(2) = {6, 15, 14, 10, 7, 8, 9, 13, 17};
   Volume(2) = {2};
 EndIf
@@ -288,4 +288,4 @@ If (tag_geom==3) // Split torus weird otherwise
 EndIf
 // Mesh.SurfaceEdges = 0;
 Mesh.VolumeEdges = 0;
-Mesh.ElementOrder = og;
\ No newline at end of file
+Mesh.ElementOrder = og;

DiffractionGratings/grating3D.pro

 // Copyright (C) 2019 Guillaume Demésy
 //
 // This file is part of the model grating3D.pro.
-// 
+//
 // This program is free software: you can redistribute it and/or modify
 // it under the terms of the GNU General Public License as published by
 // the Free Software Foundation, either version 3 of the License, or
 // (at your option) any later version.
-// 
+//
 // This program is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 // GNU General Public License for more details.
-// 
+//
 // You should have received a copy of the GNU General Public License
 // along with This program. If not, see <https://www.gnu.org/licenses/>.
 
-Include StrCat["grating3D_data_",test_case,".geo"];
-Include "grating3D_materials.pro"
+Include "grating3D_data.geo";
+Include "grating3D_materials.pro";
 
 myDir = "res3D/";
 
@@ -40,7 +40,7 @@ Group {
   SurfIntTop     = Region[301];
   SurfIntBot     = Region[302];
 
-  
+
 	SurfDirichlet  = Region[{401,402}];
   SurfBloch      = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
 
@@ -49,7 +49,7 @@ Group {
   Else
     SurfExcludeFacets  = Region[{SurfBloch}];
   EndIf
-  
+
   L_1 = Region[{L_1_temp,SurfIntTop}];
   L_6 = Region[{L_6_temp,SurfIntBot}];
 
@@ -59,7 +59,7 @@ Group {
 	Omega_super    = Region[{Scat,L_1,L_2,L_3,L_4,L_5,PMLtop}];
 	Omega_subs     = Region[{L_6,PMLbot}];
 	Omega_plot     = Region[{L_6,L_5,L_4,L_3,L_2,L_1,Scat}];
-	
+
 	// SurfNeumann    = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
 }
 
@@ -135,7 +135,7 @@ Function{
 	sz_bermutop[] = Complex[1,Damp_pml_top[]/k1norm[]];
 	sz_bermubot[] = Complex[1,Damp_pml_bot[]/k2norm[]];
 	sy            = 1.;
-  
+
 	epsr[PMLtop]  = Re[epsr1[]]*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
 	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[]];
@@ -143,7 +143,7 @@ Function{
 
   // epsr[PMLbot] = epsr2[];
   // mur[PMLbot]  = TensorDiag[1,1,1];
-	
+
 	epsr_annex[PMLbot]       = epsr[];
 	epsr_annex[PMLtop]       = epsr[];
 	epsr_annex[Omega_source] = epsr1[] * TensorDiag[1,1,1];
@@ -201,7 +201,7 @@ Function{
   H1d[Omega_subs]  = Ht[];
 
   source[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
-  
+
   // Bloch phase shifts
   skx1[] =  k1x[];
   // sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
@@ -250,7 +250,7 @@ Constraint {
 
 Jacobian {
   { Name JVol ;
-    Case { 
+    Case {
       { Region All ; Jacobian Vol ; }
     }
   }
@@ -268,9 +268,9 @@ Jacobian {
 
 Integration {
   { Name I1 ;
-    Case { 
+    Case {
       { Type Gauss ;
-        Case { 
+        Case {
           { GeoElement Point       ; NumberOfPoints   1 ; }
           { GeoElement Line        ; NumberOfPoints   4 ; }
           { GeoElement Triangle    ; NumberOfPoints  12 ; }
@@ -334,18 +334,18 @@ Resolution {
     System {
       { Name M; NameOfFormulation helmholtz_vector; Type ComplexValue; }
     }
-    Operation { 
+    Operation {
       CreateDir[Str[myDir]];
-      Generate[M]; 
-      Solve[M]; //SaveSolution[M];  
+      Generate[M];
+      Solve[M]; //SaveSolution[M];
     }
   }
 }
 
-PostProcessing {   
+PostProcessing {
   { Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
     Quantity {
-      { Name u      ; Value { Local { [ {u}       ]; In Omega; Jacobian JVol; } } }      
+      { 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 Edif   ; Value { Local { [ {u}+E1d[] ]; In Omega; Jacobian JVol; } } }
@@ -511,4 +511,4 @@ DefineConstant[
 // related to? https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2018-October/036377.html
 
 // N.B. pastix works with petsc 3.12
-// -pc_factor_mat_solver_type pastix
\ No newline at end of file
+// -pc_factor_mat_solver_type pastix

DiffractionGratings/grating3D_data.geo

+
+DefineConstant[
+  test_case = {"halfellipsoid",
+    Name "1Geometry",
+    Choices {"hole", "pyramid", "torus", "2Dlamellar", "solarcell", "conv", "skew"},
+    GmshOption "Reset", Autocheck 0
+  }
+];
+
+Include StrCat["grating3D_data_",test_case,".geo"];