From c1f9635ad4deb43ad414bd3d772295a573953c2b Mon Sep 17 00:00:00 2001
From: Guillaume Demesy <guillaume.demesy@fresnel.fr>
Date: Fri, 7 May 2021 14:27:20 +0000
Subject: [PATCH] standalone/fast T-matrix computation (no python)
---
.../scattererTmatrix.pro | 267 ++++++++++++++++++
1 file changed, 267 insertions(+)
create mode 100644 ElectromagneticScattering/scattererTmatrix.pro
diff --git a/ElectromagneticScattering/scattererTmatrix.pro b/ElectromagneticScattering/scattererTmatrix.pro
new file mode 100644
index 0000000..97da573
--- /dev/null
+++ b/ElectromagneticScattering/scattererTmatrix.pro
@@ -0,0 +1,267 @@
+///////////////////////////////
+// Author : Guillaume Demesy //
+// scattererTmatrix.pro //
+///////////////////////////////
+
+Include "scattererTmatrix_data.geo";
+myDir = "resT/";
+
+Group {
+ 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}];
+ SurfDirichlet = Region[{20}];
+}
+
+Function{
+ 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[] = Sqrt[X[]^2+Y[]^2+Z[]^2];
+ R2D[] = Sqrt[X[]^2+Y[]^2];
+ cos_theta[] = Z[]/(R3D[]+avoid_sing);
+ theta[] = Acos[cos_theta[]];
+ phi[] = Atan2[-Y[],-X[]]+Pi;
+ cart2sph[] = Tensor[X[]/R3D[], Y[]/R3D[], Z[]/R3D[],
+ X[]*Z[]/(R3D[]*R2D[]), Y[]*Z[]/(R3D[]*R2D[]), -(X[]^2+Y[]^2)/(R3D[]*R2D[]),
+ -Y[]/R2D[], X[]/R2D[],0];
+ R[] = Vector[X[]/R3D[], Y[]/R3D[], Z[]/R3D[]];
+
+ a_pml = 1.;
+ b_pml = -siwt;
+
+ sr[] = Complex[a_pml,b_pml];
+ sphi[] = Complex[1,0];
+ stheta[] = Complex[1,0];
+ r_tild[] = r_pml_in + (R3D[] - r_pml_in) * sr[];
+ theta_tild[] = theta[];
+ pml_tens_temp1[] = TensorDiag[(r_tild[]/R3D[])^2 * sphi[]*stheta[]/sr[],
+ (r_tild[]/R3D[]) * sr[]*stheta[]/sphi[],
+ (r_tild[]/R3D[]) * 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[];
+
+ 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[3,ne,me,XYZ[],k_Out];
+ Nnm_out~{pe}[] = Nnm[3,ne,me,XYZ[],k_Out];
+ Xnm~{pe}[] = Xnm[ne,me,X[],Y[],Z[]];
+ Ynm~{pe}[] = Ynm[ne,me,X[],Y[],Z[]];
+ Znm~{pe}[] = Znm[ne,me,X[],Y[],Z[]];
+ Ynm_r~{pe}[] = Ynm[ne,me,X[],Y[],Z[]] * R[];
+ source_M~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
+ source_N~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
+ SphHankelOutgoing_n~{pe}[] = (JnSph[ne ,k_Out*r_pml_in]- siwt * I[]*YnSph[ne ,k_Out*r_pml_in]) ;
+ SphHankelOutgoing_nm1~{pe}[] = (JnSph[ne-1,k_Out*r_pml_in]- siwt * I[]*YnSph[ne-1,k_Out*r_pml_in]) ;
+ dRicattiBessel~{pe}[] = (k_Out * r_pml_in * (SphHankelOutgoing_nm1~{pe}[]-((ne+1)/((k_Out*r_pml_in))) * SphHankelOutgoing_n~{pe}[]) + SphHankelOutgoing_n~{pe}[]);
+ normalize_fhnm_X~{pe}[] = 1/SphHankelOutgoing_n~{pe}[];
+ normalize_fenm_Z~{pe}[] = k_Out*r_pml_in/dRicattiBessel~{pe}[];
+ normalize_fenm_Y~{pe}[] = k_Out*r_pml_in/(SphHankelOutgoing_n~{pe}[]*Sqrt[ne*(ne+1)]);
+ EndFor
+}
+
+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 ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int_1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line2 ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 12 ; }
+ { GeoElement Triangle2 ; NumberOfPoints 12 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 17 ; }
+ { GeoElement Tetrahedron2 ; NumberOfPoints 17 ; }
+ }
+ }
+ }
+ }
+}
+
+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; }
+ If (is_FEM_o2==1)
+ { NameOfCoef un3; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un4; EntityType FacetsOf ; NameOfConstraint Dirichlet; }
+ { NameOfCoef un5; EntityType EdgesOf ; NameOfConstraint Dirichlet; }
+ EndIf
+ }
+ }
+}
+
+Formulation {
+ {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;}
+ }
+ }
+}
+
+Resolution {
+ { Name res_VPWall_helmholtz_vector;
+ System {
+ { Name T; NameOfFormulation VPWMN_helmholtz_vector; Type ComplexValue; }
+ }
+ Operation {
+ CreateDir[Str[myDir]];
+ 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
+ }
+ }
+}
+
+PostProcessing {
+ 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 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 feM~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fenm_Z~{po}[]* {u}*Conj[Znm~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+ { Name fhM~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fhnm_X~{po}[]* {u}*Conj[Xnm~{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 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; } } }
+ For po In {1:p_max}
+ { Name feN~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fenm_Z~{po}[]* {u}*Conj[Znm~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+ { Name fhN~{pe}~{po} ; Value { Integral { [ (1/r_pml_in)^2 *normalize_fhnm_X~{po}[]* {u}*Conj[Xnm~{po}[]] ]; In SurfInt; Integration Int_1 ; Jacobian JSur; } } }
+ EndFor
+ }
+ }
+ EndFor
+}
+
+PostOperation {
+ For pe In {1:p_max}
+ {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ;
+ Operation {
+ If (flag_plotcuts==1)
+ Print [ E_scat , OnGrid
+ {(r_pml_in-1*nm)*Sin[$B]*Cos[$C], (r_pml_in-1*nm)*Sin[$B]*Sin[$C], (r_pml_in-1*nm)*Cos[$B]}
+ {(r_pml_in-1*nm), {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_scat_onsphere_cart_M",Sprintf["%g.pos",pe]]],
+ Name StrCat["E_scat_onsphere_cart_M",Sprintf["%g",pe]]];
+ EndIf
+ For po In {1:p_max}
+ Print[feM~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["feM_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
+ Print[fhM~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["fhM_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
+ EndFor
+ }
+ }
+ {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
+ Operation {
+ If (flag_plotcuts==1)
+ Print [ E_scat , OnGrid
+ {(r_pml_in-1*nm)*Sin[$B]*Cos[$C], (r_pml_in-1*nm)*Sin[$B]*Sin[$C], (r_pml_in-1*nm)*Cos[$B]}
+ {(r_pml_in-1*nm), {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_scat_onsphere_cart_N",Sprintf["%g.pos",pe]]],
+ Name StrCat["E_scat_onsphere_cart_N",Sprintf["%g",pe]]];
+ EndIf
+ For po In {1:p_max}
+ Print[feN~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["feN_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
+ Print[fhN~{pe}~{po}[SurfInt] , OnGlobal, Format Table, File StrCat[myDir,StrCat[StrCat["fhN_",Sprintf["pe%g",pe]],Sprintf["po%g.dat",po]]]];
+ EndFor
+ }
+ }
+ EndFor
+}
+
+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}];
+
--
GitLab