diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 947806f0281a04d99bc4ad2d4c1c2f475e6c7d75..73cf7ebbb85e3175efe6f3a375a1a9be696afc3d 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -1,21 +1,21 @@
// 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;
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index f3e06e87b754ba70be3c3d03ecab519192cb4e63..4dde9ddb017ce164afd2bf1c8f0e474baaa33eec 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -1,22 +1,22 @@
// 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
diff --git a/DiffractionGratings/grating3D_data.geo b/DiffractionGratings/grating3D_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..d08505222fd57b476b19ca9d653d915070b2957b
--- /dev/null
+++ b/DiffractionGratings/grating3D_data.geo
@@ -0,0 +1,10 @@
+
+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"];