.

include/sctl/boundary_quadrature.hpp

@@ -2804,51 +2804,10 @@ template <class Real, Integer ORDER=10> class Stellarator {
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-      //Real alpha, beta;
-      //Vector<ElemBasis> sigma;
-      //compute_invA(sigma, alpha, beta, flux_tor, flux_pol);
-      //Vector<ElemBasis> B = compute_B(sigma, alpha, beta);
-      //{ // Write VTU
-      //  VTUData vtu;
-      //  vtu.AddElems(S.GetElemList(), sigma, ORDER);
-      //  vtu.WriteVTK("sigma", comm);
-      //}
-      //{ // Write VTU
-      //  VTUData vtu;
-      //  vtu.AddElems(S.GetElemList(), B, ORDER);
-      //  vtu.WriteVTK("B", comm);
-      //}
-
       Real alpha, beta;
       Vector<ElemBasis> sigma;
       compute_invA(sigma, alpha, beta, flux_tor, flux_pol);
-      Vector<ElemBasis> B;
-      { // Set B
-        Vector<ElemBasis> normal, area_elem;
-        compute_norm_area_elem(S, normal, area_elem);
-
-        auto compute_B0 = [&S, &Jp](const Real alpha) {
-          Vector<ElemBasis> B0;
-          S.quadrature_BS.Eval(B0, S.GetElemList(), Jp, S.BiotSavart);
-          return B0*alpha;
-        };
-        auto compute_half_n_plus_dG = [&S, &normal](const Vector<ElemBasis>& sigma) { // B = n sigma/2 + dG[sigma]
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          Vector<ElemBasis> B;
-          S.quadrature_FxdU.Eval(B, S.GetElemList(), sigma, S.Laplace_FxdU);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              for (Long k = 0; k < COORD_DIM; k++) {
-                B[i*COORD_DIM+k][j] -= 0.5*sigma[i][j]*normal[i*COORD_DIM+k][j];
-              }
-            }
-          }
-          return B;
-        };
-        B = compute_half_n_plus_dG(sigma) + compute_B0(alpha);
-      }
+      Vector<ElemBasis> B = compute_B(sigma, alpha, beta);
       { // Write VTU
         VTUData vtu;
         vtu.AddElems(S.GetElemList(), sigma, ORDER);