use solveagain - it works but still two runs instead of one

ElectromagneticScattering/scattering.pro

@@ -331,9 +331,17 @@ 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;  }
-        // 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
+        Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Mnm[1,$NE,$ME,XYZ[],k_Out], {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
+      }
+    }
+    {Name VPWN_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;  }
+        Galerkin { [ (omega0/cel)^2*(epsilonr[]-epsilonr1[])*Nnm[1,$NE,$ME,XYZ[],k_Out], {u} ]; In Scat_In; Jacobian JVol; Integration Int_1;}
       }
     }
 
@@ -396,21 +404,44 @@ 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; }
       }
       Operation {
         CreateDir[Str[myDir]];
         Evaluate[Python[]{"scattering_init.py"}];
+
         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}
           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[Python[]{"scattering_post.py"}];
+
+
+        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}];
+        EndFor
+
+        Evaluate[Python[]{"scattering_post.py"}];
       }
     }
 
@@ -506,7 +537,8 @@ PostProcessing {
           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};
+      { Name VPWN_postpro~{pe}; NameOfFormulation VPWN_helmholtz_vector~{pe};NameOfSystem B;
         Quantity {
           { Name E_scat            ; Value { Local { [{u}]; In All_domains; Jacobian JVol; } } }
           { Name E_scat_sph        ; Value { Local { [Vector[