diff --git a/ElectromagneticScattering/scattererTmatrix.geo b/ElectromagneticScattering/scattererTmatrix.geo
index 01c68c93867729f200a00f3e4258f45d23584907..422a50ce9d8a761ee88a5d74bbee5ab45ac5d35b 100644
--- a/ElectromagneticScattering/scattererTmatrix.geo
+++ b/ElectromagneticScattering/scattererTmatrix.geo
@@ -74,3 +74,9 @@ Mesh.VolumeEdges = 0;
Mesh.ElementOrder = 2;
Mesh.HighOrderOptimize = 1;
+//+
+Hide "*";
+//+
+Show {
+ Point{1}; Point{2}; Curve{1}; Curve{2}; Curve{3}; Surface{1}; Surface{2}; Volume{1};
+}
diff --git a/QuasiNormalModeExpansion/QNM_expansion_DtN.geo b/QuasiNormalModeExpansion/QNM_expansion_DtN.geo
new file mode 100644
index 0000000000000000000000000000000000000000..619575eaff68155d06eea11ef30543094621b4f1
--- /dev/null
+++ b/QuasiNormalModeExpansion/QNM_expansion_DtN.geo
@@ -0,0 +1,51 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+Include "QNM_expansion_DtN_data.geo";
+SetFactory("OpenCASCADE");
+lc_bg = lambda0/paramaille;
+lc_sc = lambda0/(paramaille*3);
+
+Rectangle(1) = {-dx_scat/2, -dy_scat/2, 0, dx_scat, dy_scat, 0};
+Disk(2) = {0, 0, 0, r_domain, r_domain};
+If (flag_PML==1)
+ Disk(3) = {0, 0, 0, r_pml, r_pml};
+EndIf
+Coherence;
+
+Physical Point("PrintPoint",100) = {1};
+
+Physical Surface("SCAT",1) = {1};
+If (flag_PML==1)
+ Physical Surface("BG",2) = {2};
+ Characteristic Length {6} = lc_bg;
+ Physical Surface("PML",3) = {3};
+ Physical Surface("SRC",4) = {};
+ Characteristic Length {2,3,4,5} = lc_sc;
+ Characteristic Length {1} = lc_bg;
+ Characteristic Length {7} = lc_bg;
+ Physical Line("SURF_SM",10) = {};
+ Else
+ Physical Surface("BG",2) = {2};
+ Physical Surface("SRC",4) = {};
+ Characteristic Length {2,3,4,5} = lc_sc;
+ Characteristic Length {1} = lc_bg;
+ Physical Line("SURF_SM",10) = {1};
+EndIf
+
+Rectangle(10) = {-dx_scat/2+domain_tot*yshift, -dy_scat/2, 0, dx_scat, dy_scat, 0};
+Disk(20) = {domain_tot*yshift, 0, 0, r_domain, r_domain};
+Delete{Surface{10,20};}
+If (flag_PML==1)
+ Disk(30) = {domain_tot*yshift, 0, 0, r_pml, r_pml};
+ Delete{Surface{30};}
+EndIf
+
+If (flag_PML==1)
+ Point(newp) = {-dx_scat/2+domain_tot*yshift, 2*r_pml, 0};
+Else
+ Point(newp) = {-dx_scat/2+domain_tot*yshift, 2*r_domain, 0};
+EndIf
+Geometry.Points = 0;
+
diff --git a/QuasiNormalModeExpansion/QNM_expansion_DtN.pro b/QuasiNormalModeExpansion/QNM_expansion_DtN.pro
new file mode 100644
index 0000000000000000000000000000000000000000..8ca575f0c1f3ccfe5d7dfcf6f2041b90449e2c57
--- /dev/null
+++ b/QuasiNormalModeExpansion/QNM_expansion_DtN.pro
@@ -0,0 +1,561 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+Include "QNM_expansion_DtN_data.geo";
+stored_exp_field_AR_PML = 10000;
+stored_exp_field_AR_mur = 20000;
+stored_exp_field_AR_BT2 = 30000;
+stored_exp_field_manual_PML = 40000;
+CreateDir["out"];
+
+Group {
+ // Domains
+ Scat = Region[1];
+ Bg = Region[2];
+ Src = Region[4];
+ If (flag_PML==0)
+ PML = Region[{}];
+ Surf = Region[10];
+ Else
+ PML = Region[{3}];
+ Surf = Region[{10}];
+ EndIf
+ Omega = Region[{Bg,Scat,PML,Src}];
+ PrintPoint = Region[100];
+}
+
+Function{
+ I[] = Complex[0.0,1.0];
+
+ k0[] = $realomega/cel;
+ omega0[] = k0[]*cel;
+
+ spml_re = 1;
+ spml_im = 1;
+
+ R_pml_in = r_domain;
+ sr[] = Complex[spml_re,spml_im];
+ rho[] = Sqrt[X[]*X[]+Y[]*Y[]];// Norm[ XYZ[ ] ];
+ phi[] = Atan2[ Y[] , X[] ];
+ rr[] = R_pml_in + (rho[] - R_pml_in) * sr[];
+ drr[] = sr[];
+ pmlcirctens[] = Rotate[TensorDiag[rr[]/(drr[]*rho[]), (drr[]*rho[])/rr[], (drr[]*rr[])/rho[]],0,0,-phi[]];
+
+ epsr[Scat] = Complex[epsr_scat_re,epsr_scat_im];
+ epsr[Bg] = epsr_bg;
+ epsr[Src] = epsr_bg;
+ epsr[PML] = pmlcirctens[];
+
+ epsr_annex[Scat] = epsr_bg;
+ epsr_annex[Bg] = epsr_bg;
+ epsr_annex[Src] = epsr_bg;
+ epsr_annex[PML] = pmlcirctens[];
+
+ mur[Scat] = 1;
+ mur[Bg] = 1;
+ mur[Src] = 1;
+ mur[PML] = pmlcirctens[];
+
+ alpha[] = k0[]*Cos[theta0];
+ beta[] = k0[]*Sin[theta0];
+ siwt = 1;
+ ui[] = Vector[0, 0, Exp[I[]*siwt*(alpha[]*X[]+beta[]*Y[])] ];
+ Hi[] = I[]/(omega0[]*mu0) * Vector[alpha[],beta[],0] /\ ui[];
+ Pinc[] = 0.5*1^2*Sqrt[eps0*epsr_bg/mu0] * 2*r_domain;
+
+ THETA[] = Atan2[Y[],X[]];
+ TG[] = Vector[-Sin[THETA[]], Cos[THETA[]],0];
+
+ source_PW[] = k0[]^2 * (epsr_annex[]-epsr[]) * ui[];
+
+ For k In {0:#tabrealomegasexpansion()-1}
+ k0~{k}[] = tabrealomegasexpansion(k)/cel;
+ alpha~{k}[] = k0~{k}[]*Cos[theta0];
+ beta~{k}[] = k0~{k}[]*Sin[theta0];
+ ui~{k}[] = Vector[0, 0, Exp[I[]*siwt*(alpha~{k}[]*X[]+beta~{k}[]*Y[])] ];
+ source_PW~{k}[] = k0~{k}[]^2 * (epsr_annex[]-epsr[]) * ui~{k}[];
+ us_exp_PML_modal_manual~{k}[] = ComplexScalarField[XYZ[]]{stored_exp_field_manual_PML+k};
+ us_exp_PML_modal~{k}[] = ComplexScalarField[XYZ[]]{stored_exp_field_AR_PML+k};
+ us_exp_mur_modal~{k}[] = ComplexScalarField[XYZ[]]{stored_exp_field_AR_mur+k};
+ us_exp_BT2_modal~{k}[] = ComplexScalarField[XYZ[]]{stored_exp_field_AR_BT2+k};
+ EndFor
+
+
+ source_G[Src] = Vector[0, 0, 1];
+ source_G[Scat] = 0;
+ source_G[Bg] = 0;
+ source_G[PML] = 0;
+
+ Rout = r_domain;
+
+ alphaBT[] = 1/(2*Rout) - I[]/(8*k0[]*Rout^2*(1+I[]/(k0[]*Rout)));
+ betaBT[] = - 1/(2*I[]*k0[]*(1+I[]/(k0[]*Rout)));
+}
+
+Constraint {
+}
+
+Jacobian {
+ { Name JVol ;
+ Case {
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name JSur ;
+ Case {
+ { Region All ; Jacobian Sur ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int_1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 12 ; }
+ }
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { Name Hgrad_P; Type Form1P;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_PerpendicularEdge_1N; Support Region[{Omega,Surf}]; Entity NodesOf[Omega]; }
+ { Name sn2; NameOfCoef un2; Function BF_PerpendicularEdge_2E; Support Region[{Omega,Surf}]; Entity EdgesOf[Omega]; }
+ }
+ Constraint {
+ }
+ }
+}
+
+Formulation {
+ {Name form_direct_PML; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ Galerkin{[ k0[]^2*epsr[]*Dof{u} , {u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ -1/mur[]*Dof{Curl u} , {Curl u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ -source_PW[] , {u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ }
+ }
+
+ {Name form_direct_mur; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ Galerkin{[ k0[]^2*epsr[]*Dof{u} , {u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ -1/mur[]*Dof{Curl u} , {Curl u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ -source_PW[] , {u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ I[]*k0[]*Dof{u} , {u}]; In Surf; Jacobian JSur; Integration Int_1;}
+ }
+ }
+
+ {Name form_direct_BT2; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ Galerkin{[ k0[]^2*epsr[]*Dof{u} , {u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ -1/mur[]*Dof{Curl u} , {Curl u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ -source_PW[] , {u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{[ I[]*k0[]*Dof{u} , {u}]; In Surf; Jacobian JSur; Integration Int_1;}
+ Galerkin{[ -alphaBT[] * Dof{u} , {u}]; In Surf; Jacobian JSur; Integration Int_1;}
+ Galerkin{[ -betaBT[] * Dof{Curl u} , {Curl u}]; In Surf; Jacobian JSur; Integration Int_1;}
+ }
+ }
+
+ {Name form_modal_linear_PML; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ Galerkin{Eig[-1/cel^2*epsr[]*Dof{u} , {u}]; Order 2 ; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{ [ -1/mur[]*Dof{Curl u} , {Curl u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ }
+ }
+
+ {Name form_modal_PML; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ Galerkin{Eig[-1/cel^2*epsr[]*Dof{u} , {u}]; Rational 1 ; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{Eig[ -1/mur[]*Dof{Curl u} , {Curl u}]; Rational 2 ; In Omega; Jacobian JVol; Integration Int_1;}
+ }
+ }
+ {Name form_modal_mur_poly; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ Galerkin{ Eig[-1/cel^2*epsr[]*Dof{u} , {u}]; Order 2 ; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{ [-1/mur[]*Dof{Curl u} , {Curl u}]; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{ Eig[ 1/cel*Dof{u} , {u}]; Order 1 ; In Surf ; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ {Name form_modal_mur; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ Galerkin{ Eig[-1/cel^2*epsr[]*Dof{u} , {u}]; Rational 1 ; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{ Eig[-1/mur[]*Dof{Curl u} , {Curl u}]; Rational 2 ; In Omega; Jacobian JVol; Integration Int_1;}
+ Galerkin{ Eig[ 1/cel*Dof{u} , {u}]; Rational 3 ; In Surf ; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ {Name form_modal_BT2; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ // {-1/cel^2,0,0}/{1} --------- (i\omega)^2/c^2
+ Galerkin{Eig[ epsr[]* Dof{u} , {u}]; Rational 1; In Omega; Jacobian JVol; Integration Int_1;}
+ // {1}/{1} 1
+ Galerkin{Eig[ -1/mur[]*Dof{Curl u} , {Curl u}]; Rational 2; In Omega; Jacobian JVol; Integration Int_1;}
+ // {1/cel,0}/{1} --------- i\omega/c
+ Galerkin{Eig[ Dof{u} , {u}]; Rational 3; In Surf ; Jacobian JSur; Integration Int_1;}
+ // {-4*Rout/cel,3}/{-8*Rout^2/cel,8*Rout} --------- \frac{-4Ri\omega/c+3}{-8R^2i\omega/c+8R}
+ Galerkin{Eig[ Dof{u} , {u}]; Rational 4; In Surf ; Jacobian JSur; Integration Int_1;}
+ // {-Rout}/{2*Rout/cel,-2} --------- \frac{-R}{2Ri\omega/c-2}
+ Galerkin{Eig[ Dof{Curl u} , {Curl u}]; Rational 5; In Surf ; Jacobian JSur; Integration Int_1;}
+ }
+ }
+ {Name form_RHS; Type FemEquation;
+ Quantity{{ Name u; Type Local; NameOfSpace Hgrad_P;}}
+ Equation{
+ // Galerkin { [ 1. * Dof{u} , { u} ]; In Omega ; Jacobian JSur; Integration Int_1;}
+ Galerkin{[-source_PW[] , {u}]; In Omega; Jacobian JVol; Integration Int_1;}
+
+ }
+ }
+}
+
+Resolution {
+ { Name res_modal_PML;
+ System {
+ { Name M1; NameOfFormulation form_modal_linear_PML; Type ComplexValue; }
+ { Name M2; NameOfFormulation form_direct_PML ; Type ComplexValue; }
+ }
+ Operation {
+ GenerateSeparate[M1] ;
+ EigenSolve[M1,neig_PML,target_re_PML,target_im_PML];
+ PostOperation[postop_modal_PML];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ Generate[M2];
+ Solve[M2];
+ PostOperation[postop_direct_PML];
+ EndFor
+ For k In {0:#tabrealomegasexpansion()-1}
+ PostOperation[postop_gmsh~{k}];
+ Sleep[0.3];
+ EndFor
+ }
+ }
+
+ { Name res_modal_AR_BT2;
+ System{
+ { Name V ; NameOfFormulation form_RHS ; Type ComplexValue; }
+ { Name M1; NameOfFormulation form_modal_BT2 ; Type ComplexValue; }
+ { Name M2; NameOfFormulation form_direct_BT2 ; Type ComplexValue; }
+
+ }
+ Operation{
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ GenerateListOfRHS[V,Omega,#tabrealomegasexpansion()];
+ EndFor
+ GenerateSeparate[M1];
+ EigenSolveAndExpand[M1,neig_BT,target_re_BT,target_im_BT,
+ {{-1/cel^2,0,0},{1},{1/cel,0},{ 4*Rout/cel,-3} ,{Rout}},
+ { {1},{1}, {1},{-8*Rout^2/cel,8*Rout},{2*Rout/cel,-2}},
+ tabrealomegasexpansion(),
+ V];
+ PostOperation[postop_modal_AR_BT2];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ Generate[M2];
+ Solve[M2];
+ PostOperation[postop_direct_BT2];
+ EndFor
+ For k In {0:#tabrealomegasexpansion()-1}
+ PostOperation[postop_gmsh~{k}];
+ Sleep[0.5];
+ EndFor
+ }
+ }
+
+ { Name res_modal_AR_PML;
+ System{
+ { Name M1; NameOfFormulation form_modal_PML ; Type ComplexValue; }
+ { Name V ; NameOfFormulation form_RHS ; Type ComplexValue; }
+ }
+ Operation{
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ GenerateListOfRHS[V,Omega,#tabrealomegasexpansion()];
+ EndFor
+ GenerateSeparate[M1];
+ EigenSolveAndExpand[M1,neig_PML,target_re_BT,target_im_BT,
+ {{1,0,0},{1}},
+ { {1},{1}},
+ tabrealomegasexpansion(),
+ V];
+ }
+ }
+ { Name res_modal_AR_mur;
+ System{
+ { Name V; NameOfFormulation form_RHS ; Type ComplexValue; }
+ { Name M1; NameOfFormulation form_modal_mur ; Type ComplexValue; }
+ { Name M2; NameOfFormulation form_direct_mur ; Type ComplexValue; }
+ }
+ Operation{
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ GenerateListOfRHS[V,Omega,#tabrealomegasexpansion()];
+ EndFor
+ GenerateSeparate[M1];
+ EigenSolveAndExpand[M1,neig_BT,target_re_BT,target_im_BT,
+ {{1,0,0},{1},{1,0}},
+ { {1},{1}, {1}},
+ tabrealomegasexpansion(),
+ V];
+ PostOperation[postop_modal_AR_mur];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ Generate[M2];
+ Solve[M2];
+ PostOperation[postop_direct_mur];
+ EndFor
+ For k In {0:#tabrealomegasexpansion()-2}
+ PostOperation[postop_gmsh~{k}];
+ Sleep[0.5];
+ EndFor
+ }
+ }
+}
+
+PostProcessing {
+ { Name postpro_direct_PML; NameOfFormulation form_direct_PML;
+ Quantity {
+ { Name us ; Value { Local{ [CompZ[{u} ]]; In Omega; Jacobian JVol; } } }
+ { Name ut ; Value { Local{ [CompZ[{u}+ui[]]]; In Omega; Jacobian JVol; } } }
+ { Name intnorm2utot ; Value { Integral { [SquNorm[{u}+ui[]]]; In Scat ; Jacobian JVol; Integration Int_1 ; } } }
+ }
+ }
+ { Name postpro_direct_mur; NameOfFormulation form_direct_mur;
+ Quantity {
+ { Name us ; Value { Local{ [CompZ[{u} ]]; In Omega; Jacobian JVol; } } }
+ { Name ut ; Value { Local{ [CompZ[{u}+ui[]]]; In Omega; Jacobian JVol; } } }
+ { Name intnorm2utot ; Value { Integral { [SquNorm[{u}+ui[]]]; In Scat ; Jacobian JVol; Integration Int_1 ; } } }
+ }
+ }
+ { Name postpro_direct_BT2; NameOfFormulation form_direct_BT2;
+ Quantity {
+ { Name us ; Value { Local{ [CompZ[{u} ]]; In Omega; Jacobian JVol; } } }
+ { Name ut ; Value { Local{ [CompZ[{u}+ui[]]]; In Omega; Jacobian JVol; } } }
+ { Name intnorm2utot ; Value { Integral { [SquNorm[{u}+ui[]]]; In Scat ; Jacobian JVol; Integration Int_1 ; } } }
+ }
+ }
+
+ // ANCHOR POSTPRO_modal_PML
+ { Name postpro_modal_PML; NameOfFormulation form_modal_PML;
+ Quantity {
+ { Name u ; Value { Local{ [CompZ[{u}] ] ; In Omega ; Jacobian JVol; } } }
+ { Name normu ; Value { Local{ [ Norm[{u}] ] ; In Omega ; Jacobian JVol; } } }
+ { Name ev_re ; Value { Local{ [$EigenvalueReal] ; In PrintPoint; Jacobian JVol; } } }
+ { Name ev_im ; Value { Local{ [$EigenvalueImag] ; In PrintPoint; Jacobian JVol; } } }
+ { Name vol ; Value { Integral { [1] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_re ; Value { Integral { [$EigenvalueReal/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_im ; Value { Integral { [$EigenvalueImag/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name ev ; Value { Integral { [Complex[$EigenvalueReal,$EigenvalueImag]/$vol] ; In PrintPoint; Jacobian JVol; Integration Int_1 ; } } }
+ For i In {0:#tabrealomegasexpansion()-1}
+ { Name overlap~{i} ; Value { Integral { [{u}*source_PW~{i}[]] ; In Scat ; Jacobian JVol; Integration Int_1 ; } } }
+ EndFor
+ { Name den ; Value { Integral { [{u}*epsr[]*{u}] ; In Omega ; Jacobian JVol; Integration Int_1 ; } } }
+ For i In {0:#tabrealomegasexpansion()-1}
+ { Name us_PML_modal_manual~{i};
+ Value {
+ For j In {0:(neig_PML-1)}
+ // cel^2/(w^2-w_k^2)
+ Local { [ CompZ[cel^2/(tabrealomegasexpansion(i)^2-$ev~{j}^2) * $overlap~{i}~{j}/$den~{j} * {u}[neig_PML-1-j]]] ; In Omega; Jacobian JVol; }
+ // // cel^2/(2*w*(w-w_k))
+ // Local { [ CompZ[cel^2/(2*tabrealomegasexpansion(i)*(tabrealomegasexpansion(i)-$ev~{j})) * $overlap~{i}~{j}/$den~{j} * {u}[neig_PML-1-j]]] ; In Omega; Jacobian JVol; }
+ // // NEP style : cel^2/(2*w_k*(w-w_k))
+ // Local { [ CompZ[cel^2/(2*$ev~{j}*(tabrealomegasexpansion(i)-$ev~{j})) * $overlap~{i}~{j}/$den~{j} * {u}[neig_PML-1-j]]] ; In Omega; Jacobian JVol; }
+ EndFor
+ }
+ }
+ { Name us_PML_modal_manual_stored~{i}; Value {Local { [us_exp_PML_modal_manual~{i}[]] ; In Omega; Jacobian JVol; } } }
+ EndFor
+ For i In {0:#tabrealomegasexpansion()-1}
+ { Name intnorm2utot~{i} ; Value { Integral { [SquNorm[us_exp_PML_modal_manual~{i}[] + CompZ[ui~{i}[]]]]; In Scat; Jacobian JVol; Integration Int_1 ; } } }
+ EndFor
+ }
+ }
+
+ { Name postpro_modal_AR_BT2; NameOfFormulation form_modal_BT2;
+ Quantity {
+ { Name u ; Value { Local{ [CompZ[{u}] ] ; In Omega ; Jacobian JVol; } } }
+ { Name normu ; Value { Local{ [ Norm[{u}] ] ; In Omega ; Jacobian JVol; } } }
+ { Name ev_re ; Value { Local{ [$EigenvalueReal] ; In PrintPoint; Jacobian JVol; } } }
+ { Name ev_im ; Value { Local{ [$EigenvalueImag] ; In PrintPoint; Jacobian JVol; } } }
+ { Name vol ; Value { Integral { [1] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_re ; Value { Integral { [$EigenvalueReal/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_im ; Value { Integral { [$EigenvalueImag/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ For i In {0:#tabrealomegasexpansion()-1}
+ { Name intnorm2utot~{i} ; Value { Integral { [SquNorm[us_exp_BT2_modal~{i}[] + CompZ[ui~{i}[]]]] ; In Scat ; Jacobian JVol; Integration Int_1 ; } } }
+ EndFor
+ }
+ }
+ { Name postpro_modal_AR_PML; NameOfFormulation form_modal_PML;
+ Quantity {
+ { Name u ; Value { Local{ [CompZ[{u}] ] ; In Omega ; Jacobian JVol; } } }
+ { Name ev_re ; Value { Local{ [$EigenvalueReal] ; In PrintPoint; Jacobian JVol; } } }
+ { Name ev_im ; Value { Local{ [$EigenvalueImag] ; In PrintPoint; Jacobian JVol; } } }
+ For i In {0:#tabrealomegasexpansion()-1}
+ { Name intnorm2utot~{i} ; Value { Integral { [SquNorm[us_exp_PML_modal~{i}[] + CompZ[ui~{i}[]]]] ; In Scat ; Jacobian JVol; Integration Int_1 ; } } }
+ EndFor
+ }
+ }
+ { Name postpro_modal_AR_mur; NameOfFormulation form_modal_mur;
+ Quantity {
+ { Name u ; Value { Local{ [CompZ[{u}] ] ; In Omega ; Jacobian JVol; } } }
+ { Name ev_re ; Value { Local{ [$EigenvalueReal] ; In PrintPoint; Jacobian JVol; } } }
+ { Name ev_im ; Value { Local{ [$EigenvalueImag] ; In PrintPoint; Jacobian JVol; } } }
+ { Name vol ; Value { Integral { [1] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_re ; Value { Integral { [$EigenvalueReal/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_im ; Value { Integral { [$EigenvalueImag/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ For i In {0:#tabrealomegasexpansion()-1}
+ { Name intnorm2utot~{i} ; Value { Integral { [SquNorm[us_exp_mur_modal~{i}[] + CompZ[ui~{i}[]]]] ; In Scat ; Jacobian JVol; Integration Int_1 ; } } }
+ EndFor
+ }
+ }
+}
+
+PostOperation {
+ { Name postop_direct_PML; NameOfPostProcessing postpro_direct_PML ;
+ Operation {
+ Print[ us , OnElementsOf Omega, File "out/us_PML_direct.pos", ChangeOfCoordinates {$X+domain_tot*yshift,$Y,$Z}, Name "us_PML_direct"];
+ Print[ intnorm2utot[Scat] , OnGlobal, File "out/ACC_PML_direct.out", Format SimpleTable];
+ If(flag_PRINTFIELDS==1)
+ Print[ us , OnElementsOf Omega, File "out/us_PML.pos", Name "us_PML"];
+ EndIf
+ }
+ }
+ { Name postop_direct_mur; NameOfPostProcessing postpro_direct_mur ;
+ Operation {
+ Print[ us , OnElementsOf Omega, File "out/us_mur_direct.pos", ChangeOfCoordinates {$X+domain_tot*yshift,$Y,$Z} , Name "us_mur_direct"];
+ Print[ intnorm2utot[Scat] , OnGlobal, File "out/ACC_mur_direct.out", Format SimpleTable];
+ If(flag_PRINTFIELDS==1)
+ Print[ us , OnElementsOf Omega, File "out/us_mur.pos", Name "us_mur"];
+ EndIf
+ }
+ }
+ { Name postop_direct_BT2; NameOfPostProcessing postpro_direct_BT2 ;
+ Operation {
+ Print[ us , OnElementsOf Omega, File "out/us_BT2_direct.pos", ChangeOfCoordinates {$X+domain_tot*yshift,$Y,$Z} , Name "us_BT2_direct"];
+ Print[ intnorm2utot[Scat] , OnGlobal, File "out/ACC_BT2_direct.out", Format SimpleTable];
+ If(flag_PRINTFIELDS==1)
+ Print[ us , OnElementsOf Omega, File "out/us_BT2.pos", Name "us_BT2"];
+ EndIf
+ }
+ }
+ { Name postop_modal_PML; NameOfPostProcessing postpro_modal_PML ;
+ Operation {
+ Print [ev_re , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesReal_PML.txt"];
+ Print [ev_im , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesImag_PML.txt"];
+ Print [vol[PrintPoint] , OnGlobal, TimeStep 0, StoreInVariable $vol,Format Table];
+ For j In {0:(neig-1)}
+ Print [eig_re[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_re"];
+ Print [eig_im[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_im"];
+ EndFor
+ Print [u , OnElementsOf Omega, File "out/u_PML.pos", Name "u_PML"];
+ If (flag_manualexp==1)
+ Print [vol[PrintPoint] , OnGlobal, TimeStep 0, StoreInVariable $vol,Format Table];
+ For i In {0:#tabrealomegasexpansion()-1}
+ For j In {0:(neig_PML-1)}
+ Print [overlap~{i}[Scat] , OnGlobal, TimeStep j, StoreInVariable $overlap~{i}~{j},Format Table];
+ EndFor
+ EndFor
+ For j In {0:(neig_PML-1)}
+ Print [den[Omega] , OnGlobal, TimeStep j, StoreInVariable $den~{j},Format Table];
+ EndFor
+ For j In {0:(neig_PML-1)}
+ Print [ev[PrintPoint] , OnGlobal, TimeStep j, StoreInVariable $ev~{j}, Format SimpleTable];
+ EndFor
+ For i In {0:(#tabrealomegasexpansion()-1)}
+ Print [us_PML_modal_manual~{i} , OnElementsOf Omega, StoreInField (stored_exp_field_manual_PML+i), File Sprintf["out/us_PML_modal_manual_%03g.pos",i], Name Sprintf["us_PML_modal_manual_%03g.pos",i], LastTimeStepOnly];
+ EndFor
+ For i In {0:(#tabrealomegasexpansion()-1)}
+ Print [intnorm2utot~{i}[Scat] , OnGlobal , Format SimpleTable,File Sprintf["out/ACC_PML_modal_manual_%03g.out",i],LastTimeStepOnly];
+ EndFor
+ EndIf
+ }
+ }
+ { Name postop_modal_AR_BT2; NameOfPostProcessing postpro_modal_AR_BT2 ;
+ Operation {
+ Print [ev_re , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesReal_BT2.txt"];
+ Print [ev_im , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesImag_BT2.txt"];
+ Print [vol[PrintPoint] , OnGlobal, TimeStep 0, StoreInVariable $vol,Format Table];
+ For j In {0:(neig-1)}
+ Print [eig_re[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_re"];
+ Print [eig_im[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_im"];
+ EndFor
+ Print [u , OnElementsOf Omega, File "out/u_BT2.pos", Name "u_BT2"];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Print [u , OnElementsOf Omega, StoreInField (stored_exp_field_AR_BT2+k), File Sprintf["out/us_BT2_modal_%03g.pos",k],Name "us_BT2_modal",TimeStep {neig_BT+k}];
+ Print [intnorm2utot~{k}[Scat] , OnGlobal, File Sprintf["out/ACC_BT2_modal_%03g.out",k], Format SimpleTable, TimeStep {neig_BT+k}];
+ EndFor
+ }
+ }
+ { Name postop_modal_AR_PML; NameOfPostProcessing postpro_modal_AR_PML ;
+ Operation {
+ Print [ev_re , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesReal_PML.txt"];
+ Print [ev_im , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesImag_PML.txt"];
+ // Print [u , OnElementsOf Omega, File "out/u_PML.pos", Name "u_BT2"];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Print [u , OnElementsOf Omega, StoreInField (stored_exp_field_AR_PML+k), File Sprintf["out/us_PML_modal_%03g.pos",k], Name "us_pml_modal",TimeStep {neig_PML+k}];
+ Print [intnorm2utot~{k}[Scat], OnGlobal , Format SimpleTable,File Sprintf["out/ACC_PML_modal_%03g.out",k],LastTimeStepOnly];
+ EndFor
+ }
+ }
+ { Name postop_modal_AR_mur; NameOfPostProcessing postpro_modal_AR_mur ;
+ Operation {
+ Print [ev_re , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesReal_mur.txt"];
+ Print [ev_im , OnElementsOf PrintPoint , Format TimeTable, File "out/EigenValuesImag_mur.txt"];
+ Print [vol[PrintPoint] , OnGlobal, TimeStep 0, StoreInVariable $vol,Format Table];
+ For j In {0:(neig-1)}
+ Print [eig_re[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_re"];
+ Print [eig_im[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_im"];
+ EndFor
+ // Print [u , OnElementsOf Omega, File "out/u_mur.pos", Name "u_BT2"];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Print [u , OnElementsOf Omega,StoreInField (stored_exp_field_AR_mur+k), File Sprintf["out/us_mur_modal_%03g.pos",k], Name "us_mur_modal",TimeStep {neig_BT+k}];
+ Print [intnorm2utot~{k}[Scat] , OnGlobal, File Sprintf["out/ACC_mur_modal_%03g.out",k], Format SimpleTable, TimeStep {neig_BT+k}];
+ EndFor
+ }
+ }
+ For k In {0:#tabrealomegasexpansion()-1}
+ { Name postop_gmsh~{k} ; NameOfPostProcessing postpro_modal_AR_BT2;
+ Operation {
+ Echo[Str["Geometry.Points = 0;","Geometry.Color.Curves = {0,0,0};","Mesh.SurfaceEdges = 0;",
+ "k=PostProcessing.NbViews-1;","View[k].Visible = 0;","View.Visible = 0;","View[0].Visible = 1;",
+ Sprintf["View[%g].Visible = 0;",k],
+ Sprintf["View[%g].Visible = 0;",nb_Omegas+k],
+ Sprintf["View[%g].Visible = 1;",k+1],
+ Sprintf["View[%g].ColormapNumber = 7;",k+1],
+ Sprintf["View[%g].IntervalsType = 3;",k+1],
+ Sprintf["View[%g].NbIso = 30;",k+1],
+ Sprintf["View[%g].Visible = 1;",nb_Omegas+k+1],
+ Sprintf["View[%g].ColormapNumber = 7;",nb_Omegas+k+1],
+ Sprintf["View[%g].IntervalsType = 3;",nb_Omegas+k+1],
+ Sprintf["View[%g].NbIso = 30;",nb_Omegas+k+1]], File "tmp0.geo"];
+ // Echo[Str["Plugin(Annotate).Text= \"TOP : DIRECT -- BOTTOM : EXPANSION\";" , "Plugin(Annotate).Run;"]];
+ }
+ }
+ EndFor
+}
+
+slepc_options_rg = Sprintf[" -rg_interval_endpoints %.8lf,%.8lf,%.8lf,%.8lf ",eig_min_re,eig_max_re,eig_min_im,eig_max_im];
+DefineConstant[
+ R_ = {"res_nleig_rat_exp", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {StrCat["-pre ",str_res," -cal -slepc ",str_slepc_solver," ",slepc_options_rg] , Name "GetDP/9ComputeCommand", Visible 1}
+];
+
+DefineConstant[
+ eig_re_ = {0, Name "GetDP/eig_re", ReadOnly 1, Graph "10", Visible 1},
+ eig_im_ = {0, Name "GetDP/eig_im", ReadOnly 1, Graph "01", Visible 1}
+];
diff --git a/QuasiNormalModeExpansion/QNM_expansion_DtN_data.geo b/QuasiNormalModeExpansion/QNM_expansion_DtN_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..3ceb72459cc1fa7112aede78d6401edf1a8df8b3
--- /dev/null
+++ b/QuasiNormalModeExpansion/QNM_expansion_DtN_data.geo
@@ -0,0 +1,85 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+pp = "Model Options/";
+DefineConstant[
+ flag_DtN = {2 , Name StrCat[pp, "0DtN"], Choices {0="PML",1="ABC",2="BT2"} , ServerAction "Reset"},
+ flag_vector_case = {0 , Name StrCat[pp, "1polarization case"], Choices {0="(0,0,Ez)",1="(Ex,Ey,0)"} },
+ flag_o2 = {1 , Name StrCat[pp, "2FE order"], Choices {0="order 1",1="order 2"} },
+ neig = {120 , Name StrCat[pp, "3number of eigenvalues"] }
+];
+
+nm = 0.1;
+paramaille = 10;
+
+neig_PML = neig;
+neig_BT = neig;
+
+a_lat = 482.5*nm;
+lambda_m = a_lat;
+
+lambda0 = 500*nm;
+dx_scat = 200*nm;
+dy_scat = 400*nm;
+epsr_scat_re = 9;
+epsr_scat_im = 0.5;
+epsr_bg = 1;
+
+r_domain = 0.5*Sqrt[dx_scat^2+dy_scat^2]+100*nm;
+r_pml = r_domain+lambda0;
+
+eps0 = 8.854187817e-3*nm;
+mu0 = 400.*Pi*nm;
+cel = 1.0/Sqrt[eps0 * mu0];
+
+norm = dy_scat/(Pi*cel);
+
+omega_target_re = 0.7/norm;
+omega_target_im = 1/norm;
+
+target_re_PML = omega_target_re^2 - omega_target_im^2;
+target_im_PML = 2*omega_target_re*omega_target_im;
+target_re_BT = omega_target_im;
+target_im_BT = omega_target_re;
+
+eig_min_re = -5/norm;
+eig_max_re = 5/norm;
+eig_min_im = -8/norm;
+eig_max_im = 8/norm;
+
+If (flag_DtN==0)
+ flag_PML = 1;
+ flag_ABC_BAYLISS = 0;
+ slepc_options_rg = Sprintf[" -rg_interval_endpoints %.8lf,%.8lf,%.8lf,%.8lf ",eig_min_re,eig_max_re,eig_min_im,eig_max_im];
+ str_res = "res_modal_PML";
+ str_slepc_solver = " ";
+ domain_tot = 2*r_pml;
+Else
+ flag_PML = 0;
+ flag_ABC_BAYLISS = 1;
+ slepc_options_rg = " ";
+ If (flag_DtN==1)
+ str_res = "res_modal_AR_mur";
+ Else
+ str_res = "res_modal_AR_BT2";
+ EndIf
+ str_slepc_solver = " -nep_type nleigs -nep_rational -nep_target_magnitude -nep_view_pre -nep_monitor_all :temp_eigenvalues.txt ";
+ domain_tot = 2*r_domain;
+EndIf
+
+flag_PRINTVECTOR = 1;
+flag_PRINTFIELDS = 0;
+
+nb_Omegas = 60;
+RealOmega_min = 0.5 / norm;
+RealOmega_max = 3 / norm;
+theta0 = 30 * Pi/180;
+
+tabrealomegasexpansion = LinSpace[RealOmega_min,RealOmega_max,nb_Omegas];
+// Freq = 0;
+
+// FIXME
+flag_manualexp = 1;
+
+yshift = 1.1;
diff --git a/QuasiNormalModeExpansion/QNM_expansion_dispersive_media.geo b/QuasiNormalModeExpansion/QNM_expansion_dispersive_media.geo
new file mode 100644
index 0000000000000000000000000000000000000000..ab4fa7271b6bb78b98c9a60a32b8ed90b403ed54
--- /dev/null
+++ b/QuasiNormalModeExpansion/QNM_expansion_dispersive_media.geo
@@ -0,0 +1,98 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+Include "QNM_expansion_dispersive_media_data.geo";
+lc_cell = lambda_m/paramaille;
+lc_sq = lc_cell/2.;
+lc_src = lc_cell;
+
+Point(1) = {-a_lat , -a_lat/2, 0. , lc_cell};
+Point(2) = {-a_lat , a_lat/2, 0. , lc_cell};
+Point(3) = { 0 , -a_lat/2, 0. , lc_sq};
+Point(4) = { 0 , a_lat/2, 0. , lc_sq};
+Point(5) = { a_lat , -a_lat/2, 0. , lc_cell};
+Point(6) = { a_lat , a_lat/2, 0. , lc_cell};
+Point(12) = {-a_lat+0.7 *2*a_lat , -a_lat/2, 0. , lc_sq};
+Point(13) = {-a_lat+0.48*2*a_lat , a_lat/2, 0. , lc_sq};
+
+Point(7) = { cx_src , cy_src , 0. , lc_src};
+Point(8) = { cx_src+r_src , cy_src , 0. , lc_src};
+Point(9) = { cx_src , cy_src+r_src , 0. , lc_src};
+Point(10) = { cx_src-r_src , cy_src , 0. , lc_src};
+Point(11) = { cx_src , cy_src-r_src , 0. , lc_src};
+
+Circle(8) = {10, 7, 9};
+Circle(9) = {9, 7, 8};
+Circle(10) = {8, 7, 11};
+Circle(11) = {11, 7, 10};
+
+Line(12) = {7, 8};
+Line(13) = {7, 9};
+Line(14) = {7, 10};
+Line(15) = {7, 11};
+
+Line Loop(1) = {8, -13, 14};
+Plane Surface(1) = {-1};
+Line Loop(2) = {13, 9, -12};
+Plane Surface(2) = {-2};
+Line Loop(3) = {10, -15, 12};
+Plane Surface(3) = {-3};
+Line Loop(4) = {14, -11, -15};
+Plane Surface(4) = {4};
+
+Line(16) = {1, 3};
+Line(17) = {3, 12};
+Line(18) = {12, 5};
+Line(19) = {5, 6};
+Line(20) = {6, 4};
+Line(21) = {4, 13};
+Line(22) = {13, 2};
+Line(23) = {2, 1};
+Line(24) = {3, 13};
+Line(25) = {12, 4};
+Curve Loop(5) = {23, 16, 24, 22};
+Plane Surface(5) = {5};
+Curve Loop(6) = {24, -21, -25, -17};
+Plane Surface(6) = {-6};
+Curve Loop(7) = {25, -20, -19, -18};
+Curve Loop(8) = {8, 9, 10, 11};
+Plane Surface(7) = {-7, 8};
+
+Physical Surface(1) = {6}; // Dispersive medium
+Physical Surface(2) = {5,7}; // Background
+Physical Surface(3) = {1,2,3,4}; // disk green
+Physical Line(10) = {16,17,18,19,20,21,22,23}; // Outer box
+Physical Line(20) = {24,25}; // bound lag
+Physical Line(30) = {12}; // Edge Source Green X
+Physical Line(31) = {13}; // Edge Source Green Y
+Physical Point(100) = {7};
+
+// dummy lines to see expansion and direct solutions side by side
+pp=newp;
+Point(pp+1) = {-a_lat+ xshift*2*a_lat , -a_lat/2 , 0. , lc_cell};
+Point(pp+2) = {-a_lat+ xshift*2*a_lat , a_lat/2 , 0. , lc_cell};
+Point(pp+3) = { 0+ xshift*2*a_lat , -a_lat/2 , 0. , lc_sq};
+Point(pp+4) = { 0+ xshift*2*a_lat , a_lat/2 , 0. , lc_sq};
+Point(pp+5) = { a_lat+ xshift*2*a_lat , -a_lat/2 , 0. , lc_cell};
+Point(pp+6) = { a_lat+ xshift*2*a_lat , a_lat/2 , 0. , lc_cell};
+Point(pp+12) = {-a_lat+0.7 *2*a_lat+ xshift*2*a_lat , -a_lat/2, 0. , lc_sq};
+Point(pp+13) = {-a_lat+0.48*2*a_lat+ xshift*2*a_lat , a_lat/2, 0. , lc_sq};
+
+c=newl;
+Line(c+16) = {pp+1, pp+3};
+Line(c+17) = {pp+3, pp+12};
+Line(c+18) = {pp+12, pp+5};
+Line(c+19) = {pp+5, pp+6};
+Line(c+20) = {pp+6, pp+4};
+Line(c+21) = {pp+4, pp+13};
+Line(c+22) = {pp+13, pp+2};
+Line(c+23) = {pp+2, pp+1};
+Line(c+24) = {pp+3, pp+13};
+Line(c+25) = {pp+12, pp+4};
+
+// dummy point to enforce window center
+pp=newp;
+Point(pp+1) = {-a_lat+ xshift*2*a_lat , 1.5*a_lat , 0. , lc_cell};
+
+Geometry.Points = 0;
diff --git a/QuasiNormalModeExpansion/QNM_expansion_dispersive_media.pro b/QuasiNormalModeExpansion/QNM_expansion_dispersive_media.pro
new file mode 100644
index 0000000000000000000000000000000000000000..1af0bc3affbb27e2b405d2cd2af3de6490310845
--- /dev/null
+++ b/QuasiNormalModeExpansion/QNM_expansion_dispersive_media.pro
@@ -0,0 +1,307 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+Include "QNM_expansion_dispersive_media_data.geo";
+
+Group {
+ // Domains
+ Om_1 = Region[1];
+ Om_src = Region[3];
+ Om_2 = Region[{2,3}];
+ Om = Region[{Om_1,Om_2}];
+ Om_but_src = Region[{Om,-Om_src}];
+
+ Surfdiri = Region[10];
+ Bound = Region[20];
+ pt_src = Region[100];
+ green_x_src = Region[30];
+ green_y_src = Region[31];
+ PrintPoint = Region[100];
+ Om_plot = Region[{Om,-Om_src}];
+}
+
+Function{
+ I[] = Complex[0,1];
+ k0[] = $realomega/cel;
+ mur[] = TensorDiag[1,1,1];
+ epsr_nod[] = TensorDiag[1,1,1];
+ om_lorentz_1[] = Complex[om_lorentz_1_re,om_lorentz_1_im];
+ mur_direct[] = 1;
+ epsr_direct[Om_1] = 1 + a_lorentz_1/($realomega-om_lorentz_1[]) - a_lorentz_1/($realomega+Conj[om_lorentz_1[]]);
+ epsr_direct[Om_2] = 1;
+ epsr_mau1[] = 1 + I[]*a_lorentz_1/($Eigenvalue-I[]*om_lorentz_1[]) - I[]*a_lorentz_1/($Eigenvalue+I[]*Conj[om_lorentz_1[]]);
+
+ If (flag_vector_case==0)
+ source_RHS[Om_src] = Vector[0,0,1];
+ source_RHS[Om_but_src] = Vector[0,0,0];
+ Else
+ source_RHS[Om_src] = Vector[1,0,0];
+ source_RHS[Om_but_src] = Vector[0,0,0];
+ EndIf
+}
+
+Constraint {
+ {Name Dirichlet; Type Assign;
+ Case {
+ { Region Surfdiri ; Value 0.; }
+ }
+ }
+
+ { Name Dirichlet_src;
+ Case {
+ { Region green_x_src;
+ Type Assign; Value 0. ;
+ }
+ }
+ }
+}
+
+Jacobian {
+ {Name JVol ; Case { { Region All ; Jacobian Vol ; } } }
+ {Name JSur ; Case { { Region All ; Jacobian Sur ; } } }
+ {Name JLin ; Case { { Region All ; Jacobian Lin ; } } }
+}
+
+Integration {
+ { Name Int_1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 12 ; }
+ }
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { Name Hgrad_perp; Type Form1P;
+ BasisFunction {
+ { Name un; NameOfCoef un; Function BF_PerpendicularEdge_1N; Support Region[Om]; Entity NodesOf[All]; }
+ { Name un2; NameOfCoef un2; Function BF_PerpendicularEdge_2E; Support Region[Om]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name un3; NameOfCoef un3; Function BF_PerpendicularEdge_2F; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name un4; NameOfCoef un4; Function BF_PerpendicularEdge_3E; Support Region[Om]; Entity EdgesOf[All]; }
+ { Name un5; NameOfCoef un5; Function BF_PerpendicularEdge_3F; Support Region[Om]; Entity FacetsOf[All]; }
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un; EntityType NodesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ If(flag_o2==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ EndIf
+ }
+ }
+
+ { Name Hcurl; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge ; Support Region[Om]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[Om]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om]; Entity EdgesOf[All]; }
+ EndIf
+ }
+ Constraint {
+ { NameOfCoef un ; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un3; EntityType FacetsOf; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType FacetsOf; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+
+ }
+ }
+
+ { Name Hcurl_direct_green_lin; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un ; Function BF_Edge ; Support Region[{Om,green_x_src}]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E; Support Region[{Om,green_x_src}]; Entity EdgesOf[All]; }
+ If(flag_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_a; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_b; Support Region[Om]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E ; Support Region[Om]; Entity EdgesOf[All]; }
+ EndIf
+ }
+ Constraint{
+ { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet_src; }
+ { NameOfCoef un3; EntityType FacetsOf; NameOfConstraint Dirichlet_src; }
+ { NameOfCoef un4; EntityType FacetsOf; NameOfConstraint Dirichlet_src; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet_src; }
+ }
+ }
+}
+
+Formulation {
+ {Name modal_helmholtz; Type FemEquation;
+ Quantity{
+ If (flag_vector_case==0)
+ { Name u; Type Local; NameOfSpace Hgrad_perp;}
+ Else
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ EndIf
+ }
+ Equation {
+ Galerkin { Eig[ 1/mur[] * Dof{Curl u}, {Curl u} ]; Rational 1; In Om ; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ epsr_nod[]/cel^2 * Dof{u}, {u} ]; Rational 2; In Om_1; Jacobian JSur; Integration Int_1;}
+ Galerkin { Eig[ epsr_nod[]/cel^2 * Dof{u}, {u} ]; Rational 3; In Om_2; Jacobian JSur; Integration Int_1;}
+ }
+ }
+
+ // if other matrices, fails to print eigs (conflict between M1 matrices)
+ {Name form_RHS; Type FemEquation;
+ Quantity{
+ If (flag_vector_case==0)
+ { Name u; Type Local; NameOfSpace Hgrad_perp;}
+ Else
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ EndIf
+ }
+ Equation{
+ Galerkin { [ 1. * Dof{u} , {u} ]; In Om; Jacobian JVol; Integration Int_1;}
+ Galerkin { [ source_RHS[] , {u} ]; In Om; Jacobian JVol; Integration Int_1; }
+ }
+ }
+
+ {Name direct_helmholtz; Type FemEquation;
+ Quantity{
+ If (flag_vector_case==0)
+ { Name u; Type Local; NameOfSpace Hgrad_perp;}
+ Else
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ EndIf
+ }
+ Equation {
+ Galerkin { [ -1/mur_direct[]*Dof{Curl u} , {Curl u}]; In Om ; Jacobian JSur; Integration Int_1; }
+ Galerkin { [ k0[]^2*epsr_direct[]*Dof{u} , {u}] ; In Om ; Jacobian JSur; Integration Int_1; }
+ Galerkin { [ -source_RHS[] , {u}] ; In Om_src; Jacobian JVol; Integration Int_1; }
+ }
+ }
+}
+
+Resolution {
+ { Name res_nleig_rat_exp;
+ System{
+ { Name V ; NameOfFormulation form_RHS ; Type ComplexValue; }
+ { Name M1; NameOfFormulation modal_helmholtz ; Type ComplexValue; }
+ { Name M2; NameOfFormulation direct_helmholtz ; Type ComplexValue; }
+ }
+ Operation{
+ CreateDir["out"];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ GenerateListOfRHS[V,Om,#tabrealomegasexpansion()];
+ EndFor
+ GenerateSeparate[M1];
+ EigenSolveAndExpand[M1,neig,eig_target_re,eig_target_im,
+ { {1}, {-1,-2*om_lorentz_1_im,2*a_lorentz_1*om_lorentz_1_re-om_lorentz_1_mod^2,0,0}, {1,0,0} } ,
+ { {1}, {-1,-2*om_lorentz_1_im,-om_lorentz_1_mod^2}, {1} } ,
+ tabrealomegasexpansion(),
+ V];
+ PostOperation[postop_lorentz_nleig_rat_exp];
+ For k In {0:#tabrealomegasexpansion()-1}
+ Evaluate[$realomega=tabrealomegasexpansion(k)];
+ Generate[M2];
+ Solve[M2];
+ PostOperation[postop_direct_helmholtz];
+ EndFor
+ For k In {0:#tabrealomegasexpansion()-1}
+ PostOperation[postop_gmsh~{k}];
+ Sleep[0.3];
+ EndFor
+ }
+ }
+}
+
+PostProcessing {
+ { Name postpro_lorentz_nleig_rat_exp; NameOfFormulation modal_helmholtz;
+ Quantity {
+ { Name intnormu ; Value { Integral { [ Norm[{u}] ] ; In Om ; Jacobian JSur; Integration Int_1 ; } } }
+ { Name intnormu_rod ; Value { Integral { [ Norm[{u}] ] ; In Om_1 ; Jacobian JSur; Integration Int_1 ;} } }
+ { Name u ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+ { Name u_im; Value { Local { [ Im[{u}] ] ; In Om; Jacobian JSur; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JSur; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JSur; } } }
+ { Name vol ; Value { Integral { [1] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_re ; Value { Integral { [$EigenvalueReal/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ { Name eig_im ; Value { Integral { [$EigenvalueImag/$vol] ; In PrintPoint ; Jacobian JVol; Integration Int_1 ; } } }
+ }
+ }
+
+ { Name postpro_direct_helmholtz; NameOfFormulation direct_helmholtz;
+ Quantity {
+ { Name eps ; Value { Local { [ epsr_direct[]] ; In Om ; Jacobian JSur; } } }
+ { Name source_RHS ; Value { Local { [ source_RHS[] ] ; In Om; Jacobian JSur; } } }
+ { Name u_direct ; Value { Local { [ {u} ] ; In Om; Jacobian JSur; } } }
+ { Name u_direct_im ; Value { Local { [ Im[{u}] ] ; In Om; Jacobian JSur; } } }
+ { Name intnormu_rod ; Value { Integral { [ Norm[{u}] ] ; In Om_1 ; Jacobian JSur; Integration Int_1 ;} } }
+ }
+ }
+}
+
+
+PostOperation {
+ { Name postop_lorentz_nleig_rat_exp; NameOfPostProcessing postpro_lorentz_nleig_rat_exp ;
+ Operation {
+ Print [vol[PrintPoint] , OnGlobal, TimeStep 0, StoreInVariable $vol,Format Table];
+ For j In {0:(neig-1)}
+ Print [eig_re[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_re"];
+ Print [eig_im[PrintPoint] , OnGlobal, TimeStep j , Format Table , SendToServer "GetDP/eig_im"];
+ EndFor
+ Print [EigenValuesReal, OnElementsOf PrintPoint, Format TimeTable, File "out/EigenValuesReal.txt"];
+ Print [EigenValuesImag, OnElementsOf PrintPoint, Format TimeTable, File "out/EigenValuesImag.txt"];
+ Print [ u , OnElementsOf Om, File Sprintf["out/u.pos"], EigenvalueLegend];
+ Print [ u , OnPoint {cx_src,-cy_src,0}, Format TimeTable, File "out/u_expand_on_probe_point.txt"];
+ For k In {0:#tabrealomegasexpansion()-1}
+ // Print [ u_im, OnElementsOf Om, File Sprintf["out/u_exp_%03g.out",k], Name Sprintf["u_expand_imag_%03g (left)",k] , EigenvalueLegend, TimeStep {neig+k}];
+ Print [ u_im, OnGrid {$A,$B, 0} {-a_lat:a_lat:a_lat/30, -a_lat/2:a_lat/2:a_lat/15, 0}, File Sprintf["out/u_exp_%03g.out",k], Name Sprintf["u_expand_imag_%03g (left)",k] , EigenvalueLegend, TimeStep {neig+k}];
+ EndFor
+ }
+ }
+ { Name postop_direct_helmholtz; NameOfPostProcessing postpro_direct_helmholtz ;
+ Operation {
+ Print[ u_direct_im, OnGrid {$A,$B, 0} {-a_lat:a_lat:a_lat/30, -a_lat/2:a_lat/2:a_lat/15, 0}, File "out/u_direct.pos" , ChangeOfCoordinates {$X+xshift*2*a_lat,$Y,$Z} , Name "u_direct_imag (right)"];
+ Print[ u_direct , OnPoint {cx_src,-cy_src,0}, Format TimeTable, File >> "out/u_direct_on_probe_point.txt"];
+ Print [intnormu_rod[Om_1], OnGlobal, Format TimeTable, File "out/intnormu_direct_on_src_point.out" ];
+ }
+ }
+ For k In {0:#tabrealomegasexpansion()-1}
+ { Name postop_gmsh~{k}; NameOfPostProcessing postpro_direct_helmholtz ;
+ Operation {
+ Echo[Str["Geometry.Points = 0;","Geometry.Color.Curves = {0,0,0};","Mesh.SurfaceEdges = 0;",
+ "k=PostProcessing.NbViews-1;","View[k].Visible = 0;",
+ "View.Visible = 0;","View[0].Visible = 1;", "View[0].ColormapNumber = 7;",
+ Sprintf["View[%g].Visible = 0;",k],
+ Sprintf["View[%g].Visible = 0;",nb_Omegas+k],
+ Sprintf["View[%g].ColormapNumber = 7;",k+1],
+ Sprintf["View[%g].ColormapNumber = 7;",nb_Omegas+k+1],
+ Sprintf["View[%g].Visible = 1;",k+1],
+ Sprintf["View[%g].Visible = 1;",nb_Omegas+k+1]], File "tmp0.geo"];
+ // Echo[Str["Plugin(Annotate).Text= "TOP : DIRECT -- BOTTOM : EXPANSION\";" , "Plugin(Annotate).Run;"]];
+ }
+ }
+ EndFor
+}
+
+
+slepc_options_rg = Sprintf[" -rg_interval_endpoints %.8lf,%.8lf,%.8lf,%.8lf ",eig_min_re,eig_max_re,eig_min_im,eig_max_im];
+DefineConstant[
+ R_ = {"res_nleig_rat_exp", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {StrCat["-pre res_nleig_rat_exp -cal -slepc -nep_type nleigs -nep_rational -nep_target_magnitude -nep_view_pre -nep_ncv 300 -nep_monitor_all :temp_eigenvalues.txt",slepc_options_rg] , Name "GetDP/9ComputeCommand", Visible 1}
+];
+
+DefineConstant[
+ eig_re_ = {0, Name "GetDP/eig_re", ReadOnly 1, Graph "10", Visible 0},
+ eig_im_ = {0, Name "GetDP/eig_im", ReadOnly 1, Graph "01", Visible 0}
+];
+
+// DefineConstant[
+// R_ = {"res_nleig_rat_exp", Name "GetDP/1ResolutionChoices", Visible 1},
+// C_ = {"-pre res_nleig_rat_exp -cal -slepc -nep_max_it 1000 -nep_type nleigs -nep_rational -nep_target_magnitude -nep_ncv 100 -rg_interval_endpoints 0.00000000,0.05763407,0.00390394,3.12315286 -nep_view_pre -nep_monitor_all :temp_eigenvalues.txt", Name "GetDP/9ComputeCommand", Visible 1}
+// ];
diff --git a/QuasiNormalModeExpansion/QNM_expansion_dispersive_media_data.geo b/QuasiNormalModeExpansion/QNM_expansion_dispersive_media_data.geo
new file mode 100644
index 0000000000000000000000000000000000000000..25e6ceb285fadd2fafcd9ca40897940b2ea8932d
--- /dev/null
+++ b/QuasiNormalModeExpansion/QNM_expansion_dispersive_media_data.geo
@@ -0,0 +1,49 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+///////////////////////////////
+
+pp = "Model Options/";
+DefineConstant[
+ flag_vector_case = {1 , Name StrCat[pp, "0polarization case"], Choices {0="(0,0,Ez)",1="(Ex,Ey,0)"} },
+ flag_o2 = {1 , Name StrCat[pp, "1FE order"], Choices {0="order 1",1="order 2"} },
+ neig = {100 , Name StrCat[pp, "0number of eigenvalues"]}
+];
+
+nm = 0.1;
+N_true_poles = 1;
+paramaille = 5;
+
+a_lat = 482.5*nm;
+lambda_m = a_lat;
+
+eps0 = 8.854187817e-3*nm;
+mu0 = 400.*Pi*nm;
+cel = 1.0/Sqrt[eps0 * mu0];
+norm = a_lat/(2.*Pi*cel);
+
+om_p_josab = .4*norm;
+om_0_josab = 0.3*norm;
+gama_josab = 0.05*norm;
+om_lorentz_1_im = -gama_josab/2;
+om_lorentz_1_re = Sqrt[om_0_josab^2-om_lorentz_1_im^2];
+a_lorentz_1 = -om_p_josab^2/(2*om_lorentz_1_re);
+om_lorentz_1_mod = Sqrt[om_lorentz_1_re^2+om_lorentz_1_im^2];
+
+eig_target_re = om_lorentz_1_im/4;
+eig_target_im = 1/norm;
+eig_max_im = 8/norm;
+eig_min_im = 0.01/norm;
+eig_min_re = 0;
+eig_max_re = 0.8*Fabs[om_lorentz_1_im];
+
+cx_src = 0.75*2*a_lat - a_lat;
+cy_src = 0.75*a_lat - a_lat/2;
+r_src = a_lat/100;
+
+nb_Omegas = 50;
+RealOmega_min = 0.1;
+RealOmega_max = 1.1;
+tabrealomegasexpansion = LinSpace[RealOmega_min,RealOmega_max,nb_Omegas];
+Freq = 0;
+
+xshift = 1.1;