only one Solve operation

ElectromagneticScattering/scattering.pro

@@ -331,7 +331,9 @@ Formulation {
       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[])*Mnm[1,$NE,$ME,XYZ[],k_Out], {u} ]; In Scat_In; 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;
@@ -404,43 +406,33 @@ Resolution {
       { Name res_VPWall_helmholtz_vector;
       System {
         { Name A; NameOfFormulation VPWM_helmholtz_vector_test; Type ComplexValue; }
-        { Name B; NameOfFormulation VPWN_helmholtz_vector_test; Type ComplexValue; }
+        // { Name B; NameOfFormulation VPWN_helmholtz_vector_test; Type ComplexValue; }
       }
       Operation {
         CreateDir[Str[myDir]];
         Evaluate[Python[]{"scattering_init.py"}];
-
+        
+        Evaluate[$isN=0];
         Evaluate[ $PE = 1 ];
         Evaluate[ $NE = Floor[Sqrt[$PE]] ];
         Evaluate[ $ME = $NE*($NE+1) - Floor[$PE] ];
         Generate[A];
         Solve[A];
-        PostOperation[VPWM_postop~{1}];
-        For pe In {2:p_max}
+        // PostOperation[VPWM_postop~{1}];
+
+        For pe In {1:p_max}
+          Evaluate[$isN=0];
           Evaluate[ $PE = pe ];
           Evaluate[ $NE = Floor[Sqrt[$PE]] ];
           Evaluate[ $ME = $NE*($NE+1) - Floor[$PE] ];
           GenerateRHS[A];
           SolveAgain[A];
-          PostOperation[VPWM_postop~{pe}];
-        EndFor
-
-
-        Evaluate[ $PE = 1 ];
-        Evaluate[ $NE = Floor[Sqrt[$PE]] ];
-        Evaluate[ $ME = $NE*($NE+1) - Floor[$PE] ];
-        Generate[B];
-        Solve[B];
-        PostOperation[VPWN_postop~{1}];
-        For pe In {2:p_max}
-          Evaluate[ $PE = pe ];
-          Evaluate[ $NE = Floor[Sqrt[$PE]] ];
-          Evaluate[ $ME = $NE*($NE+1) - Floor[$PE] ];
-          GenerateRHS[B];
-          SolveAgain[B];
-          PostOperation[VPWN_postop~{pe}];
+          Test[$isN==0]{ PostOperation[VPWM_postop~{pe}]; }
+          Evaluate[$isN=1];
+          GenerateRHS[A];
+          SolveAgain[A];
+          Test[$isN==1]{ PostOperation[VPWN_postop~{pe}]; }
         EndFor
-
         Evaluate[Python[]{"scattering_post.py"}];
       }
     }
@@ -514,9 +506,16 @@ PostProcessing {
     }
   EndIf
   If (flag_study==RES_TMAT)
-    { Name VPWM_postpro_test; NameOfFormulation VPWM_helmholtz_vector_test; NameOfSystem A;
+    { Name VPWMN_postpro_test; NameOfFormulation VPWM_helmholtz_vector_test; NameOfSystem A;
             Quantity {
-          { Name E_scat            ; Value { Local { [$PE]; In All_domains; Jacobian JVol; } } } }
+          { Name E_scat            ; Value { Local { [$PE]; 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; } } }
+    }
     }
     For pe In {1:p_max}
       // { Name VPWM_postpro~{pe}; NameOfFormulation VPWM_helmholtz_vector~{pe};NameOfSystem M~{pe};
@@ -538,7 +537,7 @@ PostProcessing {
         }
       }
       // { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem N~{pe};
-      { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem B;
+      { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem A;
         Quantity {
           { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
           { Name E_scat_sph        ; Value { Local { [Vector[
@@ -625,9 +624,16 @@ PostOperation {
     }
   EndIf
   If (flag_study==RES_TMAT)
-    {Name VPWM_postop_test ; NameOfPostProcessing VPWM_postpro_test;
+    {Name VPWMN_postop_test ; NameOfPostProcessing VPWMN_postpro_test;
       Operation {
-        Print [ E_scat   , OnElementsOf Domain     , File StrCat[myDir,"Escat.pos"]];
+        // Print [ E_scat   , OnElementsOf Domain     , File StrCat[myDir,"Escat.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[StrCat["E_scat_onsphere_sph_M",Sprintf["%g",pe]],Sprintf["_r%g.dat",kr]]], Format Table];
+          EndFor
       }
     }
     For pe In {1:p_max}