diff --git a/ElectromagneticScattering/scattering_solveagain.pro b/ElectromagneticScattering/scattering_solveagain.pro
new file mode 100644
index 0000000000000000000000000000000000000000..dbe5e4d5ebbdc833a7a9d649183db686a444b77f
--- /dev/null
+++ b/ElectromagneticScattering/scattering_solveagain.pro
@@ -0,0 +1,704 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattering.pro //
+///////////////////////////////
+
+Include "scattering_data.geo";
+myDir = "run_results/";
+
+Group {
+ If(flag_cartpml==1)
+ PMLxyz = Region[1000];
+ PMLxz = Region[1001];
+ PMLyz = Region[1002];
+ PMLxy = Region[1003];
+ PMLz = Region[1004];
+ PMLy = Region[1005];
+ PMLx = Region[1006];
+ Scat_In = Region[1008];
+ Scat_Out = Region[1007];
+ // Domains
+ Domain = Region[{Scat_In,Scat_Out}];
+ PMLs = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+ All_domains = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
+ SurfInt = Region[{}];
+ EndIf
+ If(flag_cartpml==0)
+ Scat_In = Region[1];
+ Scat_Out = Region[2];
+ Domain = Region[{Scat_In,Scat_Out}];
+ PMLs = Region[3];
+ All_domains = Region[{Scat_In,Scat_Out,PMLs}];
+ SurfInt = Region[{10}];
+ EndIf
+ PrintPoint = Region[2000];
+}
+
+Function{
+ // test={1.4,2.6,3.7};
+ // For i In {0:#test()-1}
+ // Printf("%g",test(i));
+ // EndFor
+ I[] = Complex[0,1];
+ avoid_sing = 1.e-14;
+ mu0 = 4*Pi*100.0*nm;
+ epsilon0 = 8.854187817e-3*nm;
+ cel = 1.0/Sqrt[epsilon0 * mu0];
+ Freq = cel/lambda0;
+ omega0 = 2.0*Pi*Freq;
+ k_Out = 2.0*Pi*Sqrt[epsr_Out_re]/lambda0;
+ r3D_sph[] = Sqrt[X[]^2+Y[]^2+Z[]^2];
+ r2D_sph[] = Sqrt[X[]^2+Y[]^2];
+ cos_theta[] = Z[]/(r3D_sph[]+avoid_sing);
+ theta[] = Acos[cos_theta[]];
+ phi[] = Atan2[-Y[],-X[]]+Pi;
+
+ For kr In {0:nb_cuts-1}
+ r_sph~{kr}=r_sph_min+kr*(r_sph_max-r_sph_min)/(nb_cuts-1);
+ EndFor
+
+ If (flag_study==0)
+ Ae = 1.0;
+ Ah = Ae*Sqrt[epsilon0*epsr_In_re/mu0];
+ alpha0 = k_Out*Sin[theta0]*Cos[phi0];
+ beta0 = k_Out*Sin[theta0]*Sin[phi0];
+ gamma0 = k_Out*Cos[theta0];
+ Ex0 = Ae * Cos[psi0]*Cos[theta0]*Cos[phi0] - Ae* Sin[psi0]*Sin[phi0];
+ Ey0 = Ae * Cos[psi0]*Cos[theta0]*Sin[phi0] + Ae* Sin[psi0]*Cos[phi0];
+ Ez0 = -Ae * Cos[psi0]*Sin[theta0];
+ Hx0 = -1/(omega0*mu0)*(beta0 * Ez0 - gamma0 * Ey0);
+ Hy0 = -1/(omega0*mu0)*(gamma0 * Ex0 - alpha0 * Ez0);
+ Hz0 = -1/(omega0*mu0)*(alpha0 * Ey0 - beta0 * Ex0);
+ Prop[] = Ae * Complex[ Cos[alpha0*X[]+beta0*Y[]+gamma0*Z[]] , siwt*Sin[alpha0*X[]+beta0*Y[]+gamma0*Z[]] ];
+ Einc[] = Vector[Ex0*Prop[],Ey0*Prop[],Ez0*Prop[]];
+ Hinc[] = Vector[Hx0*Prop[],Hy0*Prop[],Hz0*Prop[]];
+ Pinc = 0.5*Ae*Ae*Sqrt[epsilon0*epsr_In_re/mu0] * Cos[theta0];
+ EndIf
+ If (flag_study==2)
+ // Green Dyadic
+ normrr_p[] = Sqrt[(X[]-x_p)^2+(Y[]-y_p)^2+(Z[]-z_p)^2];
+ Gsca[] = Complex[Cos[-1.*siwt*k_Out*normrr_p[]],
+ Sin[-1.*siwt*k_Out*normrr_p[]]]
+ /(4.*Pi*normrr_p[]);
+ GDfact_diag[] = 1.+(I[]*k_Out*normrr_p[]-1.)/(k_Out*normrr_p[])^2;
+ GDfact_glob[] = (3.-3.*I[]*k_Out*normrr_p[]-(k_Out*normrr_p[])^2)/
+ (k_Out*normrr_p[]^2)^2;
+ Dyad_Green2[]=
+ Gsca[]*(
+ GDfact_diag[]*TensorDiag[1,1,1]
+ + GDfact_glob[]*
+ Tensor[(X[]-x_p)*(X[]-x_p),(X[]-x_p)*(Y[]-y_p),(X[]-x_p)*(Z[]-z_p),
+ (Y[]-y_p)*(X[]-x_p),(Y[]-y_p)*(Y[]-y_p),(Y[]-y_p)*(Z[]-z_p),
+ (Z[]-z_p)*(X[]-x_p),(Z[]-z_p)*(Y[]-y_p),(Z[]-z_p)*(Z[]-z_p)]);
+ // Getdp Func
+ Dyad_Green[] = DyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+ Einc_Green_0[] = Dyad_Green2[]*Vector[1,0,0];
+ Einc_Green_1[] = Dyad_Green2[]*Vector[0,1,0];
+ Einc_Green_2[] = Dyad_Green2[]*Vector[0,0,1];
+ CurlDyad_Green[] = CurlDyadGreenHom[x_p,y_p,z_p,XYZ[],k_Out,siwt];
+ Hinc_Green_0[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[1,0,0];
+ Hinc_Green_1[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,1,0];
+ Hinc_Green_2[] = siwt*I[]/(mu0*omega0)*CurlDyad_Green[]*Vector[0,0,1];
+ // test[] = DyadGreenHom[x_p,y_p,z_p,x_p+1e-1*nm,y_p+1e-1*nm,z_p+1e-1*nm,lambda0,Sqrt[epsr_Out_re]]*Vector[1,0,0];
+ EndIf
+
+ a_pml = 1.;
+ b_pml = -siwt;
+
+ eigtarget = (2.*Pi*cel/lambda0)^2;
+ eigfilter = 1e-4*Sqrt[eigtarget];
+ EigFilter[] = (Norm[$EigenvalueReal] > eigfilter);
+
+ If(flag_cartpml==1)
+ sx[Scat_In] = 1.;
+ sy[Scat_In] = 1.;
+ sz[Scat_In] = 1.;
+ sx[Scat_Out] = 1.;
+ sy[Scat_Out] = 1.;
+ sz[Scat_Out] = 1.;
+ sx[PMLxyz] = Complex[a_pml,b_pml];
+ sy[PMLxyz] = Complex[a_pml,b_pml];
+ sz[PMLxyz] = Complex[a_pml,b_pml];
+ sx[PMLxz] = Complex[a_pml,b_pml];
+ sy[PMLxz] = 1.0;
+ sz[PMLxz] = Complex[a_pml,b_pml];
+ sx[PMLyz] = 1.0;
+ sy[PMLyz] = Complex[a_pml,b_pml];
+ sz[PMLyz] = Complex[a_pml,b_pml];
+ sx[PMLxy] = Complex[a_pml,b_pml];
+ sy[PMLxy] = Complex[a_pml,b_pml];
+ sz[PMLxy] = 1.0;
+ sx[PMLx] = Complex[a_pml,b_pml];
+ sy[PMLx] = 1.0;
+ sz[PMLx] = 1.0;
+ sx[PMLy] = 1.0;
+ sy[PMLy] = Complex[a_pml,b_pml];
+ sz[PMLy] = 1.0;
+ sx[PMLz] = 1.0;
+ sy[PMLz] = 1.0;
+ sz[PMLz] = Complex[a_pml,b_pml];
+ Lxx[] = sy[]*sz[]/sx[];
+ Lyy[] = sz[]*sx[]/sy[];
+ Lzz[] = sx[]*sy[]/sz[];
+
+ epsilonr_In[] = Complex[epsr_In_re , epsr_In_im];
+ epsilonr_Out[] = Complex[epsr_Out_re , epsr_Out_im];
+
+ epsilonr[Scat_In] = epsilonr_In[] * TensorDiag[1.,1.,1.];
+ epsilonr[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr[PMLs] = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+ epsilonr1[Scat_In] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[PMLs] = epsr_Out_re * TensorDiag[Lxx[],Lyy[],Lzz[]];
+
+ mur[Scat_In] = TensorDiag[1.,1.,1.];
+ mur[Scat_Out] = TensorDiag[1.,1.,1.];
+ mur[PMLs] = TensorDiag[Lxx[],Lyy[],Lzz[]];
+ EndIf
+
+ If(flag_cartpml==0)
+ sr[] = Complex[a_pml,b_pml];
+ sphi[] = Complex[1,0];
+ stheta[] = Complex[1,0];
+ r_tild[] = r_pml_in + (r3D_sph[] - r_pml_in) * sr[];
+ theta_tild[] = theta[];
+ pml_tens_temp1[] = TensorDiag[(r_tild[]/r3D_sph[])^2 * sphi[]*stheta[]/sr[],
+ (r_tild[]/r3D_sph[]) * sr[]*stheta[]/sphi[],
+ (r_tild[]/r3D_sph[]) * sphi[]*sr[]/stheta[]];
+ pml_tens_temp2[] = Rotate[pml_tens_temp1[],0,-theta[]-Pi/2,0];
+ pml_tens[] = Rotate[pml_tens_temp2[],0,0,-phi[]];
+
+ epsilonr_In[] = Complex[epsr_In_re , epsr_In_im];
+ epsilonr_Out[] = Complex[epsr_Out_re , epsr_Out_im];
+
+ epsilonr[Scat_In] = epsilonr_In[] * TensorDiag[1.,1.,1.];
+ epsilonr[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr[PMLs] = pml_tens[];
+
+ epsilonr1[Scat_In] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[Scat_Out] = epsilonr_Out[] * TensorDiag[1.,1.,1.];
+ epsilonr1[PMLs] = pml_tens[];
+
+ mur[Scat_In] = TensorDiag[1.,1.,1.];
+ mur[Scat_Out] = TensorDiag[1.,1.,1.];
+ mur[PMLs] = pml_tens[];
+ EndIf
+
+ If(flag_study==RES_PW)
+ source[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc[];
+ EndIf
+ If(flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ ne = Floor[Sqrt[pe]];
+ me = ne*(ne+1) - Floor[pe];
+ Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
+ Nnm_source~{pe}[] = Nnm[1,ne,me,XYZ[],k_Out];
+ Mnm_out~{pe}[] = Mnm[4,ne,me,XYZ[],k_Out];
+ Nnm_out~{pe}[] = Nnm[4,ne,me,XYZ[],k_Out];
+ source_M~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
+ source_N~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
+ EndFor
+ EndIf
+
+ //SetVariable[1]{$pe};
+
+ f[] = ($pe = 1);
+
+ If (flag_study==RES_GREEN)
+ sourceG_0[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_0[];
+ sourceG_1[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_1[];
+ sourceG_2[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Einc_Green_2[];
+ EndIf
+}
+
+Constraint {
+ // {Name Dirichlet; Type Assign;
+ // Case {
+ // { Region SurfDirichlet; 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 Line2 ; NumberOfPoints 4 ; }
+ { GeoElement Triangle2 ; NumberOfPoints 6 ; }
+ { GeoElement Tetrahedron2 ; NumberOfPoints 15 ; }
+ // { GeoElement Line2 ; NumberOfPoints 4 ; }
+ // { GeoElement Triangle2 ; NumberOfPoints 6 ; }
+ // { GeoElement Tetrahedron2 ; NumberOfPoints 15 ; }
+ }
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { Name Hcurl; Type Form1;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Edge;
+ Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
+ { Name sn2; NameOfCoef un2; Function BF_Edge_2E;
+ Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
+ If (is_FEM_o2==1)
+ { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
+ Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
+ { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
+ Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
+ { Name sn5; NameOfCoef un5; Function BF_Edge_4E;
+ Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
+ EndIf
+ }
+ Constraint {
+ // { NameOfCoef un; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ // { NameOfCoef un2; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ }
+ }
+}
+
+Formulation {
+ If (flag_study==RES_QNM)
+ {Name QNM_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [ -1/mur[] * Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { DtDtDof[-1. * 1/cel^2 * epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndIf
+
+ If (flag_study==RES_PW)
+ {Name PW_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ // For pe In {1:p_max}
+ {Name VPWM_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ For pe In {1:p_max}
+ Galerkin { [($pe==pe ? source_M~{pe}[] : 0) , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ EndFor
+ }
+ }
+ {Name VPWN_helmholtz_vector; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [source_N~{pe}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ // EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ {Name GreenAll_helmholtz_vector~{ncomp}; Type FemEquation;
+ Quantity {
+ { Name u; Type Local; NameOfSpace Hcurl;}
+ }
+ Equation {
+ Galerkin { [-1/mur[]*Dof{Curl u} , {Curl u}]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [(omega0/cel)^2*epsilonr[]*Dof{u} , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ Galerkin { [sourceG~{ncomp}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
+ }
+ }
+ EndFor
+ EndIf
+}
+
+Resolution {
+ If (flag_study==RES_PW)
+ { Name res_PW_helmholtz_vector;
+ System {
+ { Name P; NameOfFormulation PW_helmholtz_vector; Type ComplexValue; Frequency Freq; }
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ Generate[P];
+ Solve[P];
+ PostOperation[PW_postop];
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ { Name res_VPWall_helmholtz_vector;
+ System {
+ // For pe In {1:p_max}
+ { Name M; NameOfFormulation VPWM_helmholtz_vector; Type ComplexValue; Frequency Freq; }
+ { Name N; NameOfFormulation VPWN_helmholtz_vector; Type ComplexValue; Frequency Freq; }
+ // EndFor
+ }
+ Operation {
+ Evaluate[ $c = 0.3 ];
+ While[ $c < 0.9 ]{
+ Evaluate[ $c = $c + 0.1 ];
+ }
+ Test[ $c > 1 ]{ // test performed at run-time
+ Print[{$c}, Format "$c = %g is greater than 1"]; // run-time message
+ }
+
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ Evaluate[$pe=1];
+
+ Generate[M];
+ Solve[M];
+ PostOperation[VPWM_postop];
+ Evaluate[$pe=2];
+ GenerateRHS[M];
+ SolveAgain[M];
+ // For pe In {1:p_max}
+ // GenerateRHS[M~{pe}];
+ // Solve[M~{pe}];
+ // PostOperation[VPWM_postop~{pe}];
+ // Generate[N~{pe}];
+ // Solve[N~{pe}];
+ // PostOperation[VPWN_postop~{pe}];
+ // EndFor
+
+ // For pe In {1:p_max}
+ // Generate[M~{pe}];
+ // Solve[M~{pe}];
+ // PostOperation[VPWM_postop~{pe}];
+ // Generate[N~{pe}];
+ // Solve[N~{pe}];
+ // PostOperation[VPWN_postop~{pe}];
+ // EndFor
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_GREEN)
+ { Name res_GreenAll_helmholtz_vector;
+ System {
+ For ncomp In {0:2}
+ { Name M~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp}; Type ComplexValue; Frequency Freq; }
+ EndFor
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ For ncomp In {0:2}
+ Generate[M~{ncomp}];
+ Solve[M~{ncomp}];
+ PostOperation[GreenAll_postop~{ncomp}];
+ EndFor
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+ If (flag_study==RES_QNM)
+ { Name res_QNM_helmholtz_vector;
+ System{
+ { Name E; NameOfFormulation QNM_helmholtz_vector; Type ComplexValue; }
+ }
+ Operation{
+ CreateDir[Str[myDir]];
+ Evaluate[Python[]{"scattering_init.py"}];
+ GenerateSeparate[E];
+ EigenSolve[E,neig,eigtarget,0,EigFilter[]];
+ PostOperation[QNM_postop];
+ Evaluate[Python[]{"scattering_post.py"}];
+ }
+ }
+ EndIf
+
+}
+
+PostProcessing {
+ If (flag_study==RES_QNM)
+ { Name QNM_postpro; NameOfFormulation QNM_helmholtz_vector; NameOfSystem E;
+ Quantity {
+ { Name u ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name EigenValuesReal; Value { Local{ [$EigenvalueReal]; In PrintPoint; Jacobian JVol; } } }
+ { Name EigenValuesImag; Value { Local{ [$EigenvalueImag]; In PrintPoint; Jacobian JVol; } } }
+ }
+ }
+ EndIf
+
+ If (flag_study==RES_PW)
+ { Name PW_postpro ; NameOfFormulation PW_helmholtz_vector;NameOfSystem P;
+ Quantity {
+ { Name E_tot ; Value { Local { [{u}+Einc[]]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}] ; In All_domains; Jacobian JVol; } } }
+ { Name H_tot ; Value { Local { [ siwt*I[]/(mur[]*mu0*omega0)*{Curl u} +Hinc[]]; In All_domains; Jacobian JVol; } } }
+ { Name normalized_losses ; Value { Integral { [ 0.5*omega0*epsilon0*Fabs[Im[epsilonr_In[]]]*(SquNorm[{u}+Einc[]]) ] ; In Scat_In ; Integration Int_1 ; Jacobian JVol ; } } }
+ { Name Poy_dif ; Value { Local { [ 0.5*Re[Cross[{u} , Conj[ siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]] ]; In All_domains; Jacobian JVol; } } }
+ { Name Poy_tot ; Value { Local { [ 0.5*Re[Cross[{u}+Einc[] , Conj[Hinc[]+siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]]] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector;NameOfSystem M;
+ Quantity {
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ { Name Mnm_source~{pe} ; Value { Local { [ Mnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+ For po In {1:p_max}
+ { Name a~{pe}~{po} ; Value { Integral { [ {u}*Conj[Mnm_out~{pe}[]]]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+ { Name norma~{pe}~{po} ; Value { Integral { [Mnm_out~{pe}[]*Conj[Mnm_out~{pe}[]]]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+ EndFor
+ }
+ }
+ { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector;NameOfSystem N;
+ Quantity {
+ { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph ; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}] - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ { Name Nnm_source~{pe} ; Value { Local { [ Nnm_source~{pe}[] ]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ { Name GreenAll_postpro~{ncomp}; NameOfFormulation GreenAll_helmholtz_vector~{ncomp};NameOfSystem M~{ncomp};
+ Quantity {
+ { Name E_tot~{ncomp} ; Value { Local {[{u}+Einc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+ { Name H_tot~{ncomp} ; Value { Local {[siwt*I[]/(mur[]*mu0*omega0)*{Curl u}+Hinc_Green~{ncomp}[]]; In All_domains; Jacobian JVol; } } }
+ { Name Einc_Green~{ncomp}; Value { Local {[ Einc_Green~{ncomp}[] ]; In All_domains; Jacobian JVol; } } }
+ { Name Hinc_Green~{ncomp}; Value { Local {[ Hinc_Green~{ncomp}[] ]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat~{ncomp} ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
+ { Name E_tot_im~{ncomp} ; Value { Local {[Im[{u}+Einc_Green~{ncomp}[]]]; In All_domains; Jacobian JVol; } } }
+ { Name Einc_Green_im~{ncomp}; Value { Local {[Im[ Einc_Green~{ncomp}[] ]]; In All_domains; Jacobian JVol; } } }
+ // { Name Einc_Green_im~{ncomp}; Value { Local {[test[]]; In All_domains; Jacobian JVol; } } }
+ { Name E_scat_im~{ncomp} ; Value { Local {[Im[{u} ]]; In All_domains; Jacobian JVol; } } }
+ { Name E_tot_sph~{ncomp}; Value { Local { [Vector[
+ (X[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + Y[]*CompY[{u}+Einc_Green~{ncomp}[]]
+ + Z[]*CompZ[{u}+Einc_Green~{ncomp}[]])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + Y[]*Z[]*CompY[{u}+Einc_Green~{ncomp}[]]
+ - (X[]^2+Y[]^2)*CompZ[{u}+Einc_Green~{ncomp}[]])
+ /(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}+Einc_Green~{ncomp}[]]
+ + X[]*CompY[{u}+Einc_Green~{ncomp}[]])
+ /r2D_sph[] ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name E_scat_sph~{ncomp}; Value { Local { [Vector[
+ (X[]*CompX[{u}] + Y[]*CompY[{u}] + Z[]*CompZ[{u}])/r3D_sph[],
+ (X[]*Z[]*CompX[{u}] + Y[]*Z[]*CompY[{u}]
+ - (X[]^2+Y[]^2)*CompZ[{u}])/(r3D_sph[]*r2D_sph[]),
+ (-Y[]*CompX[{u}] + X[]*CompY[{u}])/r2D_sph[]
+ ]];
+ In All_domains; Jacobian JVol; } } }
+ { Name H_scat~{ncomp} ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl u}]; In All_domains; Jacobian JVol; } } }
+ }
+ }
+ EndFor
+ EndIf
+}
+
+
+
+PostOperation {
+ If (flag_study==RES_QNM)
+ { Name QNM_postop; NameOfPostProcessing QNM_postpro ;
+ Operation {
+ Print [EigenValuesReal, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesReal.txt"]];
+ Print [EigenValuesImag, OnElementsOf PrintPoint , Format TimeTable, File StrCat[myDir,"EigenValuesImag.txt"]];
+ Print [ u , OnElementsOf All_domains, File StrCat[myDir,"eigenVectors.pos"], EigenvalueLegend];
+ }
+ }
+ EndIf
+ If (flag_study==RES_PW)
+ {Name PW_postop; NameOfPostProcessing PW_postpro ;
+ Operation {
+ Print [ E_tot , OnElementsOf Domain , File StrCat[myDir,"E_tot.pos"]];
+ Print [ E_scat , OnElementsOf All_domains, File StrCat[myDir,"E_scat.pos"]];
+ // Print [ E_scat_sph , OnElementsOf All_domains, File StrCat[myDir,"E_scat_sph.pos"]];
+ Print [ H_scat , OnElementsOf Domain , File StrCat[myDir,"H_scat.pos"]];
+ Print [ H_tot , OnElementsOf All_domains, File StrCat[myDir,"H_tot.pos"]];
+ Print [ Poy_dif , OnElementsOf All_domains, File StrCat[myDir,"Poy_dif.pos"]];
+ Print [ Poy_tot , OnElementsOf Domain , File StrCat[myDir,"Poy_tot.pos"]];
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_scat_onsphere_sph_PW",Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ }
+ }
+ EndIf
+ If (flag_study==RES_TMAT)
+ For pe In {1:p_max}
+ {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
+ Operation {
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_M",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ Print [ E_scat , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+ Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
+ For po In {1:p_max}
+ Print[a~{pe}~{po}[SurfInt], OnGlobal , Format Table, File StrCat[myDir,StrCat[StrCat["a_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
+ EndFor
+ }
+ }
+ {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
+ Operation {
+ For kr In {0:nb_cuts-1}
+ Print [ E_scat_sph , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_sph_N",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+ EndFor
+ Print [ E_scat , OnGrid
+ {r_sph~{kr}*Sin[$B]*Cos[$C], r_sph~{kr}*Sin[$B]*Sin[$C], r_sph~{kr}*Cos[$B]}
+ {r_sph~{kr}, {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat[StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]],Sprintf["_r%g.pos",kr]]],
+ Name StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]]];
+ }
+ }
+ EndFor
+ EndIf
+ If (flag_study==RES_GREEN)
+ For ncomp In {0:2}
+ {Name GreenAll_postop~{ncomp}; NameOfPostProcessing GreenAll_postpro~{ncomp} ;
+ Operation {
+ Print [ E_tot_sph~{ncomp} , OnGrid
+ {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_tot_onsphere_sph_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ Print [ E_tot~{ncomp} , OnGrid
+ {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_plot_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_plot_theta-1.0)},
+ {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_plot_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_plot_phi-1.0)} },
+ File StrCat[myDir,StrCat["E_tot_onsphere_cart_G",Sprintf["%g.pos",ncomp]]]];
+ Print [ E_scat_im~{ncomp} , OnPoint {x_p,y_p,z_p},
+ File StrCat[myDir,StrCat["E_scat_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ // Print [ Einc_Green_im~{ncomp} , OnPoint {x_p+0.1*nm,y_p+0.1*nm,z_p+0.1*nm},
+ // File StrCat[myDir,StrCat["E_inc_im_onpt_G",Sprintf["%g.dat",ncomp]]] , Format Table];
+ // If(flag_FF==1)
+ // Print [ E_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["E_tot_",Sprintf["%g.pos",ncomp]]]];
+ // Print [ H_tot~{ncomp}, OnElementsOf Scat_Out , File StrCat[myDir,StrCat["H_tot_",Sprintf["%g.pos",ncomp]]]];
+ // Print [ E_tot~{ncomp} , OnGrid
+ // {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ // {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ // {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ // File StrCat[myDir,StrCat["E_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+ // Print [ H_tot~{ncomp} , OnGrid
+ // {r_pml_in *Sin[$B]*Cos[$C], r_pml_in*Sin[$B]*Sin[$C], r_pml_in*Cos[$B]}
+ // {r_pml_in , {sph_scan : Pi-sph_scan+(Pi-2.0*sph_scan)/(10*(npts_theta-1.0)) : (Pi-2.0*sph_scan)/(npts_theta-1.0)},
+ // {sph_scan : 2.0*Pi-sph_scan+(2.0*Pi-2.0*sph_scan)/(10*(npts_phi-1.0)) : (2.0*Pi-2.0*sph_scan)/(npts_phi-1.0)} },
+ // File StrCat[myDir,StrCat["H_tot_onsphere_G",Sprintf["%g.pos",ncomp]]]];
+ //
+ // Echo[
+ // Str["Plugin(NearToFarField).Wavenumber = 2*Pi/lambda0;",
+ // "Plugin(NearToFarField).PhiStart = 0;",
+ // "Plugin(NearToFarField).PhiEnd = 2*Pi;",
+ // "Plugin(NearToFarField).NumPointsPhi = npts_phi;",
+ // "Plugin(NearToFarField).ThetaStart = 0;",
+ // "Plugin(NearToFarField).ThetaEnd = Pi;",
+ // "Plugin(NearToFarField).NumPointsTheta = npts_theta;",
+ // "Plugin(NearToFarField).EView = PostProcessing.NbViews-2;",
+ // "Plugin(NearToFarField).HView = PostProcessing.NbViews-1;",
+ // "Plugin(NearToFarField).Normalize = 0;",
+ // "Plugin(NearToFarField).dB = 0;",
+ // "Plugin(NearToFarField).NegativeTime = 0;",
+ // "Plugin(NearToFarField).RFar = r_pml_in;",
+ // "Plugin(NearToFarField).Run;"], File StrCat[myDir,"tmp1.geo"]] ;
+ // EndIf
+ }
+ }
+ EndFor
+ EndIf
+}
+
+If (flag_study==RES_QNM)
+ DefineConstant[
+ R_ = {"res_QNM_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+
+If (flag_study==RES_TMAT)
+ DefineConstant[
+ R_ = {"res_VPWall_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+If (flag_study==RES_PW)
+ DefineConstant[
+ R_ = {"res_PW_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
+If (flag_study==RES_GREEN)
+ DefineConstant[
+ R_ = {"res_GreenAll_helmholtz_vector", Name "GetDP/1ResolutionChoices", Visible 1},
+ C_ = {"-solve -pos -petsc_prealloc 200 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", Name "GetDP/9ComputeCommand", Visible 1},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}];
+EndIf
\ No newline at end of file