Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • albertpiwonski-master-patch-57409
  • quadspheres
  • fix_Tmatrix_code_epsr_background
  • albertpiwonski-master-patch-12427
  • cavity
  • c1
7 results
Blame Permalink
grating3D.pro 20.04 KiB 
// test_case = "pyramid";
// test_case = "hole";
// test_case = "torus";
// test_case = "halfellipsoid";
// test_case = "checker";
// test_case = "bisin";
Include StrCat["grating3D_data_",test_case,".geo"];
Include "grating3D_materials.pro"

myDir = "run_results3D/";

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;
	Pinc        =  0.5*Ae^2*Sqrt[eps_re_L_1*ep0/mu0] * Cos[theta0];

  // permittivities
  For i In {1:6}
    If (flag_mat~{i}==0)
      epsr[L~{i}]  = Complex[eps_re_L~{i} , eps_im_L~{i}] * TensorDiag[1,1,1];
      If (i==1)
        epsr1[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
      EndIf
      If (i==6)
        epsr2[] = Complex[eps_re_L~{i} , eps_im_L~{i}];
      EndIf
    Else
      For j In {2:nb_available_materials}
        If(flag_mat~{i}==j)
          epsr[L~{i}]  = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
          If (i==1)
            epsr1[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
          EndIf
          If (i==6)
            epsr2[] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]];
          EndIf
        EndIf
      EndFor
    EndIf
  EndFor

  If (flag_mat_scat==0)
    epsr[Scat] = Complex[eps_re_Scat , eps_im_Scat] * TensorDiag[1,1,1];
  Else
    For j In {flag_mat_scat:flag_mat_scat}
      epsr[Scat] = Complex[epsr_re_interp_mat~{j}[lambda0/nm*1e-9] , epsr_im_interp_mat~{j}[lambda0/nm*1e-9]] * TensorDiag[1,1,1];
    EndFor
  EndIf

	For i In {1:6}
		mur[L~{i}]   = TensorDiag[1,1,1];
	EndFor
	mur[Scat]  = TensorDiag[1,1,1];

  ////// PMLS
	a_pml           = 1.;
	b_pml           = 1.;
  // bermu
  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]  = 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[]];
	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] = epsr1[] * 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[]]];

  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[]);

  spol[]    = Vector[Sin[phi0],-Cos[phi0],0];
  AmplEis[] =      spol[];
  AmplErs[] = rs[]*spol[];
  AmplEts[] = ts[]*spol[];
  AmplHis[] = Sqrt[ep0*epsr1[]/mu0]      *spol[];
  AmplHrs[] = Sqrt[ep0*epsr1[]/mu0]*rp[]*spol[];
  AmplHts[] = Sqrt[ep0*epsr1[]/mu0]*tp[]*spol[];

  Eis[]     = AmplEis[] * Exp[I[]*k1[] *XYZ[]];
  Ers[]     = AmplErs[] * Exp[I[]*k1r[]*XYZ[]];
  Ets[]     = AmplEts[] * Exp[I[]*k2[] *XYZ[]];
  His[]     = AmplHis[] * Exp[I[]*k1[] *XYZ[]];
  Hrs[]     = AmplHrs[] * Exp[I[]*k1r[]*XYZ[]];
  Hts[]     = AmplHts[] * Exp[I[]*k2[] *XYZ[]];
  Eip[]     = -1/(om0*ep0*epsr1[]) * k1[]  /\ His[];
  Erp[]     = -1/(om0*ep0*epsr1[]) * k1r[] /\ Hrs[];
  Etp[]     = -1/(om0*ep0*epsr2[]) * k2[]  /\ Hts[];

  Ei[] = Ae*(Cos[psi0]*Eip[]-Sin[psi0]*Eis[]);
  Er[] = Ae*(Cos[psi0]*Erp[]-Sin[psi0]*Ers[]);
  Et[] = Ae*(Cos[psi0]*Etp[]-Sin[psi0]*Ets[]);
  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*Nmax+1;
  For i In {0:Nb_ordre-1}
    alpha~{i}[]      = -CompX[k1[]] + 2*Pi/period_x * (i-Nmax);
    expialphax~{i}[] = Exp[I[]*alpha~{i}[]*X[]];
  EndFor
  For j In {0:Nb_ordre-1}
    beta~{j}[]      = -CompY[k1[]]  + 2*Pi/period_y * (j-Nmax);
    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 I1 ;
    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 I1; }
      Galerkin { [k0^2*epsr[] * Dof{u}      ,      {u}]; In Omega       ; Jacobian JVol; Integration I1; }
      Galerkin { [source[]                  ,      {u}]; In Omega_source; Jacobian JVol; Integration I1; }
    }
  }
}

Resolution {
  { Name helmholtz_vector;
    System {
      { Name M; NameOfFormulation helmholtz_vector; Type ComplexValue; }
    }
    Operation { 
      CreateDir[Str[myDir]];
      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 H1y     ; Value { Local { [CompY[H1[]]  ]; 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 I1 ; 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 I1 ; 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 I1 ; Jacobian JSur ; } } }
          { Name int_y_t~{i}~{j} ; Value { Integral { [ CompY[{u}+E1[] ]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntBot ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_x_r~{i}~{j} ; Value { Integral { [ CompX[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; Jacobian JSur ; } } }
          { Name int_y_r~{i}~{j} ; Value { Integral { [ CompY[{u}+E1d[]]*expialphax~{i}[]*expibetay~{j}[]/(period_x*period_y) ] ; In SurfIntTop ; Integration I1 ; 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 numbering_ij~{i}~{j}   ; Value { Term{Type Global; [Vector[i-Nmax,j-Nmax,0]] ; In SurfIntBot ; } } }
        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 (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 StrCat[myDir,"u_grid.pos"], Name "u_grid"];
        Else
          Print [ u , OnElementsOf Omega, File StrCat[myDir,"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 StrCat[myDir,"Etot_grid.pos"], Name "Etot_grid"];
        Else
          Print [ Etot , OnElementsOf Omega_plot, File StrCat[myDir,"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 StrCat[myDir,"Poy_tot_grid.pos"], Name "Poy_tot_grid"];
        Else
          Print [ Poy_tot , OnElementsOf Omega_plot, File StrCat[myDir,"Poy_tot.pos"]];
        EndIf
      EndIf
      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 StrCat[myDir,"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 StrCat[myDir,"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 StrCat[myDir,"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 StrCat[myDir,"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_interpZTot} , File StrCat[myDir,"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_interpZTot} , File StrCat[myDir,"Poy_tot_cut_X=0.pos"], Name "Poy_tot_cut_X=0"];
      EndIf
      // // For DEBUG
      Print [ H1y , OnElementsOf Omega, File StrCat[myDir,"H1y.pos"]];
      Print [ Etot , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
      // Print [ epsr_xx , OnElementsOf Omega, File StrCat[myDir,"epsr_xx.pos"]];
      // Print [ epsr_0xx , OnElementsOf Scat, File StrCat[myDir,"epsr_xx.pos"]];
      // Print [ Ecm , OnElementsOf Omega_plot, File StrCat[myDir,"Ecm.pos"]];
      // Print [ E1, OnElementsOf Omega_plot, File StrCat[myDir,"E1.pos"]];
      // Print [ Damp_pml_top, OnElementsOf PMLtop, File StrCat[myDir,"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 StrCat[myDir,> "./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 StrCat[myDir,> "./Views/Edif_XYcut.out",Format Table];
      // EndFor
      For k In {2:6}
        Print[ Abs_L~{k}[L~{k}], OnGlobal, File > StrCat[myDir,Sprintf["temp-Q_L_%g.txt",k]], Format Table ];
      EndFor
      Print[ Abs_scat[Scat]  , OnGlobal, File > StrCat[myDir,"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 > StrCat[myDir, "eff_t.txt"], Format Table ];
          Print[ eff_r~{i}~{j}[SurfIntTop], OnRegion SurfIntTop, File > StrCat[myDir, "eff_r.txt"], Format Table ];
        EndFor
      EndFor
      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 ];
        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}
];