Select Git revision
topopt_conical.pro 29.17 KiB
/*
File : fictiveLamellarGratings_conical.pro
Strating from the geometry computed by fictiveLamellar_conical.geo, this program
uses the open source environment GetDP to solve the reflexion and transmission
problem on the sticks grating.
*/
Include "param.dat";
myDir_init = "res_conical/";
myDir = "../res_conical/";
lambda0 = lambda_wavelength * nm;
Group {
// Boundaries
surfLeft = Region[SURF_LEFT];
surfRight = Region[SURF_RIGHT];
surfUp = Region[SURF_UP];
surfDown = Region[SURF_DOWN];
surfCutSuperstrate = Region[SURF_INTEG_SUP];
surfCutSubstrate = Region[SURF_INTEG_SUB];
// Domains
pmlSup = Region[PML_SUP];
superstrate = Region[{SUPERSTRATE, SURF_INTEG_SUP}];
substrate = Region[{SUBSTRATE, SURF_INTEG_SUB}];
If (Flag_pml_inf == 1)
pmlInf = Region[PML_INF];
Else
pmlInf = Region[{}];
EndIf
design = Region[DESIGN];
Omega_source = Region[{design}];
Omega_nosource = Region[{pmlSup, pmlInf, superstrate, substrate}];
Omega = Region[{Omega_source, Omega_nosource}];
Omega_top = Region[{superstrate, design}];
Omega_bot = Region[{substrate}];
Omega_pml = Region[{pmlSup, pmlInf}];
Plot_domain = Region[{Omega, -pmlSup, -pmlInf}];
Plot_bnd = Region[SURF_PLOT];
Print_point = Region[PRINT_POINT];
Omega_DefTr = Region[{SUBSTRATE, SUPERSTRATE}];
Gama = Region[{surfCutSubstrate, surfCutSuperstrate}];
Tr = ElementsOf[Omega_DefTr, ConnectedTo Gama];
}
Function{
I[] = Complex[0.,1.];
EX[] = Vector[1,0,0];
EZ[] = Vector[0,0,1];
bndCol[Plot_bnd] = Complex[1.,1.];
om0 = 2.*Pi*cel/lambda0;
k0 = 2.*Pi/lambda0;
epsrSuper[] = Complex[epsilon_r_super_re, epsilon_r_super_im];
murSuper[] = Complex[mu_r_super_re, mu_r_super_im];
epsrMin[] = Complex[epsilon_r_min_re, epsilon_r_min_im];
epsrMax[] = Complex[epsilon_r_max_re, epsilon_r_max_im];
murVoxel[] = Complex[mu_r_voxel_re, mu_r_voxel_im];
epsrSubs[] = Complex[epsilon_r_subs_re, epsilon_r_subs_im];
murSubs[] = Complex[mu_r_subs_re, mu_r_subs_im];
k1norm[] = k0*n_super;
k2norm[] = k0*n_subs;
Ae = A;
theta_i = theta_incidence * deg2rad;
phi_i = phi_incidence * deg2rad; psi_i = psi_incidence * deg2rad;
alpha[] = -k0*n_super*Sin[theta_i]*Sin[phi_i];
beta1[] = -k0*n_super*Cos[theta_i];
gamma[] = -k0*n_super*Sin[theta_i]*Cos[phi_i];
beta2[] = -Sqrt[k0^2*epsrSubs[]-alpha[]^2-gamma[]^2];
k1[] = Vector[alpha[], beta1[], gamma[]];
k1r[] = Vector[alpha[],-beta1[], gamma[]];
k2[] = Vector[alpha[], beta2[], gamma[]];
rs[] = (beta1[]-beta2[])/(beta1[]+beta2[]);
ts[] = 2.*beta1[] /(beta1[]+beta2[]);
rp[] = (beta1[]*epsrSubs[]-beta2[]*epsrSuper[])/(beta1[]*epsrSubs[]+beta2[]*epsrSuper[]);
tp[] = (2.*beta1[]*epsrSubs[])/(beta1[]*epsrSubs[]+beta2[]*epsrSuper[]);
spol[] = Vector[-Cos[phi_i],0,Sin[phi_i]];
AmplEis[] = spol[];
AmplErs[] = rs[]*spol[];
AmplEts[] = ts[]*spol[];
AmplHis[] = Sqrt[ep0*epsrSuper[]/mu0] *spol[];
AmplHrs[] = Sqrt[ep0*epsrSuper[]/mu0]*rp[]*spol[];
AmplHts[] = Sqrt[ep0*epsrSuper[]/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*epsrSuper[]) * k1[] /\ His[];
Erp[] = -1/(om0*ep0*epsrSuper[]) * k1r[] /\ Hrs[];
Etp[] = -1/(om0*ep0*epsrSubs[]) * k2[] /\ Hts[];
BlochX_phase_shift[] = Exp[ I[]*alpha[]*period];
BlochX_phase_shift_adj[] = Exp[-I[]*alpha[]*period];
For i_order In {0:2*nb_orders}
alpha~{i_order}[] = CompX[k1[]] + 2*Pi/period * (i_order-nb_orders);
expmialphax~{i_order}[] = Exp[-I[]*alpha~{i_order}[]*X[]];
beta_super~{i_order}[] = -Sqrt[k0^2*epsrSuper[] - alpha~{i_order}[]^2 - gamma[]^2];
beta_subs~{i_order}[] = -Sqrt[k0^2*epsrSubs[] - alpha~{i_order}[]^2 - gamma[]^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[pmlSup] = epsilon_r_super_re * PML_Tensor[];
epsilonr[superstrate] = epsrSuper[] * TensorDiag[1,1,1];
If (Flag_design_metal == 0)
epsilonr[design] = ((epsrMax[] - epsrMin[]) * ScalarField[XYZ[]]{10000} + epsrMin[]) * TensorDiag[1,1,1];
Else
nDesign[design] = ((n_max - n_min) * ScalarField[XYZ[]]{10000} + n_min);
kDesign[design] = ((k_max - k_min) * ScalarField[XYZ[]]{10000} + k_min);
epsilonr[design] = (nDesign[] + I[]*kDesign[])^2 * TensorDiag[1,1,1];
EndIf
epsilonr[substrate] = epsrSubs[] * TensorDiag[1,1,1];
epsilonr[pmlInf] = epsilon_r_subs_re * PML_Tensor[];
epsilonr_annex[pmlSup] = epsilon_r_super_re * PML_Tensor[];
epsilonr_annex[superstrate] = epsrSuper[] * TensorDiag[1,1,1];
epsilonr_annex[Omega_source] = epsrSuper[] * TensorDiag[1,1,1];
epsilonr_annex[substrate] = epsrSubs[] * TensorDiag[1,1,1];
epsilonr_annex[pmlInf] = epsilon_r_subs_re * PML_Tensor[];
mur[pmlSup] = PML_Tensor[]; mur[superstrate] = TensorDiag[1,1,1];
mur[Omega_source] = TensorDiag[1,1,1];
mur[substrate] = TensorDiag[1,1,1];
mur[pmlInf] = PML_Tensor[];
Ei[] = Ae*(Cos[psi_i]*Eip[]-Sin[psi_i]*Eis[]);
Er[] = Ae*(Cos[psi_i]*Erp[]-Sin[psi_i]*Ers[]);
Et[] = Ae*(Cos[psi_i]*Etp[]-Sin[psi_i]*Ets[]);
Hi[] = 1/(om0*mu0*mur[]) * k1[] /\ Ei[];
Hr[] = 1/(om0*mu0*mur[]) * k1r[] /\ Er[];
Ht[] = 1/(om0*mu0*mur[]) * k2[] /\ Et[];
E1[pmlSup] = Ei[]+Er[];
E1[Omega_top] = Ei[]+Er[];
E1[Omega_bot] = Et[];
E1[pmlInf] = Et[];
E1d[pmlSup] = Er[];
E1d[Omega_top] = Er[];
E1d[Omega_bot] = Et[];
E1d[pmlInf] = Et[];
E1t[] = Vector[CompX[E1[]] , CompY[E1[]] , 0 ];
E1l[] = Vector[ 0 , 0 , CompZ[E1[]] ];
E1dt[] = Vector[CompX[E1d[]], CompY[E1d[]], 0 ];
E1dl[] = Vector[ 0 , 0 , CompZ[E1d[]]];
Etot_transfer[] = ComplexVectorField[XYZ[]]{1000};
comp_amp_deriv_x_r_transfer[] = ComplexScalarField[XYZ[]]{2000};
comp_amp_deriv_z_r_transfer[] = ComplexScalarField[XYZ[]]{2002};
H1[pmlSup] = Hi[]+Hr[];
H1[Omega_top] = Hi[]+Hr[];
H1[Omega_bot] = Ht[];
H1[pmlInf] = Ht[];
H1d[Omega_top] = Hr[];
H1d[Omega_bot] = Ht[];
Pinc[] = 0.5*Ae^2*Sqrt[epsrSuper[]*ep0/mu0] * Cos[theta_i];
}
Constraint {
{Name Bloch;
Case {
{ Region surfRight; Type LinkCplx ; RegionRef surfLeft;
Coefficient BlochX_phase_shift[]; Function Vector[$X-period,$Y,$Z] ;
}
}
}
{Name Bloch_adj;
Case {
{ Region surfRight; Type LinkCplx ; RegionRef surfLeft;
Coefficient BlochX_phase_shift_adj[]; Function Vector[$X-period,$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
{ Name GaussOnePoint ;
Case {
{ Type Gauss ;
Case {
{ GeoElement Point ; NumberOfPoints 1; }
{ GeoElement Line ; NumberOfPoints 1; }
{ GeoElement Triangle ; NumberOfPoints 1; }
}
}
}
}
}
FunctionSpace {
{ Name Hcurl; Type Form1;
BasisFunction {
{ Name se ; NameOfCoef he ; Function BF_Edge ; Support Omega; Entity EdgesOf[All]; }
{ Name se2; NameOfCoef he2; Function BF_Edge_2E ; Support Omega; Entity EdgesOf[All]; }
If (Flag_o2_interp==1)
{ Name se3; NameOfCoef he3; Function BF_Edge_3F_b; Support Omega; Entity FacetsOf[All]; }
{ Name se4; NameOfCoef he4; Function BF_Edge_3F_c; Support Omega; Entity FacetsOf[All]; }
{ Name se5; NameOfCoef he5; Function BF_Edge_4E ; Support Omega; Entity EdgesOf[All]; }
EndIf
}
Constraint {
{ NameOfCoef he ; EntityType EdgesOf ; NameOfConstraint Bloch; }
{ NameOfCoef he2; EntityType EdgesOf ; NameOfConstraint Bloch; }
If (Flag_o2_interp==1)
{ NameOfCoef he5; EntityType EdgesOf ; NameOfConstraint Bloch; }
EndIf
}
}
{ Name Hcurl_adj; Type Form1;
BasisFunction {
{ Name se ; NameOfCoef he ; Function BF_Edge ; Support Omega; Entity EdgesOf[All]; }
{ Name se2; NameOfCoef he2; Function BF_Edge_2E ; Support Omega; Entity EdgesOf[All]; }
If (Flag_o2_interp==1)
{ Name se3; NameOfCoef he3; Function BF_Edge_3F_b; Support Omega; Entity FacetsOf[All]; }
{ Name se4; NameOfCoef he4; Function BF_Edge_3F_c; Support Omega; Entity FacetsOf[All]; }
{ Name se5; NameOfCoef he5; Function BF_Edge_4E ; Support Omega; Entity EdgesOf[All]; }
EndIf
}
Constraint {
{ NameOfCoef he ; EntityType EdgesOf ; NameOfConstraint Bloch_adj; } { NameOfCoef he2; EntityType EdgesOf ; NameOfConstraint Bloch_adj; }
If (Flag_o2_interp==1)
{ NameOfCoef he5; EntityType EdgesOf ; NameOfConstraint Bloch_adj; }
EndIf
}
}
{ Name Hgrad_perp; Type Form1P;
BasisFunction {
{ Name sn ; NameOfCoef hn ; Function BF_PerpendicularEdge ; Support Omega; Entity NodesOf[All]; }
{ Name sn2; NameOfCoef hn2; Function BF_PerpendicularEdge_2E; Support Omega; Entity EdgesOf[All]; }
/* If (Flag_o2_interp==1)
{ Name sn3; NameOfCoef hn3; Function BF_PerpendicularEdge_2F; Support Omega; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef hn4; Function BF_PerpendicularEdge_3E; Support Omega; Entity EdgesOf[All]; }
{ Name sn5; NameOfCoef hn5; Function BF_PerpendicularEdge_3F; Support Omega; Entity FacetsOf[All]; }
EndIf */
// BF_PerpendicularEdge_2F, 1., QUA|HEX|PRI, 0 },
// BF_PerpendicularEdge_2V, 1., QUA|HEX, 0 },
// BF_PerpendicularEdge_3E, 2., ALL, 0 },
// BF_PerpendicularEdge_3F, 2., TRI|QUA|TET|HEX|PRI, 0 },
// BF_PerpendicularEdge_3V, 2., HEX|PRI, 0 },
}
Constraint {
{ NameOfCoef hn ; EntityType NodesOf ; NameOfConstraint Bloch; }
{ NameOfCoef hn2; EntityType EdgesOf ; NameOfConstraint Bloch; }
/* If (Flag_o2_interp==1)
{ NameOfCoef hn4; EntityType EdgesOf ; NameOfConstraint Bloch; }
EndIf */
}
}
{ Name Hgrad_perp_adj; Type Form1P;
BasisFunction {
{ Name sn ; NameOfCoef hn ; Function BF_PerpendicularEdge ; Support Omega; Entity NodesOf[All]; }
{ Name sn2; NameOfCoef hn2; Function BF_PerpendicularEdge_2E; Support Omega; Entity EdgesOf[All]; }
/* If (Flag_o2_interp==1)
{ Name sn3; NameOfCoef hn3; Function BF_PerpendicularEdge_2F; Support Omega; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef hn4; Function BF_PerpendicularEdge_3E; Support Omega; Entity EdgesOf[All]; }
{ Name sn5; NameOfCoef hn5; Function BF_PerpendicularEdge_3F; Support Omega; Entity FacetsOf[All]; }
EndIf */
// BF_PerpendicularEdge_2F, 1., QUA|HEX|PRI, 0 },
// BF_PerpendicularEdge_2V, 1., QUA|HEX, 0 },
// BF_PerpendicularEdge_3E, 2., ALL, 0 },
// BF_PerpendicularEdge_3F, 2., TRI|QUA|TET|HEX|PRI, 0 },
// BF_PerpendicularEdge_3V, 2., HEX|PRI, 0 },
}
Constraint {
{ NameOfCoef hn ; EntityType NodesOf ; NameOfConstraint Bloch_adj; }
{ NameOfCoef hn2; EntityType EdgesOf ; NameOfConstraint Bloch_adj; }
/* If (Flag_o2_interp==1)
{ NameOfCoef hn4; EntityType EdgesOf ; NameOfConstraint Bloch_adj; }
EndIf */
}
}
{ Name H_DensityField; Type Form3;
BasisFunction {
{ Name sxn ; NameOfCoef uxn ; Function BF_Volume ; Support design; Entity VolumesOf[All] ; }
}
}
}
Formulation {
// Trivial formulation to get the density field on each element of the mesh
{Name form_init_density; Type FemEquation;
Quantity {
{ Name rho; Type Local; NameOfSpace H_DensityField;}
}
Equation {
Galerkin { [ Dof{rho} , {rho} ]; In design; Jacobian JVol; Integration GaussOnePoint; }
Galerkin { [ 1 , {rho} ]; In design; Jacobian JVol; Integration GaussOnePoint; } }
}
// Trivial formulation to project the derivatives field on each element of the mesh
{Name form_project_adjoint_x; Type FemEquation;
Quantity {
{ Name drho_x; Type Local; NameOfSpace H_DensityField;}}
Equation {
Galerkin { [ Dof{drho_x} , {drho_x} ]; In design; Jacobian JVol; Integration GaussOnePoint; }
Galerkin { [ -comp_amp_deriv_x_r_transfer[] , {drho_x} ]; In design; Jacobian JVol; Integration GaussOnePoint; }
}
}
{Name form_project_adjoint_z; Type FemEquation;
Quantity {
{ Name drho_z; Type Local; NameOfSpace H_DensityField;}}
Equation {
Galerkin { [ Dof{drho_z} , {drho_z} ]; In design; Jacobian JVol; Integration GaussOnePoint; }
Galerkin { [ -comp_amp_deriv_z_r_transfer[] , {drho_z} ]; In design; Jacobian JVol; Integration GaussOnePoint; }
}
}
// {d El} = -EZ /\ grad El
{ Name helmholtz_conical; Type FemEquation;
Quantity {
{ Name Et; Type Local; NameOfSpace Hcurl; }
{ Name El; Type Local; NameOfSpace Hgrad_perp; }
}
Equation {
Galerkin {[ 1/mur[] * Dof{d Et} , {d Et} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ 1/mur[] * Dof{d El} , {d El} ] ; In Omega; Integration Int_1; Jacobian JVol; }
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; }
Galerkin {[ -k0^2 * (epsilonr[]-epsilonr_annex[]) * E1t[] , {Et} ] ; In Omega_source; Integration Int_1; Jacobian JVol; }
Galerkin {[ -k0^2 * (epsilonr[]-epsilonr_annex[]) * E1l[] , {El} ] ; In Omega_source; Integration Int_1; Jacobian JVol; }
}
}
{ Name helmholtz_conical_adj; Type FemEquation;
Quantity {
{ Name Et_adj; Type Local; NameOfSpace Hcurl_adj; }
{ Name El_adj; Type Local; NameOfSpace Hgrad_perp_adj; }
}
Equation {
Galerkin {[ 1/mur[] * Dof{d Et_adj} , {d Et_adj} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ 1/mur[] * Dof{d El_adj} , {d El_adj} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ Complex[0,gamma[]] / mur[] * Dof{d El_adj} , EZ[]*^{Et_adj} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[-Complex[0,gamma[]] / mur[] * (EZ[]*^Dof{Et_adj}) , {d El_adj} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ gamma[]^2 / mur[] * (EZ[]*^Dof{Et_adj}) , EZ[]*^{Et_adj} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ -k0^2 * epsilonr[] * Dof{Et_adj} , {Et_adj} ] ; In Omega; Integration Int_1; Jacobian JVol; }
Galerkin {[ -k0^2 * epsilonr[] * Dof{El_adj} , {El_adj} ] ; In Omega; Integration Int_1; Jacobian JVol; }
// Tertiary operator to solve the adjoint problem on x or z with SolveAgain
Galerkin {[ -expmialphax~{target_order}[]/period * ($ibasis == 0 ? EX[] : 0) , {Et_adj} ] ; In surfCutSuperstrate; Jacobian JSur; Integration Int_1; }
Galerkin {[ -expmialphax~{target_order}[]/period * ($ibasis == 0 ? EX[] : 0) , {El_adj} ] ; In surfCutSuperstrate; Jacobian JSur; Integration Int_1; }
Galerkin {[ -expmialphax~{target_order}[]/period * ($ibasis == 1 ? EZ[] : 0) , {Et_adj} ] ; In surfCutSuperstrate; Jacobian JSur; Integration Int_1; }
Galerkin {[ -expmialphax~{target_order}[]/period * ($ibasis == 1 ? EZ[] : 0) , {El_adj} ] ; In surfCutSuperstrate; Jacobian JSur; Integration Int_1; }
}
}
}
Resolution {
{ Name res_init_densities;
System {
{ Name D; NameOfFormulation form_init_density; }
} Operation {
CreateDir[Str[myDir_init]];
Generate[D]; Solve[D]; PostOperation[postop_init_density];
}
}
{ Name helmholtz_conical_direct_inverse;
System {
{ Name S ; NameOfFormulation helmholtz_conical ; Type ComplexValue; }
{ Name Sa ; NameOfFormulation helmholtz_conical_adj ; Type ComplexValue; }
{ Name Px ; NameOfFormulation form_project_adjoint_x ; Type ComplexValue; }
{ Name Pz ; NameOfFormulation form_project_adjoint_z ; Type ComplexValue; }
}
Operation {
GmshRead["rho_opt.pos",10000];
Generate[S] ; Solve[S] ; PostOperation[postop_energy];
Evaluate[$ibasis = 0]; // Adjoint problem on x
Generate[Sa] ; Solve[Sa] ; PostOperation[postop_adj_x];
Generate[Px] ; Solve[Px] ; PostOperation[postop_project_adjoint_x];
Evaluate[$ibasis = 1]; // Adjoint problem on z
GenerateRHS[Sa]; SolveAgain[Sa]; PostOperation[postop_adj_z];
Generate[Pz] ; Solve[Pz] ; PostOperation[postop_project_adjoint_z];
}
}
}
PostProcessing {
{ Name postpro_init_density; NameOfFormulation form_init_density; NameOfSystem D;
Quantity {
{ Name density; Value { Local { [ {rho} ]; In design; Jacobian JVol; } } }
}
}
{ Name postpro_project_adjoint_x; NameOfFormulation form_project_adjoint_x; NameOfSystem Px;
Quantity {
{ Name comp_amp_deriv_x_r_project; Value { Local { [ ElementVol[]*{drho_x} ]; In design; Jacobian JVol; } } }
}
}
{ Name postpro_project_adjoint_z; NameOfFormulation form_project_adjoint_z; NameOfSystem Pz;
Quantity {
{ Name comp_amp_deriv_z_r_project; Value { Local { [ ElementVol[]*{drho_z} ]; In design; Jacobian JVol; } } }
}
}
{ Name postpro_energy; NameOfFormulation helmholtz_conical; NameOfSystem S;
Quantity {
{ Name boundary ; Value { Local { [ bndCol[] ]; In Plot_bnd ; Jacobian JVol; } } }
//{ Name Etot_z ; Value { Local { [ CompZ[{Et}+{El}+E1[]]]; In Omega; Jacobian JVol; } } }
{ Name Etot ; Value { Local { [ {Et}+{El}+E1[] ]; In Omega; Jacobian JVol; } } }
//{ Name epsr ; Value { Local { [ CompZZ[epsilonr[]] ]; In Omega; Jacobian JVol; } } }
{ Name Q_tot ; Value { Integral{ [0.5*ep0*om0 * Fabs[Im[CompZZ[epsilonr[]]]] * ( SquNorm[{Et}+{El}+E1[]] ) / (Pinc[]*period) ] ; In Plot_domain ; Integration Int_1 ; Jacobian JVol ; } } }
For i_order In {0:2*nb_orders}
{ Name comp_amp_x_r~{i_order} ; Value{ Integral{ [ CompX[{Et}+E1dt[]] * expmialphax~{i_order}[]/period ] ; In surfCutSuperstrate ; Integration Int_1 ; Jacobian JSur ; } } }
{ Name comp_amp_z_r~{i_order} ; Value{ Integral{ [ CompZ[{El}+E1dl[]] * expmialphax~{i_order}[]/period ] ; In surfCutSuperstrate ; Integration Int_1 ; Jacobian JSur ; } } }
{ Name comp_amp_x_t~{i_order} ; Value{ Integral{ [ CompX[{Et}+E1t[]] * expmialphax~{i_order}[]/period ] ; In surfCutSubstrate ; Integration Int_1 ; Jacobian JSur ; } } }
{ Name comp_amp_z_t~{i_order} ; Value{ Integral{ [ CompZ[{El}+E1l[]] * expmialphax~{i_order}[]/period ] ; In surfCutSubstrate ; Integration Int_1 ; Jacobian JSur ; } } }
EndFor
For i_order In {0:2*nb_orders}
{ Name eff_r~{i_order} ; Value { Term{Type Global; [
1/(Ae^2*beta_super~{i_order}[]*beta1[]) * ((beta_super~{i_order}[]^2 + alpha~{i_order}[]^2) * SquNorm[$comp_ampx_r~{i_order}]
+ (beta_super~{i_order}[]^2 + gamma[]^2) * SquNorm[$comp_ampz_r~{i_order}]
+ 2*alpha~{i_order}[]*gamma[] * Re[$comp_ampx_r~{i_order}*Conj[$comp_ampz_r~{i_order}]]) ] ; In surfCutSuperstrate ; } } }
{ Name eff_t~{i_order} ; Value { Term{Type Global; [
1/(Ae^2*beta_subs~{i_order}[]*beta1[]) * ((beta_subs~{i_order}[]^2 + alpha~{i_order}[]^2) * SquNorm[$comp_ampx_t~{i_order}]
+ (beta_subs~{i_order}[]^2 + gamma[]^2) * SquNorm[$comp_ampz_t~{i_order}]
+ 2*alpha~{i_order}[]*gamma[] * Re[$comp_ampx_t~{i_order}*Conj[$comp_ampz_t~{i_order}]]) ] ; In surfCutSubstrate ; } } }
EndFor
} }
{ Name postpro_adj_x; NameOfFormulation helmholtz_conical_adj; NameOfSystem Sa;
Quantity {
//{ Name Eadj_x ; Value { Local { [ {Et_adj}+{El_adj} ] ; In Omega ; Jacobian JVol; } } }
//{ Name Eadj_projx ; Value { Local { [ CompX[{Et_adj}+{El_adj}] ] ; In Omega ; Jacobian JVol; } } }
If (Flag_design_metal == 0)
{ Name comp_amp_deriv_x_r ; Value { Local { [ k0^2 * (epsilon_r_max_re - epsilon_r_min_re) * Etot_transfer[] * ({Et_adj} + {El_adj}) ]; In design; Jacobian JVol; } } }
Else
{ Name comp_amp_deriv_x_r ; Value { Local { [ 2 * k0^2 * (n_max-n_min + I[]*(k_max-k_min))
* (n_min+I[]*k_min + ScalarField[XYZ[]]{10000} * (n_max-n_min + I[]*(k_max-k_min)))
* Etot_transfer[] * ({Et_adj} + {El_adj}) ]; In design; Jacobian JVol; } } }
EndIf
}
}
{ Name postpro_adj_z; NameOfFormulation helmholtz_conical_adj; NameOfSystem Sa;
Quantity {
//{ Name Eadj_z ; Value { Local { [ {Et_adj}+{El_adj} ] ; In Omega ; Jacobian JVol; } } }
//{ Name Eadj_projz ; Value { Local { [ CompZ[{Et_adj}+{El_adj}] ] ; In Omega ; Jacobian JVol; } } }
If (Flag_design_metal == 0)
{ Name comp_amp_deriv_z_r ; Value { Local { [ k0^2 * (epsilon_r_max_re - epsilon_r_min_re) * Etot_transfer[] * ({Et_adj} + {El_adj}) ]; In design; Jacobian JVol; } } }
Else
{ Name comp_amp_deriv_z_r ; Value { Local { [ 2 * k0^2 * (n_max-n_min + I[]*(k_max-k_min))
* (n_min+I[]*k_min + ScalarField[XYZ[]]{10000} * (n_max-n_min + I[]*(k_max-k_min)))
* Etot_transfer[] * ({Et_adj} + {El_adj}) ]; In design; Jacobian JVol; } } }
EndIf
}
}
}
PostOperation {
{ Name postop_init_density; NameOfPostProcessing postpro_init_density ;
Operation {
// template_density_table.pos is a template file that stores the position of each triangle of the mesh
Print[ density, OnElementsOf design, StoreInField (10000), File StrCat[myDir_init, "template_density_table.pos"]];
}
}
{ Name postop_project_adjoint_x; NameOfPostProcessing postpro_project_adjoint_x ;
Operation {
Print[ comp_amp_deriv_x_r_project, OnElementsOf design, File StrCat[myDir, Sprintf("comp_amp_deriv_x_r_project_%g.pos", i_lambda)]];
}
}
{ Name postop_project_adjoint_z; NameOfPostProcessing postpro_project_adjoint_z ;
Operation {
Print[ comp_amp_deriv_z_r_project, OnElementsOf design, File StrCat[myDir, Sprintf("comp_amp_deriv_z_r_project_%g.pos", i_lambda)]];
}
}
{ Name postop_energy; NameOfPostProcessing postpro_energy ;
Operation {
Print[ boundary, OnElementsOf Plot_bnd, File StrCat[myDir, "boundary.pos"] ] ;
//Print[ Etot, OnGrid {$A, $B, 0} {-period/2:period/2:period/10, -h_subs:h_super:(h_super+h_subs)/30, 0}, File StrCat[myDir, "Etot.pos"] ];
//Print[ Etot_z, OnElementsOf Omega, File StrCat[myDir, Sprintf("Etot_z_%g.pos", i_lambda)]];
Print[ Etot , OnElementsOf Plot_domain, StoreInField (1000), File StrCat[myDir, Sprintf("Etot_%g.pos", i_lambda)]];
//Print[ epsr , OnElementsOf Plot_domain, File StrCat[myDir, Sprintf("epsr_%g.pos", i_lambda)]];
Print[ Q_tot[Plot_domain], OnGlobal, Format FrequencyTable, File StrCat[myDir, Sprintf("absorption-Q_tot_%g.txt", i_lambda)]];
For i_order In {0:2*nb_orders}
Print[ comp_amp_x_r~{i_order}[surfCutSuperstrate], OnGlobal, StoreInVariable $comp_ampx_r~{i_order}, Format Table, File StrCat[myDir, Sprintf("comp_amp_x_r_%g_%g.txt",i_order-nb_orders,i_lambda)]];
Print[ comp_amp_z_r~{i_order}[surfCutSuperstrate], OnGlobal, StoreInVariable $comp_ampz_r~{i_order}, Format Table, File StrCat[myDir, Sprintf("comp_amp_z_r_%g_%g.txt",i_order-nb_orders,i_lambda)]];
Print[ comp_amp_x_t~{i_order}[surfCutSuperstrate], OnGlobal, StoreInVariable $comp_ampx_t~{i_order}, Format Table, File StrCat[myDir, Sprintf("comp_amp_x_t_%g_%g.txt",i_order-nb_orders,i_lambda)]];
Print[ comp_amp_z_t~{i_order}[surfCutSuperstrate], OnGlobal, StoreInVariable $comp_ampz_t~{i_order}, Format Table, File StrCat[myDir, Sprintf("comp_amp_z_t_%g_%g.txt",i_order-nb_orders,i_lambda)]];
EndFor
For i_order In {0:2*nb_orders}
Print[ eff_r~{i_order}[surfCutSuperstrate], OnRegion surfCutSuperstrate, Format Table, File StrCat[myDir, Sprintf("efficiency_r_%g_%g.txt", i_order-nb_orders,i_lambda)]];
Print[ eff_r~{i_order}[surfCutSuperstrate], OnRegion surfCutSuperstrate, Format Table, File > StrCat[myDir, Sprintf("efficiency_r_%g_%g_list.txt", i_order-nb_orders,i_lambda)]];
Print[ eff_t~{i_order}[surfCutSuperstrate], OnRegion surfCutSuperstrate, Format Table, File StrCat[myDir, Sprintf("efficiency_t_%g_%g.txt", i_order-nb_orders,i_lambda)]];
EndFor
}
}
{ Name postop_adj_x; NameOfPostProcessing postpro_adj_x ; Operation {
//Print[ Eadj_x, OnGrid {$A, $B, 0} {-period/2:period/2:period/10, -h_subs:h_super:(h_super+h_subs)/30, 0}, File StrCat[myDir, "Eadj_x.pos"] ];
//Print[ Eadj_projx, OnElementsOf Plot_domain, File StrCat[myDir, "Eadj_projx.pos"] ];
//Print[ Eadj_x , OnElementsOf Plot_domain, File StrCat[myDir, "Eadj_x.pos"] ];
Print[ comp_amp_deriv_x_r, OnElementsOf design, StoreInField (2000), File StrCat[myDir, Sprintf("comp_amp_deriv_x_r_%g.pos",i_lambda)], Depth 1];
}
}
{ Name postop_adj_z; NameOfPostProcessing postpro_adj_z ;
Operation {
//Print[ Eadj_z, OnGrid {$A, $B, 0} {-period/2:period/2:period/10, -h_subs:h_super:(h_super+h_subs)/30, 0}, File StrCat[myDir, "Eadj_z.pos"] ];
//Print[ Eadj_projz, OnElementsOf Plot_domain, File StrCat[myDir, "Eadj_projz.pos"] ];
//Print[ Eadj_z , OnElementsOf Plot_domain, File StrCat[myDir, "Eadj_z.pos"] ];
Print[ comp_amp_deriv_z_r, OnElementsOf design, StoreInField (2002), File StrCat[myDir, Sprintf("comp_amp_deriv_z_r_%g.pos",i_lambda)], Depth 1];
}
}
}