Select Git revision
grating2D_scalar.pro
-
Guillaume Demesy authoredGuillaume Demesy authored
grating2D_scalar.pro 25.24 KiB
// 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/>.
myDir = "res2D/";
DefineConstant[
lambda0 = {lambda_min , Min lambda_min, Max lambda_max, Step (lambda_max-lambda_min)/(nb_lambdas-1), Name StrCat[pp2, "1wavelength [nm]"] , Loop 1, Highlight Str[colorpp2],Graph "200000200020", ServerAction "Reset GetDP/T0, GetDP/R0, GetDP/Lambda_step, GetDP/total absorption"}
];
lambda0 = lambda0 * nm;
lambda_min = lambda_min * nm;
lambda_max = lambda_max * nm;
A = 1.;
nb_orders = 2;
Group {
// Boundaries
SurfBlochLeft = Region[SURF_BLOCH_X_LEFT];
SurfBlochRight = Region[SURF_BLOCH_X_RIGHT];
SurfCutSubs1 = Region[SURF_INTEG_SUB1];
SurfCutSuper1 = Region[SURF_INTEG_SUP1];
SurfCutSubs2 = Region[SURF_INTEG_SUB2];
SurfCutSuper2 = Region[SURF_INTEG_SUP2];
// Domains
pmlbot = Region[PMLBOT];
sub = Region[{SUB,SURF_INTEG_SUB1,SURF_INTEG_SUB2}];
layer_dep = Region[LAYERDEP];
rod_out = Region[RODOUT];
For i In {0:N_rods-1:1}
rod~{i} = Region[{ROD~{i}}];
rods += Region[{ROD~{i}}];
EndFor
layer_cov = Region[LAYERCOV];
sup = Region[{SUP,SURF_INTEG_SUP1,SURF_INTEG_SUP2}];
pmltop = Region[PMLTOP];
Omega_source = Region[{layer_dep,rod_out,rods,layer_cov}];
Omega_nosource = Region[{pmltop,pmlbot,sup,sub}];
Omega = Region[{Omega_source,Omega_nosource}];
Omega_top = Region[{layer_dep,rod_out,rods,layer_cov,sup}];
Omega_bot = Region[{sub}];
Omega_pml = Region[{pmltop,pmlbot}];
Plot_domain = Region[{Omega,-pmltop,-pmlbot}];
Plot_bnd = Region[SURF_PLOT];
Print_point = Region[PRINT_POINT];
}
Function{
I[] = Complex[0.0,1.0];
bndCol[Plot_bnd] = Complex[1,1];
epsr_re_interp_mat_1[] = InterpolationLinear[$1]{ListAlt[lambdamat_1 ,epsr_mat_re_1 ]};
epsr_im_interp_mat_1[] = InterpolationLinear[$1]{ListAlt[lambdamat_1 ,epsr_mat_im_1 ]};
epsr_re_interp_mat_2[] = InterpolationLinear[$1]{ListAlt[lambdamat_2 ,epsr_mat_re_2 ]};
epsr_im_interp_mat_2[] = InterpolationLinear[$1]{ListAlt[lambdamat_2 ,epsr_mat_im_2 ]};
epsr_re_interp_mat_3[] = InterpolationLinear[$1]{ListAlt[lambdamat_3 ,epsr_mat_re_3 ]};
epsr_im_interp_mat_3[] = InterpolationLinear[$1]{ListAlt[lambdamat_3 ,epsr_mat_im_3 ]};
epsr_re_interp_mat_4[] = InterpolationLinear[$1]{ListAlt[lambdamat_4 ,epsr_mat_re_4 ]};
epsr_im_interp_mat_4[] = InterpolationLinear[$1]{ListAlt[lambdamat_4 ,epsr_mat_im_4 ]};
epsr_re_interp_mat_5[] = InterpolationLinear[$1]{ListAlt[lambdamat_5 ,epsr_mat_re_5 ]};
epsr_im_interp_mat_5[] = InterpolationLinear[$1]{ListAlt[lambdamat_5 ,epsr_mat_im_5 ]};
epsr_re_interp_mat_6[] = InterpolationLinear[$1]{ListAlt[lambdamat_6 ,epsr_mat_re_6 ]};
epsr_im_interp_mat_6[] = InterpolationLinear[$1]{ListAlt[lambdamat_6 ,epsr_mat_im_6 ]};
epsr_re_interp_mat_7[] = InterpolationLinear[$1]{ListAlt[lambdamat_7 ,epsr_mat_re_7 ]};
epsr_im_interp_mat_7[] = InterpolationLinear[$1]{ListAlt[lambdamat_7 ,epsr_mat_im_7 ]};
epsr_re_interp_mat_8[] = InterpolationLinear[$1]{ListAlt[lambdamat_8 ,epsr_mat_re_8 ]};
epsr_im_interp_mat_8[] = InterpolationLinear[$1]{ListAlt[lambdamat_8 ,epsr_mat_im_8 ]};
epsr_re_interp_mat_9[] = InterpolationLinear[$1]{ListAlt[lambdamat_9 ,epsr_mat_re_9 ]};
epsr_im_interp_mat_9[] = InterpolationLinear[$1]{ListAlt[lambdamat_9 ,epsr_mat_im_9 ]};
epsr_re_interp_mat_10[] = InterpolationLinear[$1]{ListAlt[lambdamat_10,epsr_mat_re_10]};
epsr_im_interp_mat_10[] = InterpolationLinear[$1]{ListAlt[lambdamat_10,epsr_mat_im_10]};
epsr_re_interp_mat_11[] = InterpolationLinear[$1]{ListAlt[lambdamat_11,epsr_mat_re_11]};
epsr_im_interp_mat_11[] = InterpolationLinear[$1]{ListAlt[lambdamat_11,epsr_mat_im_11]};
epsr_re_interp_mat_12[] = InterpolationLinear[$1]{ListAlt[lambdamat_12,epsr_mat_re_12]};
epsr_im_interp_mat_12[] = InterpolationLinear[$1]{ListAlt[lambdamat_12,epsr_mat_im_12]};
epsr_re_interp_mat_13[] = InterpolationLinear[$1]{ListAlt[lambdamat_13,epsr_mat_re_13]};
epsr_im_interp_mat_13[] = InterpolationLinear[$1]{ListAlt[lambdamat_13,epsr_mat_im_13]};
epsr_re_interp_mat_14[] = InterpolationLinear[$1]{ListAlt[lambdamat_14,epsr_mat_re_14]};
epsr_im_interp_mat_14[] = InterpolationLinear[$1]{ListAlt[lambdamat_14,epsr_mat_im_14]};
epsr_re_interp_mat_15[] = InterpolationLinear[$1]{ListAlt[lambdamat_15,epsr_mat_re_15]};
epsr_im_interp_mat_15[] = InterpolationLinear[$1]{ListAlt[lambdamat_15,epsr_mat_im_15]};
epsr_re_interp_mat_16[] = InterpolationLinear[$1]{ListAlt[lambdamat_16,epsr_mat_re_16]};
epsr_im_interp_mat_16[] = InterpolationLinear[$1]{ListAlt[lambdamat_16,epsr_mat_im_16]};
For i In {1:5}
For j In {1:nb_available_materials}
If(flag_mat~{i}==j-1)
epsr_re_dom~{i}[] = epsr_re_interp_mat~{j}[lambda0/nm*1e-9];
epsr_im_dom~{i}[] = epsr_im_interp_mat~{j}[lambda0/nm*1e-9];
EndIf
EndFor
EndFor
For k In {0:nb_available_lossless_materials-1}
If(flag_mat_6==lossless_material_list(k)-1)
epsr_re_dom_6[] = epsr_re_interp_mat~{lossless_material_list(k)}[lambda0/nm*1e-9];
epsr_im_dom_6[] = epsr_im_interp_mat~{lossless_material_list(k)}[lambda0/nm*1e-9];
EndIf
EndFor
epsr2_re[] = epsr_re_dom_1[];
epsr2_im[] = epsr_im_dom_1[];
epsr_layer_dep_re[] = epsr_re_dom_2[];
epsr_layer_dep_im[] = epsr_im_dom_2[];
epsr_rod_out_re[] = epsr_re_dom_3[];
epsr_rod_out_im[] = epsr_im_dom_3[];
epsr_rods_re[] = epsr_re_dom_4[];
epsr_rods_im[] = epsr_im_dom_4[];
epsr_layer_cov_re[] = epsr_re_dom_5[];
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[]];
epsr2[] = Complex[epsr2_re[],epsr2_im[]];
n1[] = Sqrt[epsr1[]];
n2[] = Sqrt[epsr2[]];
k1norm[]= k0*n1[];
k2norm[]= k0*n2[];
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];
If (flag_polar==1)
r[] = (CompY[k1[]]*epsr2[]-CompY[k2[]]*epsr1[])/(CompY[k1[]]*epsr2[]+CompY[k2[]]*epsr1[]);
t[] = (2.*CompY[k1[]]*epsr2[])/(CompY[k1[]]*epsr2[]+CompY[k2[]]*epsr1[]);
Pinc[] = 0.5*A^2*Sqrt[mu0/(ep0*epsr1_re[])] * Cos[theta_deg*deg2rad];
EndIf
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];
EndIf
BlochX_phase_shift[] = Exp[I[]*alpha[]*d];
For i In {0:2*nb_orders}
alpha~{i}[] = CompX[k1[]] + 2.*Pi/d*(i-nb_orders);
expmialpha~{i}[] = Exp[-I[]*alpha~{i}[]*X[]];
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];
sz = 1.;
PML_Tensor[] = TensorDiag[sz*sy[]/sx,sx*sz/sy[],sx*sy[]/sz];
epsilonr[sup] = epsr1[] * TensorDiag[1,1,1];
epsilonr[sub] = epsr2[] * TensorDiag[1,1,1];
epsilonr[layer_dep] = Complex[epsr_layer_dep_re[],epsr_layer_dep_im[]] * TensorDiag[1,1,1];
epsilonr[rod_out] = Complex[epsr_rod_out_re[],epsr_rod_out_im[] ] * TensorDiag[1,1,1];
If (Flag_aniso==0)
epsilonr[rods] = Complex[epsr_rods_re[],epsr_rods_im[]] * TensorDiag[1,1,1];
Else
epsilonr[rods] = Tensor[Complex[epsr_custom_anisoXX_re, epsr_custom_anisoXX_im],
Complex[epsr_custom_anisoXY_re, epsr_custom_anisoXY_im],
0,
Complex[epsr_custom_anisoXY_re,-epsr_custom_anisoXX_im],
Complex[epsr_custom_anisoYY_re,epsr_custom_anisoYY_im],
0,
0,
0,
Complex[epsr_custom_anisoZZ_re,epsr_custom_anisoZZ_im]
];
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_annex[sub] = epsr2[] * TensorDiag[1,1,1];
epsilonr_annex[sup] = epsr1[] * TensorDiag[1,1,1];
epsilonr_annex[layer_dep] = epsr1[] * TensorDiag[1,1,1];
epsilonr_annex[rod_out] = epsr1[] * TensorDiag[1,1,1];
epsilonr_annex[rods] = epsr1[] * TensorDiag[1,1,1];
epsilonr_annex[layer_cov] = epsr1[] * TensorDiag[1,1,1];
epsilonr_annex[pmltop] = epsr1_re[]*PML_Tensor[];
epsilonr_annex[pmlbot] = epsr2_re[]*PML_Tensor[];
mur[pmltop] = PML_Tensor[];
mur[pmlbot] = PML_Tensor[];
mur[sub] = TensorDiag[1,1,1];
mur[sup] = TensorDiag[1,1,1];
mur[layer_dep] = TensorDiag[1,1,1];
mur[rods] = TensorDiag[1,1,1];
mur[rod_out] = TensorDiag[1,1,1];
mur[layer_cov] = TensorDiag[1,1,1];
ui[pmltop] = 0.;
ui[pmlbot] = 0.;
ui[sup] = A*Exp[I[]*k1[]*XYZ[]];
ui[layer_cov] = A*Exp[I[]*k1[]*XYZ[]];
ui[rod_out] = A*Exp[I[]*k1[]*XYZ[]];
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[]];
ur[layer_dep] = r[]*Exp[I[]*k1r[]*XYZ[]];
ur[rod_out] = r[]*Exp[I[]*k1r[]*XYZ[]];
ur[rods] = r[]*Exp[I[]*k1r[]*XYZ[]];
ur[layer_cov] = r[]*Exp[I[]*k1r[]*XYZ[]];
ur[sub] = 0;
ut[pmltop] = 0.;
ut[pmlbot] = 0.;
ut[sup] = 0.;
ut[layer_dep] = 0.;
ut[rod_out] = 0.;
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[]];
E1t[] = -1/(om0*ep0*epsilonr_annex[]) * k2[] /\ Vector[0,0,ut[]];
E1[] = E1i[]+E1r[]+E1t[];
Ex1[] = CompX[E1[]];
Ey1[] = CompY[E1[]];
detepst[] = CompXX[epsilonr[]]*CompYY[epsilonr[]]-CompXY[epsilonr[]]*Conj[CompYX[epsilonr[]]];
detepst_annex[] = CompXX[epsilonr_annex[]]*CompYY[epsilonr_annex[]]-CompXY[epsilonr_annex[]]*Conj[CompYX[epsilonr_annex[]]];
xsi[] = Transpose[epsilonr[]]/detepst[];
xsi_annex[] = Transpose[epsilonr_annex[]]/detepst_annex[];
chi[] = CompZZ[mur[]];
source_xsi_r[] = (xsi[]-xsi_annex[]) * -k1[] * I[] * ui[];
source_xsi_i[] = (xsi[]-xsi_annex[]) * -k1r[] * I[] * ur[];
source_chi_r[] = 0;
source_chi_i[] = 0;
Else
detmut[] = CompXX[mur[]]*CompYY[mur[]]-CompXY[mur[]]*Conj[CompYX[mur[]]];
xsi[] = Transpose[mur[]]/detmut[];
chi[] = CompZZ[epsilonr[]];
chi_annex[] = CompZZ[epsilonr_annex[]];
source_xsi_r[] = Vector[0.,0.,0.];
source_xsi_i[] = Vector[0.,0.,0.];
source_chi_r[] = k0^2 * (chi[]-chi_annex[]) * ur[];
source_chi_i[] = k0^2 * (chi[]-chi_annex[]) * ui[];
EndIf
}
Constraint {
{Name Bloch;
Case {
{ Region SurfBlochRight; Type LinkCplx ; RegionRef SurfBlochLeft;
Coefficient BlochX_phase_shift[]; Function Vector[$X-d,$Y,$Z] ;
}
}
}
}
Jacobian {
{ Name JVol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name JSur ;
Case {
{ Region All ; Jacobian Sur ; }
}
}
}
Integration {
If (Flag_o2_geom==0)
{ Name Int_1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 12 ; }
}
}
}
}
EndIf
If (Flag_o2_geom==1)
{ Name Int_1 ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line2 ; NumberOfPoints 4 ; }
{ GeoElement Triangle2 ; NumberOfPoints 12 ; }
}
}
}
}
EndIf
}
FunctionSpace {
{ Name Hgrad; Type Form0;
BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Node; Support Region[Omega]; Entity NodesOf[Omega]; }
{ Name sn2; NameOfCoef un2; Function BF_Node_2E; Support Region[Omega]; Entity EdgesOf[Omega]; }
If (Flag_o2_interp==1)
{ Name un3; NameOfCoef un3; Function BF_Node_2F; Support Region[Omega]; Entity FacetsOf[Omega]; }
{ 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; }
{ NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Bloch; }
If (Flag_o2_interp==1)
{ NameOfCoef un3; EntityType EdgesOf ; NameOfConstraint Bloch; }
EndIf
}
}
}
Formulation {
{Name helmoltz_scalar; Type FemEquation;
Quantity {
{ Name u2d; Type Local; NameOfSpace Hgrad;}}
Equation {
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_chi_r[] , {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
Galerkin { [source_chi_i[] , {u2d}]; In Omega; Jacobian JVol; Integration Int_1; }
}
}
}
Resolution {
{ Name helmoltz_scalar;
System {
{ Name S; NameOfFormulation helmoltz_scalar; Type ComplexValue; }
}
Operation {
CreateDir[Str[myDir]];
Printf["%f",lambda0/nm];
Generate[S] ;
Solve[S] ;
}
}
}
PostProcessing {
{ Name postpro_energy; NameOfFormulation helmoltz_scalar;
Quantity {
{ Name u ; Value { Local { [ {u2d} ]; In Omega; Jacobian JVol; } } }
{ Name u_diff ; Value { Local { [ {u2d}+u1d[] ]; In Omega; Jacobian JVol; } } }
{ Name u_tot ; Value { Local { [ {u2d}+u1[] ]; In Omega; Jacobian JVol; } } }
{ Name u1 ; Value { Local { [ u1[] ]; In Omega; Jacobian JVol; } } }
{ Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega; Jacobian JVol; } } }
{ Name boundary ; Value { Local { [ bndCol[] ]; In Plot_bnd ; Jacobian JVol ; } } }
For i In {-nb_plot_periods:nb_plot_periods}
{ Name u_tot~{i} ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*i*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
EndFor
For i In {0:2*nb_orders}
{ Name s_r~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSuper1 ; Jacobian JSur ; Integration Int_1 ; } } }
{ Name s_t~{i} ; Value { Integral{ [ expmialpha~{i}[] * ({u2d}+u1d[])/d ] ; In SurfCutSubs1 ; Jacobian JSur ; Integration Int_1 ; } } }
{ Name order_t_angle~{i} ; Value { Local{ [-Atan2[Re[alpha~{i}[]],Re[beta_subs~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
{ Name order_r_angle~{i} ; Value { Local{ [ Atan2[Re[alpha~{i}[]],Re[beta_super~{i}[]]]/deg2rad ] ; In Omega; Jacobian JVol; } } }
EndFor
If (flag_polar==1)
For i In {0:2*nb_orders}
{ Name eff_r~{i} ; Value { Term{ Type Global; [ SquNorm[#i]*beta_super~{i}[]/beta1[] ] ; In SurfCutSuper1 ; } } }
{ Name eff_t~{i} ; Value { Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(beta_subs~{i}[]/beta1[])*(epsr1[]/epsr2[]) ] ; In SurfCutSubs1 ; } } }
EndFor
For i In {0:N_rods-1:1}
{ Name Q_rod~{i} ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_rods_im[]] * ( SquNorm[CompY[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ex1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] + SquNorm[-CompX[{Grad u2d}]*I[]/(om0*ep0*CompXX[epsilonr[]])+Ey1[]/CompXX[epsilonr[]]*CompXX[epsilonr_annex[]]] ) / (Pinc[]*d) ] ; In rod~{i} ; Integration Int_1 ; Jacobian JVol ; } } }
EndFor
{ Name Q_subs ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr2_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 sub ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name Q_rod_out ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_rod_out_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 rod_out ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name Q_layer_dep ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_dep_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_dep ; Integration Int_1 ; Jacobian JVol ; } } }
{ 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)
For i In {0:2*nb_orders}
{ Name eff_r~{i} ; Value { Term{ Type Global; [ SquNorm[#i]*beta_super~{i}[]/beta1[] ] ; In SurfCutSuper1 ; } } }
{ Name eff_t~{i} ; Value { Term{ Type Global; [ SquNorm[#(2*nb_orders+1+i)]*(beta_subs~{i}[]/beta1[])] ; In SurfCutSubs1 ; } } }
EndFor
For i In {0:N_rods-1:1}
{ Name Q_rod~{i} ; Value { Integral { [0.5 * ep0*om0*Fabs[epsr_rods_im[]] * (SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod~{i} ; Integration Int_1 ; Jacobian JVol ; } } }
EndFor
{ Name Q_subs ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr2_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In sub ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name Q_rod_out ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_rod_out_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In rod_out ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name Q_layer_dep ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_dep_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_dep ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name Q_layer_cov ; Value { Integral { [ 0.5 * ep0*om0*Fabs[epsr_layer_cov_im[]] *(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In layer_cov ; Integration Int_1 ; Jacobian JVol ; } } }
{ Name Q_tot ; Value { Integral { [ 0.5 * ep0*om0*Fabs[Im[CompXX[epsilonr[]]]]*(SquNorm[{u2d}+u1[]]) / (Pinc[]*d) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
EndIf
}
}
}
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" ] ;
For i In {0:2*nb_orders}
Print[ s_r~{i}[SurfCutSuper1], OnGlobal, Store i , Format Table , File > StrCat[myDir, Sprintf("temp_s_r_%g.txt", i-nb_orders)]];
Print[ s_t~{i}[SurfCutSubs1] , OnGlobal, Store (2*nb_orders+1+i), Format Table , File > StrCat[myDir, Sprintf("temp_s_t_%g.txt", i-nb_orders)]];
EndFor
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)]];
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_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"]];
Print[ Q_layer_dep[layer_dep] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_dep.txt"]];
Print[ Q_layer_cov[layer_cov] , OnGlobal, Format FrequencyTable, File > StrCat[myDir, "absorption-Q_layer_cov.txt"]];
If (flag_polar==1)
Print[ u_tot , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
EndIf
If (flag_polar==0)
Print[ u_tot , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Ez_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Ez_tot_%.2fnm.pos", lambda0/nm)];
// Print[ u, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u2d_%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u2d_%.2fnm.pos", lambda0/nm)];
EndIf
If(multiplot)
Echo[ Str["For i In {PostProcessing.NbViews-1:0:-1}",
" If(!StrCmp(View[i].Name, 'boundary') || !StrCmp(View[i].Name, 'boundary_Combine'))",
" Delete View[i];",
" EndIf",
"EndFor"], File StrCat[myDir,"tmp1.geo"]] ;
For i In {-nb_plot_periods:nb_plot_periods}
Print [ u_tot~{i}, OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("u_tot_lambda%.2fnm_%g.pos", lambda0/nm,i+nb_plot_periods)], ChangeOfCoordinates {$X+i*d,$Y,$Z}, Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm) ] ;
Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir,Sprintf("boundary_%g.pos",i+nb_plot_periods)], ChangeOfCoordinates {$X+i*d,$Y,$Z} ] ;
EndFor
Echo[ Str["Combine ElementsByViewName;",
"Hide {",
"Point{1,2,7,8,9,10,20,22};",
"Line{1,7,8,9,10,30,32,34,2,3,4,5,6,12,16,20,24,28};",
"Surface{36,48};}",
"Geometry.Color.Lines = {0,0,0};",
"l=PostProcessing.NbViews-1; View[l].ColorTable={Black}; ",
"View[l-1].Visible=1; View[l-1].ShowScale=0;",
"View[l].ShowScale=0; View[l].LineWidth=1.5; View[l].LineType=1;Geometry.LineWidth=0;"],
File StrCat[myDir,"tmp3.geo" ]] ;
EndIf
}
}
}
DefineConstant[
R_ = {"helmoltz_scalar", Name "GetDP/1ResolutionChoices", Visible 1},
C_ = {"-solve -pos -petsc_prealloc 100 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
P_ = {"postop_energy", Name "GetDP/2PostOperationChoices", Visible 1}
];
If(plotRTgraphs)
DefineConstant[
refl_ = {0, Name "GetDP/R0", ReadOnly 1, Graph "02000000", Visible 1},
abs_ = {0, Name "GetDP/total absorption", ReadOnly 1, Graph "00000002", Visible 1},
trans_ = {0, Name "GetDP/T0", ReadOnly 1, Graph "000000000002", Visible 1}
];
EndIf