Skip to content
Snippets Groups Projects
Commit 5b9e2008 authored by Guillaume Demesy's avatar Guillaume Demesy
Browse files

conical

parent 4cf738ee
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DiffractionGratings/grating2D.geo
+ 2
− 2
View file @ 5b9e2008
...@@ -16,11 +16,11 @@ ...@@ -16,11 +16,11 @@
// along with This program. If not, see <https://www.gnu.org/licenses/>. // along with This program. If not, see <https://www.gnu.org/licenses/>.
/** Select the config. (See also grating2D.pro) **/ /** Select the config. (See also grating2D.pro) **/
// Include "grating2D_data_AnisotropicGrating.geo"; Include "grating2D_data_AnisotropicGrating.geo";
// Include "grating2D_data_LamellarGrating.geo"; // Include "grating2D_data_LamellarGrating.geo";
// Include "grating2D_data_PhotonicCrystalSlab.geo"; // Include "grating2D_data_PhotonicCrystalSlab.geo";
// Include "grating2D_data_ResonantGrating.geo"; // Include "grating2D_data_ResonantGrating.geo";
Include "grating2D_data_plasmonics.geo"; // Include "grating2D_data_plasmonics.geo";
paramaille_rods = lambda_min*nm/(paramaille*paramaille_scale_rods); paramaille_rods = lambda_min*nm/(paramaille*paramaille_scale_rods);
lc_rod_out = lambda_min*nm/(paramaille*paramaille_scale_rod_out); lc_rod_out = lambda_min*nm/(paramaille*paramaille_scale_rod_out);
......
DiffractionGratings/grating2D_conical.pro 0 → 100644
+ 581
− 0
View file @ 5b9e2008
File added
DiffractionGratings/grating2D_data.geo 0 → 100644
+ 11
− 0
View file @ 5b9e2008
DefineConstant[
test_case = {"LamellarGrating",
Name "0Test case",
Choices {"AnisotropicGrating", "LamellarGrating", "PhotonicCrystalSlab", "ResonantGrating", "plasmonics"},
GmshOption "Reset", Autocheck 0
}
];
Include StrCat["grating2D_data_",test_case,".geo"];
DiffractionGratings/grating2D_scalar.pro 0 → 100644
+ 536
− 0
View file @ 5b9e2008
// Copyright (C) 2019 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/>.
/** Select the config. (See also grating2D.geo) **/
// Include "grating2D_data_AnisotropicGrating.geo";
// Include "grating2D_data_LamellarGrating.geo";
// Include "grating2D_data_PhotonicCrystalSlab.geo";
// Include "grating2D_data_ResonantGrating.geo";
// Include "grating2D_data_plasmonics.geo";
// Include "grating2D_materials.pro";
myDir = "res2D/";
DefineConstant[
lambda0 = {lambda_min , Min lambda_min, Max lambda_max, Step (lambda_max-lambda_min)/(nb_lambdas-1), Name StrCat[pp2, "0wavelength [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 {
If (flag_polar==1)
{ Name debr ; Value { Local { [ r[] ]; In Omega; Jacobian JVol; } } }
{ Name debt ; Value { Local { [ t[] ]; In Omega; Jacobian JVol; } } }
{ Name u ; Value { Local { [ {u2d} ]; In Omega; Jacobian JVol; } } }
{ Name epsr ; Value { Local { [ CompZZ[epsilonr[]] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_diff ; Value { Local { [ {u2d}+u1d[] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_tot ; Value { Local { [ {u2d}+u1[] ]; In Omega; Jacobian JVol; } } }
{ Name NormHz_tot ; Value { Local { [ Norm[{u2d}+u1[]] ]; In Omega; Jacobian JVol; } } }
{ Name E1 ; Value { Local { [ E1[] ]; In Omega; Jacobian JVol; } } }
{ Name u1 ; Value { Local { [ u1[] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totp1 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totp2 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totp3 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totp4 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totm1 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totm2 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totm3 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Hz_totm4 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name boundary ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
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
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 ; } } }
{ Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
EndIf
If (flag_polar==0)
{ Name u ; Value { Local { [ {u2d} ]; In Omega; Jacobian JVol; } } }
{ Name epsr ; Value { Local { [ CompZZ[epsilonr[]] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_diff ; Value { Local { [ {u2d}+u1d[] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_tot ; Value { Local { [ {u2d}+u1[] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totp1 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totp2 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totp3 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totp4 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]* 4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totm1 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-1*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totm2 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-2*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totm3 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-3*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name Ez_totm4 ; Value { Local { [ ({u2d}+u1[])*Exp[I[]*-4*CompX[k1[]]*d] ]; In Omega; Jacobian JVol; } } }
{ Name boundary ; Value { Local { [ bndCol[] ] ; In Plot_bnd ; Jacobian JVol ; } } }
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
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 ; } } }
{ Name lambda_step ; Value { Local { [ lambda0/nm ]; In Omega ; Jacobian JVol; } } }
EndIf
}
}
}
PostOperation {
{ Name postop_energy; NameOfPostProcessing postpro_energy ;
Operation {
If (flag_polar==1)
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"]];
Print[ Hz_tot , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("Hz_tot_lambda%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("Hz_tot_%.2fnm.pos", lambda0/nm)];
Print[ u, OnElementsOf Omega, File StrCat[myDir, Sprintf("u2d_%.2fnm_1.pos", lambda0/nm)] , Name Sprintf("u2d_%.2fnm.pos", lambda0/nm)];
EndIf
If (flag_polar==0)
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_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"] ];
Print[ Ez_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
}
}
}
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment