Select Git revision
onelab.html
-
Christophe Geuzaine authoredChristophe Geuzaine authored
scattering.pro 33.78 KiB
///////////////////////////////
// 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
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~{pe}; 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_M~{pe}[] , {u} ]; In All_domains; Jacobian JVol; Integration Int_1; }
}
}
{Name VPWN_helmholtz_vector~{pe}; 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
{Name VPWMN_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 { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*
($isN ? Nnm[1,$NE,$ME,XYZ[],k_Out] : Mnm[1,$NE,$ME,XYZ[],k_Out])
, {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
}
}
{Name VPWN_helmholtz_vector_test; 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 { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm[1,$NE,$ME,XYZ[],k_Out], {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
}
}
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 {
{ Name T; NameOfFormulation VPWMN_helmholtz_vector; Type ComplexValue; }
}
Operation {
CreateDir[Str[myDir]];
Evaluate[Python[]{"scattering_init.py"}];
For pe In {1:p_max}
Evaluate[ $isN = 0 ];
Evaluate[ $PE = pe ];
Evaluate[ $NE = Floor[Sqrt[$PE]] ];
Evaluate[ $ME = $NE*($NE+1) - Floor[$PE] ];
If (pe==1)
Generate[T];
Solve[T];
EndIf
GenerateRHS[T];
SolveAgain[T];
PostOperation[VPWM_postop~{pe}];
Evaluate[$isN=1];
GenerateRHS[T];
SolveAgain[T];
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 VPWMN_helmholtz_vector; NameOfSystem T;
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~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
{ Name norma~{pe}~{po} ; Value { Integral { [ Mnm_out~{pe}[]*Conj[Mnm_out~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
EndFor
}
}
{ Name VPWN_postpro~{pe}; NameOfFormulation VPWMN_helmholtz_vector; NameOfSystem T;
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]]]];
Print[norma~{pe}~{po}[SurfInt], OnGlobal , Format Table, File StrCat[myDir,StrCat[StrCat["norma_",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