diff --git a/DiffractionGratings/grating2D_conical.pro b/DiffractionGratings/grating2D_conical.pro
index 0c9fada46a7f37a6a5a02ab43a512c5d139c4f51..4250dcbc4aebfc721471e2b8ae9fd19ad6f3d799 100644
--- a/DiffractionGratings/grating2D_conical.pro
+++ b/DiffractionGratings/grating2D_conical.pro
@@ -1,17 +1,17 @@
// Copyright (C) 2020 Guillaume Demésy
//
// This file is part of the model grating2D.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/>.
@@ -59,7 +59,7 @@ Group {
Plot_bnd = Region[SURF_PLOT];
Print_point = Region[PRINT_POINT];
- Omega_DefTr = Region[{SUB,SUP}];
+ Omega_DefTr = Region[{SUB,SUP}];
Gama = Region[{SurfCutSubs1,SurfCutSuper1}];
Tr = ElementsOf[Omega_DefTr, ConnectedTo Gama];
}
@@ -130,7 +130,7 @@ Function{
epsr_layer_cov_im[] = epsr_im_dom_5[];
epsr1_re[] = epsr_re_dom_6[];
epsr1_im[] = epsr_im_dom_6[];
-
+
om0 = 2.*Pi*cel/lambda0;
k0 = 2.*Pi/lambda0;
epsr1[] = Complex[epsr1_re[],epsr1_im[]];
@@ -191,7 +191,7 @@ Function{
a_pml = 1.;
b_pml = 1.;
sx = 1.;
- sy[] = Complex[a_pml,b_pml];
+ sy[] = Complex[a_pml,b_pml];
sz = 1.;
PML_Tensor[] = TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz];
@@ -215,8 +215,8 @@ Function{
EndIf
epsilonr[layer_cov] = Complex[epsr_layer_cov_re[],epsr_layer_cov_im[]] * TensorDiag[1,1,1];
epsilonr[pmltop] = epsr1_re[]*PML_Tensor[];
- epsilonr[pmlbot] = epsr2_re[]*PML_Tensor[];
-
+ epsilonr[pmlbot] = epsr2_re[]*PML_Tensor[];
+
epsilonr_annex[sub] = epsr2[] * TensorDiag[1,1,1];
epsilonr_annex[sup] = epsr1[] * TensorDiag[1,1,1];
epsilonr_annex[layer_dep] = epsr1[] * TensorDiag[1,1,1];
@@ -242,11 +242,11 @@ Function{
Hr[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
Ht[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
-
+
E1[pmltop] = Ei[]+Er[];
E1[Omega_top] = Ei[]+Er[];
- E1[Omega_bot] = Et[];
- E1[pmlbot] = Et[];
+ E1[Omega_bot] = Et[];
+ E1[pmlbot] = Et[];
E1d[pmltop] = Er[];
E1d[Omega_top] = Er[];
E1d[Omega_bot] = Et[];
@@ -278,7 +278,7 @@ Constraint {
Jacobian {
{ Name JVol ;
- Case {
+ Case {
{ Region All ; Jacobian Vol ; }
}
}
@@ -292,9 +292,9 @@ Jacobian {
Integration {
If (Flag_o2_geom==0)
{ Name Int_1 ;
- Case {
+ Case {
{ Type Gauss ;
- Case {
+ Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 12 ; }
@@ -305,9 +305,9 @@ Integration {
EndIf
If (Flag_o2_geom==1)
{ Name Int_1 ;
- Case {
+ Case {
{ Type Gauss ;
- Case {
+ Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line2 ; NumberOfPoints 4 ; }
{ GeoElement Triangle2 ; NumberOfPoints 12 ; }
@@ -386,7 +386,7 @@ Formulation {
Galerkin {[-Complex[0,gamma[]] /mur[] * Dof{d El} , EZ[]*^{Et} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ Complex[0,gamma[]] /mur[] * (EZ[]*^Dof{Et}) , {d El} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ gamma[]^2 /mur[] * (EZ[]*^Dof{Et}) , EZ[]*^{Et} ] ; In Omega; Integration Int_1; Jacobian JVol; }
-
+
Galerkin {[ -k0^2 * epsilonr[] * Dof{Et} , {Et} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ -k0^2 * epsilonr[] * Dof{El} , {El} ] ; In Omega; Integration Int_1; Jacobian JVol; }
@@ -459,7 +459,7 @@ PostProcessing {
{ Name eff_r~{i} ; Value { Term{Type Global; [
beta_super~{i}[]/(Ae^2*-beta1[]) * ( SquNorm[$int_x_r~{i}]+
SquNorm[$int_y_r~{i}]+
- SquNorm[$int_z_r~{i}] ) ] ; In SurfCutSuper1 ; } } }
+ SquNorm[$int_z_r~{i}] ) ] ; In SurfCutSuper1 ; } } }
// // BUGGY
// { Name eff_t~{i} ; Value{ Term{Type Global; [
// 1/(Ae^2*beta_subs~{i}[]*-beta1[]) * ((beta_subs~{i}[]^2+gamma[]^2 )*SquNorm[$int_z_t~{i}]+
@@ -477,7 +477,7 @@ PostProcessing {
PostOperation {
{ Name postop_energy; NameOfPostProcessing postpro_energy ;
Operation {
- Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
+ Print[ lambda_step, OnPoint{0,0,0}, Format ValueOnly, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
Print [ Etot , OnElementsOf Omega, File "Etot.pos" ];
Print [ Etot , OnPlane { {-d/2,-h_sub,0} {d/2,-h_sub,0} {-d/2,h_layer_dep+dy*N_rods+h_layer_cov+h_sup,0} } {25,50} , File StrCat[myDir,"Etot_grid.pos"]];
For i In {0:2*nb_orders}
diff --git a/DiffractionGratings/grating2D_scalar.pro b/DiffractionGratings/grating2D_scalar.pro
index 44cdaaee668eb49701a56544867c817c2fd76387..8ed72406441df10c9723279941de804d70722bbe 100644
--- a/DiffractionGratings/grating2D_scalar.pro
+++ b/DiffractionGratings/grating2D_scalar.pro
@@ -1,17 +1,17 @@
// Copyright (C) 2020 Guillaume Demésy
//
// This file is part of the model grating2D.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/>.
@@ -122,7 +122,7 @@ Function{
epsr_layer_cov_im[] = epsr_im_dom_5[];
epsr1_re[] = epsr_re_dom_6[];
epsr1_im[] = epsr_im_dom_6[];
-
+
om0 = 2.*Pi*cel/lambda0;
k0 = 2.*Pi/lambda0;
epsr1[] = Complex[epsr1_re[],epsr1_im[]];
@@ -134,7 +134,7 @@ Function{
alpha[] = k0*n1[]*Sin[theta_deg*deg2rad];
beta1[] = -k0*n1[]*Cos[theta_deg*deg2rad];
beta2[] = -Sqrt[k0^2*epsr2[]-alpha[]^2];
-
+
k1[] = Vector[alpha[], beta1[],0];
k1r[]= Vector[alpha[],-beta1[],0];
k2[] = Vector[alpha[], beta2[],0];
@@ -147,7 +147,7 @@ Function{
If (flag_polar==0)
r[] = (CompY[k1[]]-CompY[k2[]])/(CompY[k1[]]+CompY[k2[]]);
t[] = 2.*CompY[k1[]] /(CompY[k1[]]+CompY[k2[]]);
- Pinc[] = 0.5*A^2*Sqrt[ep0*epsr1_re[]/mu0] * Cos[theta_deg*deg2rad];
+ Pinc[] = 0.5*A^2*Sqrt[ep0*epsr1_re[]/mu0] * Cos[theta_deg*deg2rad];
EndIf
BlochX_phase_shift[] = Exp[I[]*alpha[]*d];
@@ -158,11 +158,11 @@ Function{
beta_super~{i}[] = -Sqrt[k0^2*epsr1[]-alpha~{i}[]^2];
beta_subs~{i}[] = -Sqrt[k0^2*epsr2[]-alpha~{i}[]^2];
EndFor
-
+
a_pml = 1.;
b_pml = 1.;
sx = 1.;
- sy[] = Complex[a_pml,b_pml];
+ sy[] = Complex[a_pml,b_pml];
sz = 1.;
PML_Tensor[] = TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz];
@@ -186,8 +186,8 @@ Function{
EndIf
epsilonr[layer_cov] = Complex[epsr_layer_cov_re[],epsr_layer_cov_im[]] * TensorDiag[1,1,1];
epsilonr[pmltop] = epsr1_re[]*PML_Tensor[];
- epsilonr[pmlbot] = epsr2_re[]*PML_Tensor[];
-
+ epsilonr[pmlbot] = epsr2_re[]*PML_Tensor[];
+
epsilonr_annex[sub] = epsr2[] * TensorDiag[1,1,1];
epsilonr_annex[sup] = epsr1[] * TensorDiag[1,1,1];
epsilonr_annex[layer_dep] = epsr1[] * TensorDiag[1,1,1];
@@ -214,7 +214,7 @@ Function{
ui[rods] = A*Exp[I[]*k1[]*XYZ[]];
ui[layer_dep] = A*Exp[I[]*k1[]*XYZ[]];
ui[sub] = 0.;
-
+
ur[pmltop] = 0.;
ur[pmlbot] = 0.;
ur[sup] = r[]*Exp[I[]*k1r[]*XYZ[]];
@@ -232,10 +232,10 @@ Function{
ut[rods] = 0.;
ut[layer_cov] = 0.;
ut[sub] = t[]*Exp[I[]*k2[]*XYZ[]];
-
+
u1[] = ui[]+ur[]+ut[];
u1d[] = ur[]+ut[];
-
+
If (flag_polar==1)
E1i[] = -1/(om0*ep0*epsilonr_annex[]) * k1[] /\ Vector[0,0,ui[]];
E1r[] = -1/(om0*ep0*epsilonr_annex[]) * k1r[] /\ Vector[0,0,ur[]];
@@ -278,7 +278,7 @@ Constraint {
Jacobian {
{ Name JVol ;
- Case {
+ Case {
{ Region All ; Jacobian Vol ; }
}
}
@@ -292,9 +292,9 @@ Jacobian {
Integration {
If (Flag_o2_geom==0)
{ Name Int_1 ;
- Case {
+ Case {
{ Type Gauss ;
- Case {
+ Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 12 ; }
@@ -305,9 +305,9 @@ Integration {
EndIf
If (Flag_o2_geom==1)
{ Name Int_1 ;
- Case {
+ Case {
{ Type Gauss ;
- Case {
+ Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line2 ; NumberOfPoints 4 ; }
{ GeoElement Triangle2 ; NumberOfPoints 12 ; }
@@ -328,7 +328,7 @@ FunctionSpace {
{ Name un4; NameOfCoef un4; Function BF_Node_3E; Support Region[Omega]; Entity EdgesOf[Omega]; }
{ Name un5; NameOfCoef un5; Function BF_Node_3F; Support Region[Omega]; Entity FacetsOf[Omega]; }
EndIf
-
+
}
Constraint {
{ NameOfCoef un; EntityType NodesOf ; NameOfConstraint Bloch; }
@@ -348,7 +348,7 @@ Formulation {
Galerkin { [-xsi[]*Dof{Grad u2d} , {Grad u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
Galerkin { [k0^2*chi[]*Dof{u2d} , {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
Galerkin { [source_xsi_i[] , {Grad u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
- Galerkin { [source_xsi_r[] , {Grad u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source_xsi_r[] , {Grad u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
Galerkin { [source_chi_r[] , {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
Galerkin { [source_chi_i[] , {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
}
@@ -402,7 +402,7 @@ PostProcessing {
{ Name Q_layer_cov ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_cov_im[]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In layer_cov ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name Q_tot ; Value { Integral { [ 0.5 * ep0*om0*Fabs[Im[CompZZ[epsilonr[]]]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompYY[epsilonr[]])+Ey1[]/CompYY[epsilonr[]]*CompYY[epsilonr_annex[]] ] ) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
EndIf
- If (flag_polar==0)
+ If (flag_polar==0)
For i In {0:2*nb_orders}
{ Name eff_r~{i} ; Value { Term{ Type Global; [ SquNorm[$sr~{i}]*beta_super~{i}[]/beta1[] ] ; In SurfCutSuper1 ; } } }
{ Name eff_t~{i} ; Value { Term{ Type Global; [ SquNorm[$st~{i}]*(beta_subs~{i}[]/beta1[])] ; In SurfCutSubs1 ; } } }
@@ -423,8 +423,8 @@ 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" ] ;
+ Print[ SSize, OnPoint{0,0,0}, Format ValueOnly, File > StrCat[myDir, "SSize.txt"]] ;
+ Print[ lambda_step, OnPoint{0,0,0}, Format ValueOnly, 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)]];
Print[ s_t~{i}[SurfCutSubs1] , OnGlobal, StoreInVariable $st~{i}, Format Table , File > StrCat[myDir, Sprintf("temp_s_t_%g.txt", i-nb_orders)]];
@@ -432,15 +432,15 @@ PostOperation {
For i In {0:2*nb_orders}
Print[ eff_r~{i}[SurfCutSuper1], OnRegion SurfCutSuper1,Store (4*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_r_%g.txt", i-nb_orders)]];
Print[ eff_t~{i}[SurfCutSubs1] , OnRegion SurfCutSubs1 ,Store (6*nb_orders+1+i), Format FrequencyTable, File > StrCat[myDir, Sprintf("efficiency_t_%g.txt", i-nb_orders)]];
- Print[ order_r_angle~{i} , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_r_angle_%g.txt", i-nb_orders)]];
- Print[ order_t_angle~{i} , OnPoint{0,0,0}, Format Table , File > StrCat[myDir, Sprintf("order_t_angle_%g.txt", i-nb_orders)]];
+ Print[ order_r_angle~{i} , OnPoint{0,0,0}, Format ValueOnly , File > StrCat[myDir, Sprintf("order_r_angle_%g.txt", i-nb_orders)]];
+ Print[ order_t_angle~{i} , OnPoint{0,0,0}, Format ValueOnly , File > StrCat[myDir, Sprintf("order_t_angle_%g.txt", i-nb_orders)]];
EndFor
Print[ eff_r~{nb_orders}[SurfCutSuper1], OnRegion SurfCutSuper1, Format Table, SendToServer "GetDP/R0", File StrCat[myDir, "temp_R0.txt"]];
Print[ eff_t~{nb_orders}[SurfCutSubs1] , OnRegion SurfCutSubs1 , Format Table, SendToServer "GetDP/T0", File StrCat[myDir, "temp_T0.txt"]];
Print[ Q_tot[Plot_domain] , OnGlobal , Format FrequencyTable ,SendToServer "GetDP/total absorption", File > StrCat[myDir, "absorption-Q_tot.txt"]];
For i In {0:N_rods-1:1}
- Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];
- EndFor
+ Print[ Q_rod~{i}[rod~{i}] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, Sprintf("absorption-Q_rod_%g.txt", i+1) ]];
+ EndFor
Print[ Q_tot[Plot_domain] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_tot.txt"]];
Print[ Q_subs[sub] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_subs.txt"]];
Print[ Q_rod_out[rod_out] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_rod_out.txt"]];
diff --git a/DiffractionGratings/grating3D.pro b/DiffractionGratings/grating3D.pro
index 4bd19e1d18adac76588566a938a58f621a496a42..2247a03a4763eafcc0cb2ba621f36489a588b7f4 100644
--- a/DiffractionGratings/grating3D.pro
+++ b/DiffractionGratings/grating3D.pro
@@ -17,9 +17,9 @@
Include "grating3D_data.geo";
Include "grating3D_materials.pro";
-
+
myDir = "res3D/";
-
+
Group {
// SubDomains
PMLbot = Region[1];
@@ -31,7 +31,7 @@
L_1 = Region[7];
PMLtop = Region[8];
Scat = Region[9];
-
+
// Boundaries
SurfBlochXm = Region[101];
SurfBlochXp = Region[102];
@@ -39,41 +39,41 @@
SurfBlochYp = Region[202];
SurfIntTop = Region[301];
SurfIntBot = Region[302];
-
-
+
+
SurfDirichlet = Region[{401,402}];
SurfBloch = Region[{SurfBlochXm,SurfBlochXp,SurfBlochYm,SurfBlochYp}];
-
+
If (FlagLinkFacets==1)
SurfExcludeFacets = Region[{}];
Else
SurfExcludeFacets = Region[{SurfBloch}];
EndIf
-
+
// 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}];
-
+
// get normal comp of E field on integ surfaces
Gama = Region[{SurfIntBot,SurfIntTop}];
Tr = ElementsOf[Omega_plot, ConnectedTo Gama];
}
-
-
-
+
+
+
Function{
I[] = Complex[0,1];
zhat[] = Vector[0,0,1];
-
+
ispecular = Nmax;
jspecular = Nmax;
-
+
small_delta = 0.0*nm;
mu0 = 4*Pi*1.e2*nm;
ep0 = 8.854187817e-3*nm;
@@ -82,7 +82,7 @@
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)
@@ -107,7 +107,7 @@
EndFor
EndIf
EndFor
-
+
If (flag_mat_scat==0)
epsr[Scat] = Complex[eps_re_Scat , eps_im_Scat] * TensorDiag[1,1,1];
Else
@@ -115,13 +115,13 @@
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];
mur[SurfIntTop] = TensorDiag[1,1,1];
-
+
////// PMLS
a_pml = 1.;
b_pml = 1.;
@@ -130,14 +130,14 @@
n2[] = Sqrt[epsr2[]];
k1norm[] = k0*n1[];
k2norm[] = k0*n2[];
-
+
Zmax = PML_top_hh;
Zmin = hh_L_6;
Damp_pml_top[] = 1/(Zmax + PML_top - Fabs[Z[]]) - 1/(PML_top);
Damp_pml_bot[] = 1/(Zmin + PML_top - Fabs[Z[]]) - 1/(PML_bot);
Sigma_top[] = 0.5*(Damp_pml_top[] + Fabs[Damp_pml_top[]]);
Sigma_bot[] = 0.5*(Damp_pml_bot[] + Fabs[Damp_pml_bot[]]);
-
+
If (PML_TYPE==0)
sz[] = Complex[a_pml,b_pml];
Else
@@ -146,26 +146,26 @@
EndIf
sx = 1.;
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] = Re[epsr1[]]*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] = Re[epsr2[]]*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[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];
epsr_annex[L_1] = epsr[];
epsr_annex[L_6] = epsr[];
-
+
//// Reference Field solution of annex problem (simple diopter)
k1x[] = -k0*n1[]*Sin[theta0]*Cos[phi0];
k1y[] = -k0*n1[]*Sin[theta0]*Sin[phi0];
@@ -176,23 +176,23 @@
k1[] = Vector[k1x[],k1y[], k1z[]];
k2[] = Vector[k2x[],k2y[], k2z[]];
k1r[] = Vector[k1x[],k1y[],-k1z[]];
-
+
rs[] = (k1z[]-k2z[])/(k1z[]+k2z[]);
ts[] = 2.*k1z[] /(k1z[]+k2z[]);
rp[] = (k1z[]*epsr2[]-k2z[]*epsr1[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
tp[] = (2.*k1z[]*epsr2[])/(k1z[]*epsr2[]+k2z[]*epsr1[]);
-
+
spol[] = Vector[Sin[phi0],-Cos[phi0],0];
ppol_r[] = (k1r[]/Norm[k1r[]]) /\ spol[];
ppol_t[] = (k2[] /Norm[k2[]] ) /\ spol[];
-
+
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[]];
@@ -202,14 +202,14 @@
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[SurfIntTop] = Ei[]+Er[];
E1[Omega_super] = Ei[]+Er[];
E1[Omega_subs] = Et[];
@@ -218,22 +218,22 @@
E1d[Omega_super] = Er[];
E1d[Omega_subs] = Et[];
E1d[SurfIntBot] = Et[];
-
+
H1[Omega_super] = Hi[]+Hr[];
H1[Omega_subs] = Ht[];
H1d[Omega_super] = Hr[];
H1d[Omega_subs] = Ht[];
-
+
source_vol_scat[] = (om0/cel)^2*(epsr[]-epsr_annex[])*E1[];
source_surf_tot[] = -2*I[]*om0*mu0*Vector[0,0,1] /\ Hi[];
// Bloch phase shifts
skx1[] = k1x[];
// sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
sky1[] = -k0*n1[]*Sin[theta0]*Sin[phi0+xsi];
-
+
dephX[] = Exp[I[]*skx1[]*period_x];
dephY[] = Exp[I[]*sky1[]*period_y];
-
+
// Fourier coefficients variables
Nb_ordre = 2*Nmax+1;
For i In {0:Nb_ordre-1}
@@ -249,9 +249,9 @@
gammat~{i}~{j}[] = Sqrt[k0^2*epsr2[] - alpha~{i}~{j}[]^2 - beta~{i}~{j}[]^2];
EndFor
EndFor
-
+
}
-
+
Constraint {
{ Name Dirichlet; Type Assign;
Case {
@@ -271,13 +271,13 @@
}
}
}
-
+
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 {
@@ -294,7 +294,7 @@
}
}
}
-
+
FunctionSpace {
{ Name Hcurl; Type Form1;
BasisFunction {
@@ -335,7 +335,7 @@
}
}
}
-
+
Formulation {
{ Name helmholtz_vector; Type FemEquation;
Quantity {
@@ -355,7 +355,7 @@
}
}
}
-
+
Resolution {
{ Name helmholtz_vector;
System {
@@ -368,7 +368,7 @@
}
}
}
-
+
PostProcessing {
{ Name postpro_helmholtz_vector; NameOfFormulation helmholtz_vector; NameOfSystem M;
Quantity {
@@ -380,7 +380,7 @@
{ Name Poy_inc; Value { Local { [ 0.5*Re[Cross[ Ei[] , Conj[Hi[]]]] ]; In Omega; Jacobian JVol; } } }
{ Name E1 ; Value { Local { [ E1[] ]; In Omega; Jacobian JVol; } } }
{ Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
-
+
If (FLAG_TOTAL==0)
{ 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; } } }
@@ -400,7 +400,7 @@
EndFor
{ Name Abs_scat ; Value { Integral { [ ep0*om0 * 0.5*Im[CompXX[epsr[]]]*(SquNorm[{u}]) / (Pinc*period_x*dyc) ] ; In Scat ; Integration I1 ; Jacobian JVol ; } } }
EndIf
-
+
For i In {0:Nb_ordre-1}
For j In {0:Nb_ordre-1}
If (FLAG_TOTAL==0)
@@ -442,26 +442,26 @@
// Is it better to compute it at the surface on which the scatterer is relying? (so that if there is no scatterer,
// it just corresponds to the usual definition of rs/rp for a simple diopter). Maybe... uncomment phasor if necessary.
// For the Mmatrix, we do not care about the phase.
- { Name er_specular ; Value { Term{ Type Global; [
- // Exp[I[]*k1z[]*(thick_L_1+thick_L_2+thick_L_3)] *
+ { Name er_specular ; Value { Term{ Type Global; [
+ // Exp[I[]*k1z[]*(thick_L_1+thick_L_2+thick_L_3)] *
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; [
+ { 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_r[] * $er_specular] ; In SurfIntTop ; } } }
{ Name rs ; Value { Term{ Type Global; [ spol[] * $er_specular] ; In SurfIntTop ; } } }
{ Name tp ; Value { Term{ Type Global; [ ppol_t[] * $et_specular] ; In SurfIntBot ; } } }
{ Name ts ; Value { Term{ Type Global; [ spol[] * $et_specular] ; In SurfIntBot ; } } }
-
+
}
}
}
-
+
PostOperation {
{ Name postop_helmholtz_vector; NameOfPostProcessing postpro_helmholtz_vector ;
Operation {
@@ -470,13 +470,13 @@
// Print [ uz , OnElementsOf SurfIntBot, File StrCat[myDir,"uz_ZM.pos"]];
// Print [ CompZu , OnPlane { {-period_x/2,-period_y/2,hh_L_6+0.5*nm} { period_x/2,-period_y/2,hh_L_6+0.5*nm} {-period_x/2, period_y/2,hh_L_6+0.5*nm} } {npts_interpX,npts_interpY} , File StrCat[myDir,"u_cut_ZM.pos"]];
// Print [ CompZu , OnPlane { {-period_x/2,-period_y/2,hh_L_1+thick_L_1-0.5*nm} { period_x/2,-period_y/2,hh_L_1+thick_L_1-0.5*nm} {-period_x/2, period_y/2,hh_L_1+thick_L_1-0.5*nm} } {npts_interpX,npts_interpY} , File StrCat[myDir,"u_cut_ZP.pos"]];
- // // Debug : print opto-geometric parameters
+ // // Debug : print opto-geometric parameters
// Print [ epsr_xx , OnElementsOf Omega, File StrCat[myDir,"epsr_xx.pos"]];
// // Debug : print raw u and Etot
// Print [ u , OnElementsOf Omega, File StrCat[myDir,"Edif.pos"]];
// Print [ Etot , OnElementsOf Omega, File StrCat[myDir,"Etot.pos"]];
-
-
+
+
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"];
@@ -516,7 +516,7 @@
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
-
+
Print [ Poy_tot , OnPlane { {0.5*(-period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
{0.5*( period_x-dys), -dyc/2,(hh_L_6+hh_L_5)/2}
{0.5*(-period_x+dys), dyc/2,(hh_L_6+hh_L_5)/2} }
@@ -529,23 +529,23 @@
{0.5*( period_x-dys), -dyc/2, hh_L_1+thick_L_1/2}
{0.5*(-period_x+dys), dyc/2, hh_L_1+thick_L_1/2} }
{npts_checkpoyX-1,npts_checkpoyY-1} , File StrCat[myDir,"Poy_inc_gd.pos"], Format Table];
-
+
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_z_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_z_t~{i}~{j}, Format Table];
+ Print[ int_z_t~{i}~{j}[SurfIntBot], OnGlobal, StoreInVariable $int_z_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];
Print[ int_z_r~{i}~{j}[SurfIntTop], OnGlobal, StoreInVariable $int_z_r~{i}~{j}, Format Table];
EndFor
EndFor
-
+
For i In {0:Nb_ordre-1}
For j In {0:Nb_ordre-1}
Print[ eff_t1~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir, "eff_t1.txt"], Format Table ];
@@ -565,14 +565,13 @@
Print[ numbering_ij~{i}~{j}[SurfIntBot], OnRegion SurfIntBot, File > StrCat[myDir,"numbering_ij.txt"], Format Table ];
EndFor
EndFor
- Print[ lambda_step, OnPoint{0,0,0}, Format Table, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
+ Print[ lambda_step, OnPoint{0,0,0}, Format ValueOnly, File > StrCat[myDir, "temp_lambda_step.txt"], SendToServer "GetDP/Lambda_step" ] ;
}
}
}
-
+
DefineConstant[
R_ = {"helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
C_ = {"-solve -pos -petsc_prealloc 500 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -ksp_error_if_not_converged", Name "GetDP/9ComputeCommand", Visible 1},
P_ = {"postop_helmholtz_vector", Name "GetDP/2PostOperationChoices", Visible 1}
];
-
\ No newline at end of file
diff --git a/Magnetometer/magnetometer.pro b/Magnetometer/magnetometer.pro
index 12b7cf5072177f6004388b76d7a8d420e72e8ff8..8f35090e6ea4563f8f81fa3c8d9d924c91ad324d 100644
--- a/Magnetometer/magnetometer.pro
+++ b/Magnetometer/magnetometer.pro
@@ -231,7 +231,7 @@ PostOperation {
Print[ um, OnPoint {l/2,a/2,b/2}, Format ValueOnly, File "res/um_middle.txt",
SendToServer "Output/Middle diplacement [m]", Color "LightYellow" ];
Echo[ Str["l=PostProcessing.NbViews-1; View[l].VectorType=5; ",
- "View[l].DisplacementFactor = 1e10;"],
+ "View[l].DisplacementFactor = 1e07;"],
File "res/tmp.geo", LastTimeStepOnly] ;
EndIf
}
diff --git a/Magnetostriction/choke.pro b/Magnetostriction/choke.pro
index c7009ca18311c6d3f07ca7195a6af06a992ff349..184e98b380db8a42c92a70ad9ec20b44a44c5483 100644
--- a/Magnetostriction/choke.pro
+++ b/Magnetostriction/choke.pro
@@ -184,15 +184,15 @@ PostOperation {
// Print [ eps_N, OnElementsOf Domain_Disp, File StrCat [Res_dir,"eps_N.pos" ]] ;
// Print [ Fmaxwell, OnElementsOf Domain_Disp, File StrCat [Res_dir,"Fmaxwell.pos" ]] ;
Print[ bn, OnElementsOf Domain_Mag, File StrCat [Res_dir,"bn.pos" ]] ;
- Print [ u_N, OnPoint {x_measurement,y_measurement-1e-3,0},
+ Print [ u_N, OnPoint {x_measurement,y_measurement-1e-3,0}, Format ValueOnly,
LastTimeStepOnly,
File StrCat [Res_dir,"u_measurement_N.dat"],
SendToServer "Output/displacement |u| [m]", Color "Ivory"];
- Print [ u_x, OnPoint {x_measurement,y_measurement-1e-3,0},
+ Print [ u_x, OnPoint {x_measurement,y_measurement-1e-3,0}, Format ValueOnly,
LastTimeStepOnly,
File StrCat [Res_dir,"u_measurement_x.dat" ],
SendToServer "Output/Displacement X |u_x| [m]", Color "Ivory"];
- Print [ u_y, OnPoint {x_measurement,y_measurement-1e-3,0},
+ Print [ u_y, OnPoint {x_measurement,y_measurement-1e-3,0}, Format ValueOnly,
LastTimeStepOnly,
File StrCat [Res_dir,"u_measurement_y.dat" ],
SendToServer "Output/displacement Y |u_y| [m]", Color "Ivory"];
diff --git a/Shielding/formulations.pro b/Shielding/formulations.pro
index 8f59580d2bd8ef6fbf8fb18959361f3ca4f72b59..cdeb6c5363ea4f55d77f647956ad31e8ef0f9363 100644
--- a/Shielding/formulations.pro
+++ b/Shielding/formulations.pro
@@ -168,7 +168,7 @@ PostOperation {
}
{ Name Get_ShieldingEffectiveness ; NameOfPostProcessing Microwave_e ;
Operation {
- Print[ SE, OnPoint {0,0,0}, Format Table, File StrCat[myDir,"temp",ExtGnuplot],
+ Print[ SE, OnPoint {0,0,0}, Format ValueOnly, File StrCat[myDir,"temp",ExtGnuplot],
SendToServer StrCat(po,"0Shielding effectiveness [dB]")];
}
}