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

solveagain version

parent 791e9998
No related branches found
No related tags found
No related merge requests found
Pipeline #4452 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
ElectromagneticScattering/scattering_solveagain.pro 0 → 100644
+ 704
− 0
View file @ 33747375
///////////////////////////////
// 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment