grating3D initial commit

DiffractionGratings/grating3D.geo

+Include "grating3D_data_torus.geo";
+// Include "grating3D_data_hole.geo";
+// Include "grating3D_data_pyramid.geo";
+// Include "grating3D_data_halfellipsoid.geo";
+
+SetFactory("OpenCASCADE");
+// Geometry.LineWidth = 5;
+// Geometry.Points = 0;
+// Geometry.Color.Lines = {255,66,55};
+
+// General.BackgroundImageFileName="struct2.png";
+
+Macro SetPBCs
+    e = 1e-6;
+    masterX() = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e,-period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+    slaveX()  = Surface In BoundingBox{ period_x/2-e,-period_y/2-e, PML_bot_hh-e, period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+    masterY() = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e, period_x/2+e,-period_y/2+e, PML_top_hh+PML_top+e};
+    slaveY()  = Surface In BoundingBox{-period_x/2-e, period_y/2-e, PML_bot_hh-e, period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+    Periodic Surface{slaveX()} = {masterX()} Translate{period_x,        0,        0};
+    Periodic Surface{slaveY()} = {masterY()} Translate{       0, period_y,        0};
+Return
+
+
+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);
+Scat_n   = 2 ;
+
+PML_lc     = lambda_m/(paramaille*.5);
+L1_lc      = lambda_m/(paramaille*L1_n);
+L2_lc      = lambda_m/(paramaille*L2_n);
+L3_lc      = lambda_m/(paramaille*L3_n);
+L4_lc      = lambda_m/(paramaille*L4_n);
+L5_lc      = lambda_m/(paramaille*L5_n);
+L6_lc      = lambda_m/(paramaille*L6_n);
+
+Box(1) = {-period_x/2,-period_y/2, PML_bot_hh,period_x,period_y, PML_bot  };
+Box(2) = {-period_x/2,-period_y/2,     hh_L_6,period_x,period_y, thick_L_6};
+Box(3) = {-period_x/2,-period_y/2,     hh_L_5,period_x,period_y, thick_L_5};
+Box(4) = {-period_x/2,-period_y/2,     hh_L_4,period_x,period_y, thick_L_4};
+Box(5) = {-period_x/2,-period_y/2,     hh_L_3,period_x,period_y, thick_L_3};
+Box(6) = {-period_x/2,-period_y/2,     hh_L_2,period_x,period_y, thick_L_2};
+Box(7) = {-period_x/2,-period_y/2,     hh_L_1,period_x,period_y, thick_L_1};
+Box(8) = {-period_x/2,-period_y/2, PML_top_hh,period_x,period_y, PML_top  };
+w_arm1 = 35*nm;
+w_arm2 = 35*nm;
+w_arm3 = 35*nm;
+
+// tag_geom = 1; // half ellipsoid
+// tag_geom = 1; // pyramid
+
+If (tag_geom==1)
+    p=newp;
+    Point(p) = {0,0,hh_L_3+rz};
+    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};
+    Plane Surface(51) = {51};
+    Curve Loop(52) = {99, -46, 97};
+    Plane Surface(52) = {52};
+    Surface Loop(9) = {52, 51, 50, 49, 24};
+    Volume(9) = {9};
+EndIf
+
+If (tag_geom==2)
+    Cylinder(9) = {0,0,hh_L_3,0,0,thick_L_3,rx};
+EndIf
+
+If (tag_geom==3)
+    Torus(9) = {0,0,hh_L_3+ry-2,rx,ry};
+    Dilate { { 0,0,hh_L_3 }, { 1, 1, rz/ry } } { Volume{9}; }
+    BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+EndIf
+
+If (tag_geom==4)
+    Sphere(9) = {0,0,hh_L_3,rx};
+    Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, rz/rx } } { Volume{9}; }
+    BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+EndIf
+
+
+Coherence;
+
+// // // Physical Surface("ENDPML",770) = {69,37};
+
+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_L1_SUPER",7) = {7};
+Physical Volume("PMLTOP"        ,8) = {8};
+Physical Volume("SCAT"          ,9) = {9};
+
+Physical Surface("BXM" ,101 ) = masterX();
+Physical Surface("BXP" ,102 ) = slaveX();
+Physical Surface("BYM" ,201 ) = masterY();
+Physical Surface("BYP" ,202 ) = slaveY();
+Physical Surface("STOP",301 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_top_hh-e,period_x/2+e, period_y/2+e, PML_top_hh+e};
+Physical Surface("SBOT",302 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, hh_L_6-e,period_x/2+e, period_y/2+e, hh_L_6+e};
+Physical Surface("SPMLTOP",401 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_top_hh+PML_top-e,period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+Physical Surface("SPMLBOT",402 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e,period_x/2+e, period_y/2+e, PML_bot_hh+e};
+
+
+// Physical Volume("DEPOSIT"       ,10000) = {11};
+
+pts_PMLBOT()   = PointsOf{ Physical Volume{1}; };
+pts_LAYER_L6() = PointsOf{ Physical Volume{2}; };
+pts_LAYER_L5() = PointsOf{ Physical Volume{3}; };
+pts_LAYER_L4() = PointsOf{ Physical Volume{4}; };
+pts_LAYER_L3() = PointsOf{ Physical Volume{5}; };
+pts_LAYER_L2() = PointsOf{ Physical Volume{6}; };
+pts_LAYER_L1() = PointsOf{ Physical Volume{7}; };
+pts_PMLTOP()   = PointsOf{ Physical Volume{8}; };
+pts_ROD()      = PointsOf{ Physical Volume{9}; };
+// pts_DEPOSIT()  = PointsOf{ Physical Volume{10000}; };
+// For k In {0:#pts_PMLBOT()-1}
+//     Printf("blabal %g %f PML_lc;",pts_PMLBOT(k),PML_lc);
+// EndFor
+
+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;
+If (tag_geom==3)
+    Mesh.Algorithm = 6;
+EndIf
+Mesh.Lines = 0;

DiffractionGratings/grating3D.pro

+Include "grating3D_data_torus.geo";
+// Include "grating3D_data_hole.geo";
+// Include "grating3D_data_pyramid.geo";
+// Include "grating3D_data_halfellipsoid.geo";
+
+Include "gratings_common_materials.pro"
+Group {
+
+	// SubDomains
+	PMLbot   = Region[1];
+	L_6_temp = Region[2];
+	L_5      = Region[3];
+	L_4      = Region[4];
+	L_3      = Region[5];
+	L_2      = Region[6];
+	L_1_temp = Region[7];
+	PMLtop   = Region[8];
+	Scat     = Region[9];
+
+	// Boundaries
+	SurfBlochXm    = Region[101];
+	SurfBlochXp    = Region[102];
+	SurfBlochYm    = Region[201];
+	SurfBlochYp    = Region[202];
+  SurfIntTop     = Region[301];
+  SurfIntBot     = Region[302];
+
+  
+	SurfDirichlet  = Region[{401,402}];
+  SurfBloch      = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
+
+  L_1 = Region[{L_1_temp,SurfIntTop}];
+  L_6 = Region[{L_6_temp,SurfIntBot}];
+
+	Omega          = Region[{PMLbot,L_6,L_5,L_4,L_3,L_2,L_1,PMLtop,Scat}];
+	Omega_nosource = Region[{PMLbot,L_6,L_1,PMLtop}];
+	Omega_source   = Region[{Scat,L_2,L_3,L_4,L_5}];
+	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}];
+}
+
+
+
+Function{
+  I[] = Complex[0,1];
+  zhat[] = Vector[0,0,1];
+
+  small_delta = 0.0*nm;
+	mu0         = 4*Pi*1.e2*nm;
+	ep0         = 8.854187817e-3*nm;
+	cel         = 1/Sqrt[ep0 * mu0];
+	om0         = 2*Pi*cel/lambda0;
+	k0          = 2.*Pi/lambda0;
+	Ae          = 1;
+	Ah          = Ae*Sqrt[ep0/mu0];
+	Pinc        =  0.5*Ae^2*Sqrt[eps_re_L_1*ep0/mu0] * Cos[theta0];
+
+	For i In {1:6}
+		epsr[L~{i}]  = Complex[eps_re_L~{i} , eps_im_L~{i}] * TensorDiag[1,1,1];
+		mur[L~{i}]   = TensorDiag[1,1,1];
+	EndFor
+  epsr[Scat] = Complex[eps_re_Scat , eps_im_Scat] * TensorDiag[1,1,1];
+	mur[Scat]  = TensorDiag[1,1,1];
+
+  ////// PMLS
+	a_pml           = 1.;
+	b_pml           = 1.;
+  // bermu
+  epsr1[] = Complex[eps_re_L_1 , eps_im_L_1];
+ 	epsr2[] = Complex[eps_re_L_6 , eps_im_L_6];
+  n1[]    = Sqrt[epsr1[]];
+  n2[]    = Sqrt[epsr2[]];
+  k1norm[]= k0*n1[];
+  k2norm[]= k0*n2[];
+  Zmax = PML_top_hh;
+  Zmin = hh_L_6;
+  Damp_pml_top[] = 1/Abs[Zmax + PML_top - Z[]];
+  Damp_pml_bot[] = 1/Abs[Zmin - PML_bot - Z[]];
+  // DampingProfileX[] = 1/(Xmax + SizePMLX - Fabs[X[]]) - 1/(SizePMLX);
+
+  sx            = 1.;
+	sz[]          = Complex[a_pml,b_pml];
+	sz_bermutop[] = Complex[1,Damp_pml_top[]/k1norm[]];
+	sz_bermubot[] = Complex[1,Damp_pml_bot[]/k2norm[]];
+	sy            = 1.;
+  
+	epsr[PMLtop]  = eps_re_L_1*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	mur[PMLtop]   =            TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	epsr[PMLbot]  = eps_re_L_6*TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	mur[PMLbot]   =            TensorDiag[sz[]*sy/sx,sx*sz[]/sy,sx*sy/sz[]];
+	// 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[]];
+	// mur[PMLbot]      = TensorDiag[sz_bermubot[]*sy/sx,sx*sz_bermubot[]/sy,sx*sy/sz_bermubot[]];
+	
+	epsr_annex[PMLbot]       = epsr[];
+	epsr_annex[PMLtop]       = epsr[];
+	epsr_annex[Omega_source] = Complex[eps_re_L_1 , eps_im_L_1] * TensorDiag[1,1,1];
+	epsr_annex[L_1]          = epsr[];
+	epsr_annex[L_6]          = epsr[];
+
+  //// Reference Field solution of annex problem (simple diopter)  
+  k1[]    = k0*n1[]*Vector[-Sin[theta0]*Cos[phi0], -Sin[theta0]*Sin[phi0], -Cos[theta0]];
+  k2[] = Vector[
+            CompX[k1[]],
+            CompY[k1[]],
+           -Sqrt[k0^2*epsr2[]-CompX[k1[]]^2-CompY[k1[]]^2]
+        ];
+  k1r[]  = Vector[CompX[k1[]], CompY[k1[]], -CompZ[k1[]]];
+  k1s[]  = (Abs[theta0]>1e-12) ?  k1[]/\zhat[] / Norm[ k1[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+  k2s[]  = (Abs[theta0]>1e-12) ?  k2[]/\zhat[] / Norm[ k2[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+  k1rs[] = (Abs[theta0]>1e-12) ? k1r[]/\zhat[] / Norm[k1r[]/\zhat[]] : Vector[-Sin[phi0], Cos[phi0], 0.] ;
+  k1p[]  = (Abs[theta0]>1e-12) ?  k1[]/\k1s[]  / Norm[ k1[]/\k1s[]]  : Vector[ Cos[phi0], Sin[phi0], 0.] ;
+  k2p[]  = (Abs[theta0]>1e-12) ?  k2[]/\k2s[]  / Norm[ k2[]/\k2s[]]  : Vector[ Cos[phi0], Sin[phi0], 0.] ;
+  k1rp[] = (Abs[theta0]>1e-12) ? k1r[]/\k1rs[] / Norm[k1r[]/\k1rs[]] : Vector[-Cos[phi0],-Sin[phi0], 0.] ;
+
+  rs[] = (CompZ[k1[]]-CompZ[k2[]])/(CompZ[k1[]]+CompZ[k2[]]);
+  ts[] =           2.*CompZ[k1[]] /(CompZ[k1[]]+CompZ[k2[]]);
+  rp[] = (CompZ[k1[]]*epsr2[]-CompZ[k2[]]*epsr1[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]);
+  tp[] =                  (2.*CompZ[k1[]]*epsr2[])/(CompZ[k1[]]*epsr2[]+CompZ[k2[]]*epsr1[]) * Sqrt[epsr1[]/epsr2[]];
+
+  AmplEi[] =      Cos[psi0]* k1p[] +      Sin[psi0]* k1s[];
+  AmplEr[] = rp[]*Cos[psi0]*k1rp[] + rs[]*Sin[psi0]*k1rs[];
+  AmplEt[] = tp[]*Cos[psi0]* k2p[] + ts[]*Sin[psi0]* k2s[];
+  
+  Ei[]     = AmplEi[] * Exp[I[]*k1[]*XYZ[]];
+  Et[]     = AmplEt[] * Exp[I[]*k2[]*XYZ[]];
+  Er[]     = AmplEr[] * Exp[I[]*k1r[]*XYZ[]];
+  Hi[]     = 1/(om0*mu0*mur[]) * k1[]  /\ Ei[];
+  Hr[]     = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
+  Ht[]     = 1/(om0*mu0*mur[]) * k2[]  /\ Et[];
+
+  E1[Omega_super]  = Ei[]+Er[];
+  E1[Omega_subs]   = Et[];
+  E1d[Omega_super] = Er[];
+  E1d[Omega_subs]  = Et[];
+  H1[Omega_super]  = Hi[]+Hr[];
+  H1[Omega_subs]   = Ht[];
+  H1d[Omega_super] = Hr[];
+  H1d[Omega_subs]  = Ht[];
+
+  source[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
+  
+  // Bloch phase shifts
+  dephX[] = Exp[I[]*CompX[k1[]]*period_x];
+	dephY[] = Exp[I[]*CompY[k1[]]*period_y];
+
+	// Fourier coefficients variables
+  Nb_ordre = 2*N_d_order+1;
+  For i In {0:Nb_ordre-1}
+    alpha~{i}[]      = -CompX[k1[]] + 2*Pi/period_x * (i-N_d_order);
+    expialphax~{i}[] = Exp[I[]*alpha~{i}[]*X[]];
+  EndFor
+  For j In {0:Nb_ordre-1}
+    beta~{j}[]      = -CompY[k1[]]  + 2*Pi/period_y * (j-N_d_order);
+    expibetay~{j}[] = Exp[I[]*beta~{j}[]*Y[]];
+  EndFor
+  For i In {0:Nb_ordre-1}
+    For j In {0:Nb_ordre-1}
+      gammar~{i}~{j}[] = Sqrt[k0^2*epsr1[] - alpha~{i}[]^2 - beta~{j}[]^2];
+      gammat~{i}~{j}[] = Sqrt[k0^2*epsr2[] - alpha~{i}[]^2 - beta~{j}[]^2];
+    EndFor
+  EndFor
+
+}
+
+Constraint {
+  { Name Dirichlet; Type Assign;
+    Case {
+      { Region SurfDirichlet; Value 0.; }
+    }
+  }
+  { Name BlochX;
+    Case {
+    { Region SurfBlochXp; Type LinkCplx ; RegionRef SurfBlochXm;
+      Coefficient dephX[]; Function Vector[$X-period_x,$Y,$Z] ; }
+    }
+  }
+  { Name BlochY;
+    Case {
+      { Region SurfBlochYp; Type LinkCplx ; RegionRef SurfBlochYm;
+      Coefficient dephY[]; Function Vector[$X,$Y-period_y,$Z] ; }
+    }
+  }
+}
+
+Jacobian {
+  { Name JVol ;
+    Case { 
+      { Region All ; Jacobian Vol ; }
+    }
+  }
+  { Name JSur ;
+    Case {
+      { Region All ; Jacobian Sur ; }
+    }
+  }
+  { Name JLin ;
+    Case {
+      { Region All ; Jacobian Lin ; }
+    }
+  }
+}
+
+Integration {
+  { Name Int_1 ;
+    Case { 
+      { Type Gauss ;
+        Case { 
+          { GeoElement Point       ; NumberOfPoints   1 ; }
+          { GeoElement Line        ; NumberOfPoints   4 ; }
+          { GeoElement Triangle    ; NumberOfPoints  12 ; }
+          { GeoElement Triangle2   ; NumberOfPoints  12 ; }
+          { GeoElement Tetrahedron ; NumberOfPoints  16 ; }
+          { GeoElement Tetrahedron2; NumberOfPoints  16 ; }
+        }
+      }
+    }
+  }
+}
+
+FunctionSpace {
+  { Name Hcurl; Type Form1;
+    BasisFunction {
+      { Name sn; NameOfCoef un; Function BF_Edge; Support Region[{Omega}]; Entity EdgesOf[All]; }
+      { Name sn2; NameOfCoef un2; Function BF_Edge_2E;Support Region[{Omega}]; Entity EdgesOf[All]; }
+      { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b; Support Region[Omega]; Entity FacetsOf[Omega, Not SurfBloch]; }
+      { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c; Support Region[Omega]; Entity FacetsOf[Omega, Not SurfBloch]; }
+      { Name sn5; NameOfCoef un5; Function BF_Edge_4E  ; Support Region[Omega]; Entity EdgesOf[Omega, Not SurfBloch]; }
+    }
+    Constraint {
+      { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochX; }
+      { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint BlochY; }
+      { NameOfCoef un;  EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochX; }
+      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint BlochY; }
+      { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+     // { 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; }
+      }
+  }
+}
+
+Formulation {
+  { Name helmholtz_vector; Type FemEquation;
+    Quantity {
+      { Name u; Type Local; NameOfSpace Hcurl;}
+    }
+		Equation{
+      Galerkin { [-1/mur[]    * Dof{Curl u} , {Curl u}]; In Omega       ; Jacobian JVol; Integration Int_1; }
+      Galerkin { [k0^2*epsr[] * Dof{u}      ,      {u}]; In Omega       ; Jacobian JVol; Integration Int_1; }
+      Galerkin { [source[]                  ,      {u}]; In Omega_source; Jacobian JVol; Integration Int_1; }
+    }
+  }
+}
+
+Resolution {
+  { Name helmholtz_vector;
+    System {
+      { Name M; NameOfFormulation helmholtz_vector; Type ComplexValue; }
+    }
+    Operation { 
+      Generate[M]; Solve[M]; //SaveSolution[M];  
+    }
+  }
+}
+
+PostProcessing {   
+  { Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
+    Quantity {
+      { Name u      ; Value { Local { [ {u}       ]; In Omega; Jacobian JVol; } } }
+      { 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 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; } } }
+      { Name Damp_pml_top; Value { Local { [Damp_pml_top[]  ]; In Omega; Jacobian JVol; } } }
+      { Name Poy_tot; Value { Local { [  0.5*Re[Cross[{u}+E1[] , Conj[ H1[]-I[]/(mur[]*mu0*om0)*{Curl u}]]]  ]; In Omega; Jacobian JVol; } } }
+      
+      For k In {2:6}
+        { Name Abs_L~{k} ; Value { Integral { [ ep0*om0 * 0.5*Fabs[eps_im_L~{k}]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In L~{k} ; Integration Int_1 ; Jacobian JVol ; } } }
+      EndFor
+      { Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Fabs[eps_im_Scat]*(SquNorm[{u}+E1[]]) / (Pinc*period_x*period_y) ] ; In Scat ; Integration Int_1 ; Jacobian JVol ; } } }
+
+      For i In {0:Nb_ordre-1}
+        For j In {0:Nb_ordre-1}
+          { Name int_x_t~{i}~{j} ; Value { Integral { [ CompX[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration Int_1 ; Jacobian JSur ; } } }
+          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration Int_1 ; Jacobian JSur ; } } }
+          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration Int_1 ; Jacobian JSur ; } } }
+          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration Int_1 ; Jacobian JSur ; } } }
+          { Name eff_t~{i}~{j}   ; Value { Term{Type Global; [
+                 1/(gammat~{i}~{j}[]*-CompZ[k1[]]) *((gammat~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_t~{i}~{j}]+
+                                                       (gammat~{i}~{j}[]^2+ beta~{j}[]^2)*SquNorm[$int_y_t~{i}~{j}]+
+                                                       2*alpha~{i}[]*beta~{j}[]*Re[$int_x_t~{i}~{j}*Conj[$int_y_t~{i}~{j}]]) ] ; In SurfIntBot ; } } }
+          { Name eff_r~{i}~{j}   ; Value { Term{Type Global; [
+                 1/(gammar~{i}~{j}[]*-CompZ[k1[]])*((gammar~{i}~{j}[]^2+alpha~{i}[]^2)*SquNorm[$int_x_r~{i}~{j}]+
+                                                      (gammar~{i}~{j}[]^2+ beta~{j}[]^2)*SquNorm[$int_y_r~{i}~{j}]+
+                                                      2*alpha~{i}[]*beta~{j}[]*Re[$int_x_r~{i}~{j}*Conj[$int_y_r~{i}~{j}]]) ] ; In SurfIntTop ; } } }
+          { Name gammat~{i}~{j}   ; Value { Term{Type Global; [gammat~{i}~{j}[]] ; In SurfIntBot ; } } }
+          { Name gammar~{i}~{j}   ; Value { Term{Type Global; [gammar~{i}~{j}[]] ; In SurfIntTop ; } } }          
+        EndFor
+      EndFor
+
+      For i In {0:Nb_ordre-1}
+        { Name alpha~{i}       ; Value { Term{Type Global; [alpha~{i}[]] ; In SurfIntBot ; } } }
+      EndFor
+      For j In {0:Nb_ordre-1}
+        { Name beta~{j}   ; Value { Term{Type Global; [beta~{j}[]] ; In SurfIntBot ; } } }
+      EndFor
+    }
+  }
+}
+
+
+
+
+PostOperation {
+  { Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
+    Operation {
+      If (FlagOutEscaCuts==1)
+        Print [ u , OnPlane { {-period_x/2,0,hh_L_6-PML_bot} {period_x/2,0,hh_L_6-PML_bot} {-period_x/2,0,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpZSca} , File "u_cut_Y=0.pos", Name "u_cut_Y=0"];
+        Print [ u , OnPlane { {0,-period_y/2,hh_L_6-PML_bot} {0,period_y/2,hh_L_6-PML_bot} {0,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpY,npts_interpZSca} , File "u_cut_X=0.pos", Name "u_cut_X=0"];
+      EndIf
+      If (FlagOutEtotCuts==1)
+        Print [ Etot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZTot} , File "Etot_cut_Y=0.pos", Name "Etot_cut_Y=0"];
+        Print [ Etot , 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 "Etot_cut_X=0.pos", Name "Etot_cut_X=0"];
+      EndIf
+      If (FlagOutPoyCut==1)
+        Print [ Poy_tot , OnPlane { {-period_x/2,0,hh_L_6} {period_x/2,0,hh_L_6} {-period_x/2,0,hh_L_1+thick_L_1} } {npts_interpX,npts_interpZSca} , File "Poy_tot_cut_Y=0.pos", Name "Poy_tot_cut_Y=0"];
+        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_interpZSca} , File "Poy_tot_cut_X=0.pos", Name "Poy_tot_cut_X=0"];
+      EndIf
+      If (FlagOutEscaFull==1)
+        If (Flag_interp_cubic==1)
+          Print [ u , OnBox { {-period_x/2,-period_y/2,hh_L_6-PML_bot} {period_x/2,-period_y/2,hh_L_6-PML_bot} {-period_x/2,period_y/2,hh_L_6-PML_bot} {-period_x/2,-period_y/2,hh_L_1+thick_L_1+PML_top} } {npts_interpX,npts_interpY,npts_interpZSca} , File "u_grid.pos", Name "u_grid"];
+        Else
+          Print [ u , OnElementsOf Omega, File "Etot.pos"];
+        EndIf
+      EndIf
+      If (FlagOutEtotFull==1)
+        If (Flag_interp_cubic==1)
+          Print [ Etot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File "Etot_grid.pos", Name "Etot_grid"];
+        Else
+          Print [ Etot , OnElementsOf Omega_plot, File "Etot.pos"];
+        EndIf
+      EndIf
+      If (FlagOutPoyFull==1)
+        If (Flag_interp_cubic==1)
+          Print [ Poy_tot , OnBox { {-period_x/2,-period_y/2,hh_L_6} {period_x/2,-period_y/2,hh_L_6} {-period_x/2,period_y/2,hh_L_6} {-period_x/2,-period_y/2,hh_L_1+thick_L_1} } {npts_interpX,npts_interpY,npts_interpZTot} , File "Poy_tot_grid.pos", Name "Poy_tot_grid"];
+        Else
+          Print [ Poy_tot , OnElementsOf Omega_plot, File "Poy_tot.pos"];
+        EndIf
+      EndIf
+      // // For DEBUG
+      // Print [ epsr_xx , OnElementsOf Omega, File "epsr_xx.pos"];
+      // Print [ Ecm , OnElementsOf Omega_plot, File "Ecm.pos"];
+      // Print [ E1, OnElementsOf Omega_plot, File "E1.pos"];
+      // Print [ Damp_pml_top, OnElementsOf PMLtop, File "Damp_pml_top.pos"];
+      // For i In {0:nb_slice-1}
+      //   Print [ Etot , OnPlane { {-period_x/2+small_delta,-period_y/2+small_delta,zcut_sub~{i}} {-period_x/2+small_delta,period_y/2-small_delta,zcut_sub~{i}} {period_x/2-small_delta,-period_y/2+small_delta,zcut_sub~{i}} } {ninterv_integ,ninterv_integ} , File > "./Views/Etot_XYcut.out",Format Table];
+      //   Print [ Edif , OnPlane { {-period_x/2+small_delta,-period_y/2+small_delta,zcut_sup~{i}} {-period_x/2+small_delta,period_y/2-small_delta,zcut_sup~{i}} {period_x/2-small_delta,-period_y/2+small_delta,zcut_sup~{i}} } {ninterv_integ,ninterv_integ} , File > "./Views/Edif_XYcut.out",Format Table];
+      // EndFor
+      For k In {2:6}
+        Print[ Abs_L~{k}[L~{k}], OnGlobal, File Sprintf["temp-Q_L_%g.txt",k], Format Table ];
+      EndFor
+      Print[ Abs_scat[Scat]  , OnGlobal, File "temp-Q_scat.txt", Format Table ];
+
+      For i In {0:Nb_ordre-1}
+        For j In {0:Nb_ordre-1}
+          Print[ int_x_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_x_t~{i}~{j}, Format Table];
+          Print[ int_y_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_y_t~{i}~{j}, Format Table];
+          Print[ int_x_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_x_r~{i}~{j}, Format Table];
+          Print[ int_y_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_y_r~{i}~{j}, Format Table];
+        EndFor
+      EndFor
+
+      For i In {0:Nb_ordre-1}
+        For j In {0:Nb_ordre-1}
+          Print[ eff_t~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > "eff_t.txt", Format Table ];
+          Print[ eff_r~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File > "eff_r.txt", Format Table ];
+        EndFor
+      EndFor
+      For i In {0:Nb_ordre-1}
+        For j In {0:Nb_ordre-1}
+          Print[ gammat~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > "gammat.txt", Format Table ];
+          Print[ gammar~{i}~{j}[SurfIntBot], OnRegion SurfIntTop, File > "gammar.txt", Format Table ];
+        EndFor
+      EndFor
+    }
+  }
+}
+
+DefineConstant[
+  R_ = {"helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+  C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+  P_ = {"postop_helmholtz_vector", Name "GetDP/2PostOperationChoices", Visible 1}
+];

DiffractionGratings/grating3D_data_torus.geo

+nm  = 1;
+pp1 = "1Incident Plane Wave";
+pp2 = "2Layers Thicknesses";
+pp3 = "3Scatterer Properties";
+pp4 = "4Layer Materials";
+pp5 = "5Computational Paramameters";
+pp6 = "6Output";
+DefineConstant[
+    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     = {200    , Name StrCat[pp2,"/4thickness layer 2 [nm]"]},
+    thick_L_3     = {100   , Name StrCat[pp2,"/5thickness layer 3 [nm]"]},
+    thick_L_4     = {200    , 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)"]},
+    
+    tag_geom      = {  3   , Name StrCat[pp3,"/0Shape"], Choices {1="Pyramid",2="Cylindrical Hole",3="Torus",4="HalfEllipspoid"}},
+    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 eps_re_Scat"]},
+    eps_im_Scat   = { 20   , Name StrCat[pp3,"/6Custom eps_im_Scat"]},
+
+    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"]},
+    
+    paramaille    = {10         , Name StrCat[pp5,"/1Number of mesh elements per wavelength [-]"]},
+    scat_lc       = {20         , 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]"]},
+    N_d_order     = {1          , Name StrCat[pp5,"/6Number of non specular order to output [-]"]},
+    refine_mesh_L1= {1          , Name StrCat[pp5,"/7refine layers/1refine mesh layer 1 [-]"]},
+    refine_mesh_L2= {1          , Name StrCat[pp5,"/7refine layers/2refine mesh layer 2 [-]"]},
+    refine_mesh_L3= {1          , Name StrCat[pp5,"/7refine layers/3refine mesh layer 3 [-]"]},
+    refine_mesh_L4= {1          , Name StrCat[pp5,"/7refine layers/4refine mesh layer 4 [-]"]},
+    refine_mesh_L5= {1          , Name StrCat[pp5,"/7refine layers/5refine mesh layer 5 [-]"]},
+    refine_mesh_L6= {1          , Name StrCat[pp5,"/7refine layers/6refine mesh layer 6 [-]"]},
+    
+    InterpSampling     = { 30   , Name StrCat[pp6,"/0Interpolation grid step [nm]"]},
+    Flag_interp_cubic  = { 0    , Name StrCat[pp6,"/1Interpolate on cubic grid?"], Choices {0,1} },
+    FlagOutEtotCuts    = { 1    , Name StrCat[pp6,"/2Output Total Electric Field cuts?"] , Choices {0,1} },
+    FlagOutEscaCuts    = { 1    , Name StrCat[pp6,"/3Output Scattered Electric Field cuts?"] , Choices {0,1} },
+    FlagOutPoyCut      = { 1    , Name StrCat[pp6,"/4Output Poynting cuts?"] , Choices {0,1} },
+    FlagOutEtotFull    = { 1    , Name StrCat[pp6,"/5Total Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutEscaFull    = { 1    , Name StrCat[pp6,"/6Scattered Electric Field Full Output?"] , Choices {0,1} },
+    FlagOutPoyFull     = { 1    , Name StrCat[pp6,"/7Poynting Full Output?"] , Choices {0,1} }
+];
+
+lambda_m = lambda0;
+
+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;
+
+theta0 = thetadeg*Pi/180;
+phi0 = phideg*Pi/180;
+psi0 = psideg*Pi/180;
+
+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_interpZSca = DomainZsizeSca/InterpSampling;
+npts_interpZTot = DomainZsizeTot/InterpSampling;