From 9de4957cff9378e68598142e3f4579e8993b1f86 Mon Sep 17 00:00:00 2001 From: Guillaume Demesy <guillaume.demesy@fresnel.fr> Date: Mon, 2 Sep 2019 08:05:09 +0200 Subject: [PATCH] runtine variables + SolveAgain --- ElectromagneticScattering/scattering.pro | 86 ++- .../scattering_solveagain.pro | 704 ------------------ 2 files changed, 59 insertions(+), 731 deletions(-) delete mode 100644 ElectromagneticScattering/scattering_solveagain.pro diff --git a/ElectromagneticScattering/scattering.pro b/ElectromagneticScattering/scattering.pro index 18f13e4..1379de6 100644 --- a/ElectromagneticScattering/scattering.pro +++ b/ElectromagneticScattering/scattering.pro @@ -39,6 +39,7 @@ Function{ // 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; @@ -323,6 +324,19 @@ Formulation { } } EndFor + {Name VPWM_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; } + // For pe In {1:p_max} + Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm[$PE,ne,me,XYZ[],k_Out], {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;} + // EndFor + } + } + EndIf If (flag_study==RES_GREEN) For ncomp In {0:2} @@ -357,41 +371,49 @@ Resolution { } EndIf If (flag_study==RES_TMAT) - { Name res_VPWall_helmholtz_vector; + // { Name res_VPWall_helmholtz_vector; + // System { + // For pe In {1:p_max} + // { Name M~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; } + // { Name N~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; } + // EndFor + // } + // Operation { + // CreateDir[Str[myDir]]; + // Evaluate[Python[]{"scattering_init.py"}]; + + // 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"}]; + // } + // } + { Name res_VPWall_helmholtz_vector; System { - For pe In {1:p_max} - { Name M~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; } - { Name N~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe}; Type ComplexValue; Frequency Freq; } - EndFor + { Name A; NameOfFormulation VPWM_helmholtz_vector_test; Type ComplexValue; } } Operation { CreateDir[Str[myDir]]; Evaluate[Python[]{"scattering_init.py"}]; - - - // Generate[M_1]; - // Solve[M_1]; - // PostOperation[VPWM_postop_1]; - For pe In {1:p_max} - // GenerateRHS[M~{pe}]; - Solve[M~{pe}]; + Evaluate[ $PE = 1 ]; + Generate[A]; + Solve[A]; + PostOperation[VPWM_postop~{1}]; + For pe In {2:p_max} + Evaluate[ $PE = pe ]; + GenerateRHS[A]; + SolveAgain[A]; 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"}]; + // Evaluate[Python[]{"scattering_post.py"}]; } } + EndIf If (flag_study==RES_GREEN) { Name res_GreenAll_helmholtz_vector; @@ -461,8 +483,13 @@ PostProcessing { } EndIf If (flag_study==RES_TMAT) + { Name VPWM_postpro_test; NameOfFormulation VPWM_helmholtz_vector_test; NameOfSystem A; + Quantity { + { Name E_scat ; Value { Local { [$PE]; In All_domains; Jacobian JVol; } } } } + } For pe In {1:p_max} - { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe};NameOfSystem M~{pe}; + // { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe};NameOfSystem M~{pe}; + { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector_test; NameOfSystem A; Quantity { { Name E_scat ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } } { Name E_scat_sph ; Value { Local { [Vector[ @@ -566,6 +593,11 @@ PostOperation { } EndIf If (flag_study==RES_TMAT) + {Name VPWM_postop_test ; NameOfPostProcessing VPWM_postpro_test; + Operation { + Print [ E_scat , OnElementsOf Domain , File StrCat[myDir,"Escat.pos"]]; + } + } For pe In {1:p_max} {Name VPWM_postop~{pe}; NameOfPostProcessing VPWM_postpro~{pe} ; Operation { diff --git a/ElectromagneticScattering/scattering_solveagain.pro b/ElectromagneticScattering/scattering_solveagain.pro deleted file mode 100644 index dbe5e4d..0000000 --- a/ElectromagneticScattering/scattering_solveagain.pro +++ /dev/null @@ -1,704 +0,0 @@ -/////////////////////////////// -// 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 -- GitLab