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

order 2

parent b3257b60
No related branches found
No related tags found
No related merge requests found
Pipeline #4400 passed
Stage: test
Show whitespace changes
Inline Side-by-side
ElectromagneticScattering/scattering.geo
+ 2
− 0
View file @ 972a768c
...@@ -723,6 +723,7 @@ If(flag_cartpml==0) ...@@ -723,6 +723,7 @@ If(flag_cartpml==0)
Physical Volume(2) = {1073}; // Out - air Physical Volume(2) = {1073}; // Out - air
Physical Volume(3) = {1091}; // pml Physical Volume(3) = {1091}; // pml
Physical Point(2000) = {1}; // PrintPoint Physical Point(2000) = {1}; // PrintPoint
Physical Surface(10) = {1055,1057,1059,1062,1064,1066,1068,1070};
EndIf EndIf
// Mesh.Algorithm = 1; // // 1=MeshAdapt, 5=Delaunay, 6=Frontal // Mesh.Algorithm = 1; // // 1=MeshAdapt, 5=Delaunay, 6=Frontal
...@@ -752,3 +753,4 @@ Printf(Sprintf("n_max = %g;" ,n_max)) >> "scattering_tmp.py"; ...@@ -752,3 +753,4 @@ Printf(Sprintf("n_max = %g;" ,n_max)) >> "scattering_tmp.py";
Printf(Sprintf("p_max = %g;" ,p_max)) >> "scattering_tmp.py"; Printf(Sprintf("p_max = %g;" ,p_max)) >> "scattering_tmp.py";
Printf(Sprintf("siwt = %g;" ,siwt)) >> "scattering_tmp.py"; Printf(Sprintf("siwt = %g;" ,siwt)) >> "scattering_tmp.py";
Mesh.ElementOrder = 2;
ElectromagneticScattering/scattering.pro
+ 38
− 10
View file @ 972a768c
...@@ -21,6 +21,7 @@ Group { ...@@ -21,6 +21,7 @@ Group {
Domain = Region[{Scat_In,Scat_Out}]; Domain = Region[{Scat_In,Scat_Out}];
PMLs = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}]; PMLs = Region[{PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
All_domains = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}]; All_domains = Region[{Scat_In,Scat_Out,PMLxyz,PMLxy,PMLxz,PMLyz,PMLx,PMLy,PMLz}];
SurfInt = Region[{}];
EndIf EndIf
If(flag_cartpml==0) If(flag_cartpml==0)
Scat_In = Region[1]; Scat_In = Region[1];
...@@ -28,6 +29,7 @@ Group { ...@@ -28,6 +29,7 @@ Group {
Domain = Region[{Scat_In,Scat_Out}]; Domain = Region[{Scat_In,Scat_Out}];
PMLs = Region[3]; PMLs = Region[3];
All_domains = Region[{Scat_In,Scat_Out,PMLs}]; All_domains = Region[{Scat_In,Scat_Out,PMLs}];
SurfInt = Region[{10}];
EndIf EndIf
PrintPoint = Region[2000]; PrintPoint = Region[2000];
} }
...@@ -192,6 +194,8 @@ Function{ ...@@ -192,6 +194,8 @@ Function{
me = ne*(ne+1) - Floor[pe]; me = ne*(ne+1) - Floor[pe];
Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out]; Mnm_source~{pe}[] = Mnm[1,ne,me,XYZ[],k_Out];
Nnm_source~{pe}[] = Nnm[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_M~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm_source~{pe}[];
source_N~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[]; source_N~{pe}[] = (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm_source~{pe}[];
EndFor EndFor
...@@ -236,9 +240,12 @@ Integration { ...@@ -236,9 +240,12 @@ Integration {
{ Type Gauss ; { Type Gauss ;
Case { Case {
{ GeoElement Point ; NumberOfPoints 1 ; } { GeoElement Point ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 4 ; } { GeoElement Line2 ; NumberOfPoints 4 ; }
{ GeoElement Triangle ; NumberOfPoints 6 ; } { GeoElement Triangle2 ; NumberOfPoints 6 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 15 ; } { GeoElement Tetrahedron2 ; NumberOfPoints 15 ; }
// { GeoElement Line2 ; NumberOfPoints 4 ; }
// { GeoElement Triangle2 ; NumberOfPoints 6 ; }
// { GeoElement Tetrahedron2 ; NumberOfPoints 15 ; }
} }
} }
} }
...@@ -249,16 +256,16 @@ FunctionSpace { ...@@ -249,16 +256,16 @@ FunctionSpace {
{ Name Hcurl; Type Form1; { Name Hcurl; Type Form1;
BasisFunction { BasisFunction {
{ Name sn; NameOfCoef un; Function BF_Edge; { Name sn; NameOfCoef un; Function BF_Edge;
Support Region[All_domains]; Entity EdgesOf[All]; } Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
{ Name sn2; NameOfCoef un2; Function BF_Edge_2E; { Name sn2; NameOfCoef un2; Function BF_Edge_2E;
Support Region[All_domains]; Entity EdgesOf[All]; } Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
If (is_FEM_o2==1) If (is_FEM_o2==1)
{ Name sn3; NameOfCoef un3; Function BF_Edge_3F_b; { Name sn3; NameOfCoef un3; Function BF_Edge_3F_b;
Support Region[All_domains]; Entity FacetsOf[All]; } Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
{ Name sn4; NameOfCoef un4; Function BF_Edge_3F_c; { Name sn4; NameOfCoef un4; Function BF_Edge_3F_c;
Support Region[All_domains]; Entity FacetsOf[All]; } Support Region[{All_domains,SurfInt}]; Entity FacetsOf[All]; }
{ Name sn5; NameOfCoef un5; Function BF_Edge_4E; { Name sn5; NameOfCoef un5; Function BF_Edge_4E;
Support Region[All_domains]; Entity EdgesOf[All]; } Support Region[{All_domains,SurfInt}]; Entity EdgesOf[All]; }
EndIf EndIf
} }
Constraint { Constraint {
...@@ -360,14 +367,28 @@ Resolution { ...@@ -360,14 +367,28 @@ Resolution {
Operation { Operation {
CreateDir[Str[myDir]]; CreateDir[Str[myDir]];
Evaluate[Python[]{"scattering_init.py"}]; Evaluate[Python[]{"scattering_init.py"}];
Generate[M_1];
Solve[M_1];
PostOperation[VPWM_postop_1];
For pe In {1:p_max} For pe In {1:p_max}
Generate[M~{pe}]; GenerateRHS[M~{pe}];
Solve[M~{pe}]; Solve[M~{pe}];
PostOperation[VPWM_postop~{pe}]; PostOperation[VPWM_postop~{pe}];
Generate[N~{pe}]; Generate[N~{pe}];
Solve[N~{pe}]; Solve[N~{pe}];
PostOperation[VPWN_postop~{pe}]; PostOperation[VPWN_postop~{pe}];
EndFor 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"}]; Evaluate[Python[]{"scattering_post.py"}];
} }
} }
...@@ -452,6 +473,10 @@ PostProcessing { ...@@ -452,6 +473,10 @@ PostProcessing {
In All_domains; Jacobian JVol; } } } In All_domains; Jacobian JVol; } } }
{ Name H_scat ; Value { Local { [siwt*I[]/(mur[]*mu0*omega0)*{Curl 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; } } } { 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~{pe};NameOfSystem N~{pe}; { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem N~{pe};
...@@ -557,6 +582,9 @@ PostOperation { ...@@ -557,6 +582,9 @@ PostOperation {
{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)} }, {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]]], 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]]]; 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} ; {Name VPWN_postop~{pe}; NameOfPostProcessing VPWN_postpro~{pe} ;
......
ElectromagneticScattering/scattering_init.py
+ 2
− 0
View file @ 972a768c
...@@ -218,6 +218,8 @@ def field_VSH_expansion(post_filename): ...@@ -218,6 +218,8 @@ def field_VSH_expansion(post_filename):
orth3 = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) orth3 = np.trapz(np.trapz((np.sin(theta_sph)*ZnmconjXnm).transpose(),theta_sph[:,0]),phi_sph[0,:])
aM_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_aM_nm2 aM_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjMnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_aM_nm2
bN_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_bN_nm2 bN_nm[int(po)-1] = np.trapz(np.trapz((np.sin(theta_sph)*EdotconjNnm).transpose(),theta_sph[:,0]),phi_sph[0,:]) / normalize_bN_nm2
print(np.loadtxt('run_results/a_pe%gpo%g.dat'%(1,int(po))))
print(aM_nm)
return [fhnm_X,fenm_Y,fenm_Z,aM_nm,bN_nm,FF_Xnm_t,FF_Xnm_p,FF_erCrossXnm_t,FF_erCrossXnm_p] return [fhnm_X,fenm_Y,fenm_Z,aM_nm,bN_nm,FF_Xnm_t,FF_Xnm_p,FF_erCrossXnm_t,FF_erCrossXnm_p]
########################### ###########################
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment