.

include/sctl/boundary_quadrature.hpp

@@ -2074,9 +2074,9 @@ template <class Real, Integer ORDER=10> class Stellarator {
         ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
         ValueType rinv = (r2>1e-16 ? 1/sqrt<ValueType>(r2) : 0);
         ValueType rinv3 = rinv * rinv * rinv;
-        u[0][0] =   (0) * rinv3; u[0][1] = -r[2] * rinv3; u[0][2] =  r[1] * rinv3;
-        u[1][0] =  r[2] * rinv3; u[1][1] =   (0) * rinv3; u[1][2] = -r[0] * rinv3;
-        u[2][0] = -r[1] * rinv3; u[2][1] =  r[0] * rinv3; u[2][2] =   (0) * rinv3;
+        u[0][0] =   (0) * rinv3; u[1][0] =  r[2] * rinv3; u[2][0] = -r[1] * rinv3;
+        u[0][1] = -r[2] * rinv3; u[1][1] =   (0) * rinv3; u[2][1] =  r[0] * rinv3;
+        u[0][2] =  r[1] * rinv3; u[1][2] = -r[0] * rinv3; u[2][2] =   (0) * rinv3;
       }
     };
 
@@ -2092,16 +2092,16 @@ template <class Real, Integer ORDER=10> class Stellarator {
         ValueType rinv5 = rinv2 * rinv3;
 
         u[0][0] =                                      0; u[1][0] =              - 3 * r[2] * r[0] * rinv5; u[2][0] =                3 * r[1] * r[0] * rinv5;
-        u[0][1] =                                      0; u[1][0] =              - 3 * r[2] * r[1] * rinv5; u[2][0] = -(1) * rinv3 + 3 * r[1] * r[1] * rinv5;
-        u[0][2] =                                      0; u[1][0] =  (1) * rinv3 - 3 * r[2] * r[2] * rinv5; u[2][0] =                3 * r[1] * r[2] * rinv5;
+        u[0][1] =                                      0; u[1][1] =              - 3 * r[2] * r[1] * rinv5; u[2][1] = -(1) * rinv3 + 3 * r[1] * r[1] * rinv5;
+        u[0][2] =                                      0; u[1][2] =  (1) * rinv3 - 3 * r[2] * r[2] * rinv5; u[2][2] =                3 * r[1] * r[2] * rinv5;
 
-        u[0][3] =                3 * r[2] * r[0] * rinv5; u[1][1] =                                      0; u[2][1] =  (1) * rinv3 - 3 * r[0] * r[0] * rinv5;
-        u[0][4] =                3 * r[2] * r[1] * rinv5; u[1][1] =                                      0; u[2][1] =              - 3 * r[0] * r[1] * rinv5;
-        u[0][5] = -(1) * rinv3 + 3 * r[2] * r[2] * rinv5; u[1][1] =                                      0; u[2][1] =              - 3 * r[0] * r[2] * rinv5;
+        u[0][3] =                3 * r[2] * r[0] * rinv5; u[1][3] =                                      0; u[2][3] =  (1) * rinv3 - 3 * r[0] * r[0] * rinv5;
+        u[0][4] =                3 * r[2] * r[1] * rinv5; u[1][4] =                                      0; u[2][4] =              - 3 * r[0] * r[1] * rinv5;
+        u[0][5] = -(1) * rinv3 + 3 * r[2] * r[2] * rinv5; u[1][5] =                                      0; u[2][5] =              - 3 * r[0] * r[2] * rinv5;
 
-        u[0][6] =              - 3 * r[1] * r[0] * rinv5; u[1][2] = -(1) * rinv3 + 3 * r[0] * r[0] * rinv5; u[2][2] =                                      0;
-        u[0][7] =  (1) * rinv3 - 3 * r[1] * r[1] * rinv5; u[1][2] =                3 * r[0] * r[1] * rinv5; u[2][2] =                                      0;
-        u[0][8] =              - 3 * r[1] * r[2] * rinv5; u[1][2] =                3 * r[0] * r[2] * rinv5; u[2][2] =                                      0;
+        u[0][6] =              - 3 * r[1] * r[0] * rinv5; u[1][6] = -(1) * rinv3 + 3 * r[0] * r[0] * rinv5; u[2][6] =                                      0;
+        u[0][7] =  (1) * rinv3 - 3 * r[1] * r[1] * rinv5; u[1][7] =                3 * r[0] * r[1] * rinv5; u[2][7] =                                      0;
+        u[0][8] =              - 3 * r[1] * r[2] * rinv5; u[1][8] =                3 * r[0] * r[2] * rinv5; u[2][8] =                                      0;
       }
     };
 
@@ -2351,43 +2351,46 @@ template <class Real, Integer ORDER=10> class Stellarator {
 
       Long Nelem = S.GetElemList().NElem();
       Vector<ElemBasis> Jt(Nelem*COORD_DIM), Jp(Nelem*COORD_DIM);
-      for (Long k = 0; k < S.Nsurf(); k++) {
-        Long Nt = S.NTor(k)*ORDER, Np = S.NPol(k)*ORDER;
-        const auto& X_ = S.GetElemList().ElemVector();
-        Vector<ElemBasis> X(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator<ElemBasis>)X_.begin()+S.ElemDsp(k)*COORD_DIM, false);
-
-        biest::Surface<Real> SS(Nt, Np);
-        biest::SurfaceOp<Real> surf_op(comm, Nt, Np);
-        SS.Coord() = cheb2grid(X, S.NTor(k), S.NPol(k), Nt, Np);
-
-        Vector<Real> dX, d2X;
-        surf_op.Grad2D(dX, SS.Coord());
-        surf_op.Grad2D(d2X, dX);
-
-        sctl::Vector<Real> Jt_(COORD_DIM * Nt * Np);
-        sctl::Vector<Real> Jp_(COORD_DIM * Nt * Np);
-        { // Set Jt_, Jp_
-          Vector<Real> DivV, InvLapDivV, GradInvLapDivV;
-          for (sctl::Long i = 0; i < Nt*Np; i++) { // Set V
-            for (sctl::Long k =0; k < COORD_DIM; k++) {
-              Jt_[k * Nt*Np + i] = dX[(k*2+0) * Nt*Np + i];
-              Jp_[k * Nt*Np + i] = dX[(k*2+1) * Nt*Np + i];
+      auto compute_harmonic_vector_potentials = [&S,&comm,&cheb2grid,&grid2cheb,max_iter,gmres_tol](Vector<ElemBasis>& Jt, Vector<ElemBasis>& Jp) {
+        for (Long k = 0; k < S.Nsurf(); k++) {
+          Long Nt = S.NTor(k)*ORDER, Np = S.NPol(k)*ORDER;
+          const auto& X_ = S.GetElemList().ElemVector();
+          Vector<ElemBasis> X(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator<ElemBasis>)X_.begin()+S.ElemDsp(k)*COORD_DIM, false);
+
+          biest::Surface<Real> SS(Nt, Np);
+          biest::SurfaceOp<Real> surf_op(comm, Nt, Np);
+          SS.Coord() = cheb2grid(X, S.NTor(k), S.NPol(k), Nt, Np);
+
+          Vector<Real> dX, d2X;
+          surf_op.Grad2D(dX, SS.Coord());
+          surf_op.Grad2D(d2X, dX);
+
+          sctl::Vector<Real> Jt_(COORD_DIM * Nt * Np);
+          sctl::Vector<Real> Jp_(COORD_DIM * Nt * Np);
+          { // Set Jt_, Jp_
+            Vector<Real> DivV, InvLapDivV, GradInvLapDivV;
+            for (sctl::Long i = 0; i < Nt*Np; i++) { // Set V
+              for (sctl::Long k =0; k < COORD_DIM; k++) {
+                Jt_[k * Nt*Np + i] = dX[(k*2+0) * Nt*Np + i];
+                Jp_[k * Nt*Np + i] = dX[(k*2+1) * Nt*Np + i];
+              }
             }
-          }
 
-          surf_op.SurfDiv(DivV, dX, Jt_);
-          surf_op.GradInvSurfLap(GradInvLapDivV, dX, d2X, DivV, gmres_tol * max_norm(Jt_) / max_norm(DivV), max_iter, 1.5);
-          Jt_ = Jt_ - GradInvLapDivV;
+            surf_op.SurfDiv(DivV, dX, Jt_);
+            surf_op.GradInvSurfLap(GradInvLapDivV, dX, d2X, DivV, gmres_tol * max_norm(Jt_) / max_norm(DivV), max_iter, 1.5);
+            Jt_ = Jt_ - GradInvLapDivV;
 
-          surf_op.SurfDiv(DivV, dX, Jp_);
-          surf_op.GradInvSurfLap(GradInvLapDivV, dX, d2X, DivV, gmres_tol * max_norm(Jp_) / max_norm(DivV), max_iter, 1.5);
-          Jp_ = Jp_ - GradInvLapDivV;
+            surf_op.SurfDiv(DivV, dX, Jp_);
+            surf_op.GradInvSurfLap(GradInvLapDivV, dX, d2X, DivV, gmres_tol * max_norm(Jp_) / max_norm(DivV), max_iter, 1.5);
+            Jp_ = Jp_ - GradInvLapDivV;
+          }
+          Vector<ElemBasis> Jt__(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator<ElemBasis>)Jt.begin()+S.ElemDsp(k)*COORD_DIM, false);
+          Vector<ElemBasis> Jp__(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator<ElemBasis>)Jp.begin()+S.ElemDsp(k)*COORD_DIM, false);
+          Jt__ = grid2cheb(Jt_, Nt, Np, S.NTor(k), S.NPol(k));
+          Jp__ = grid2cheb(Jp_, Nt, Np, S.NTor(k), S.NPol(k));
         }
-        Vector<ElemBasis> Jt__(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator<ElemBasis>)Jt.begin()+S.ElemDsp(k)*COORD_DIM, false);
-        Vector<ElemBasis> Jp__(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator<ElemBasis>)Jp.begin()+S.ElemDsp(k)*COORD_DIM, false);
-        Jt__ = grid2cheb(Jt_, Nt, Np, S.NTor(k), S.NPol(k));
-        Jp__ = grid2cheb(Jp_, Nt, Np, S.NTor(k), S.NPol(k));
-      }
+      };
+      compute_harmonic_vector_potentials(Jt, Jp);
 
       Vector<ElemBasis> normal, area_elem;
       auto compute_dot_prod = [](const Vector<ElemBasis>& A, const Vector<ElemBasis>& B) {
@@ -2467,16 +2470,19 @@ template <class Real, Integer ORDER=10> class Stellarator {
       };
       compute_norm_area_elem(S.GetElemList(), normal, area_elem);
 
-      S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
-      S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
-      S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
-      S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
+      S.quadrature_BS  .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavart    , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+      S.quadrature_dBS .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavartGrad, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+      S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU   , order_singular, order_direct, -1.0, comm);
+      S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU   , order_singular, order_direct, -1.0, comm);
+      S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU  , order_singular, order_direct, -1.0, comm);
+      S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF  , order_singular, order_direct, -1.0, comm);
 
-      GenericKernel<BiotSavart3D> BiotSavart;
-      GenericKernel<BiotSavartGrad3D> BiotSavartGrad;
-      Quadrature<Real> quadrature_BS, quadrature_dBS;
-      quadrature_BS .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), BiotSavart    , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
-      quadrature_dBS.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), BiotSavartGrad, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+      //Quadrature<Real> quadrature_Fxd2U;
+      //Quadrature<Real> quadrature_dUxF;
+      //Quadrature<Real> quadrature_dUxD;
+      //quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+      //quadrature_dUxF .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
+      //quadrature_dUxD .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
 
       auto compute_B0_deprecated = [&S](const Real alpha) { // alpha/|r| \hat{\theta}
         const Vector<ElemBasis> X = S.GetElemList().ElemVector();
@@ -2523,12 +2529,12 @@ template <class Real, Integer ORDER=10> class Stellarator {
             x(2) = X[i*COORD_DIM+2][j];
             Real R2inv = 1 / (x(0)*x(0) + x(1)*x(1));
 
-            dB0[(i*COORD_DIM+0)*COORD_DIM+0][j] = alpha * (2*x(0)*x(1) * R2inv*R2inv);
-            dB0[(i*COORD_DIM+0)*COORD_DIM+1][j] = alpha * (-R2inv + 2*x(1)*x(1) * R2inv*R2inv);
+            dB0[(i*COORD_DIM+0)*COORD_DIM+0][j] = -alpha * (2*x(0)*x(1) * R2inv*R2inv);
+            dB0[(i*COORD_DIM+0)*COORD_DIM+1][j] = -alpha * (-R2inv + 2*x(1)*x(1) * R2inv*R2inv);
             dB0[(i*COORD_DIM+0)*COORD_DIM+2][j] = 0;
 
-            dB0[(i*COORD_DIM+1)*COORD_DIM+0][j] = alpha * (R2inv - 2*x(0)*x(0) * R2inv*R2inv);
-            dB0[(i*COORD_DIM+1)*COORD_DIM+1][j] = alpha * (-2*x(0)*x(1) * R2inv*R2inv);
+            dB0[(i*COORD_DIM+1)*COORD_DIM+0][j] = -alpha * (R2inv - 2*x(0)*x(0) * R2inv*R2inv);
+            dB0[(i*COORD_DIM+1)*COORD_DIM+1][j] = -alpha * (-2*x(0)*x(1) * R2inv*R2inv);
             dB0[(i*COORD_DIM+1)*COORD_DIM+2][j] = 0;
 
             dB0[(i*COORD_DIM+2)*COORD_DIM+0][j] = 0;
@@ -2538,14 +2544,14 @@ template <class Real, Integer ORDER=10> class Stellarator {
         }
         return dB0;
       };
-      auto compute_B0 = [&S, &Jp,&BiotSavart,&quadrature_BS](const Real alpha) {
+      auto compute_B0 = [&S, &Jp](const Real alpha) {
         Vector<ElemBasis> B0;
-        quadrature_BS.Eval(B0, S.GetElemList(), Jp, BiotSavart);
+        S.quadrature_BS.Eval(B0, S.GetElemList(), Jp, S.BiotSavart);
         return B0*alpha;
       };
-      auto compute_dB0 = [&S, &Jp,&BiotSavartGrad,&quadrature_dBS](const Real alpha) {
+      auto compute_dB0 = [&S, &Jp](const Real alpha) {
         Vector<ElemBasis> dB0;
-        quadrature_dBS.Eval(dB0, S.GetElemList(), Jp, BiotSavartGrad);
+        S.quadrature_dBS.Eval(dB0, S.GetElemList(), Jp, S.BiotSavartGrad);
         return dB0*alpha;
       };
       auto compute_half_n_plus_dG = [&S, &normal](const Vector<ElemBasis>& sigma) { // B = n sigma/2 + dG[sigma]
@@ -2973,17 +2979,17 @@ template <class Real, Integer ORDER=10> class Stellarator {
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
               dg_dnu2[i][j] = 0;
-              dg_dnu2[i][j] += dB0[i*9+0][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+0][j];
-              dg_dnu2[i][j] += dB0[i*9+1][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+0][j];
-              dg_dnu2[i][j] += dB0[i*9+2][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+0][j];
+              dg_dnu2[i][j] -= dB0[i*9+0][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+0][j];
+              dg_dnu2[i][j] -= dB0[i*9+1][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+0][j];
+              dg_dnu2[i][j] -= dB0[i*9+2][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+0][j];
 
-              dg_dnu2[i][j] += dB0[i*9+3][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+1][j];
-              dg_dnu2[i][j] += dB0[i*9+4][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+1][j];
-              dg_dnu2[i][j] += dB0[i*9+5][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+1][j];
+              dg_dnu2[i][j] -= dB0[i*9+3][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+1][j];
+              dg_dnu2[i][j] -= dB0[i*9+4][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+1][j];
+              dg_dnu2[i][j] -= dB0[i*9+5][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+1][j];
 
-              dg_dnu2[i][j] += dB0[i*9+6][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+2][j];
-              dg_dnu2[i][j] += dB0[i*9+7][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+2][j];
-              dg_dnu2[i][j] += dB0[i*9+8][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+2][j];
+              dg_dnu2[i][j] -= dB0[i*9+6][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+2][j];
+              dg_dnu2[i][j] -= dB0[i*9+7][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+2][j];
+              dg_dnu2[i][j] -= dB0[i*9+8][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+2][j];
             }
           }
 
@@ -3254,17 +3260,17 @@ template <class Real, Integer ORDER=10> class Stellarator {
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
               Real n_n_dB0 = 0;
-              n_n_dB0 += dB0[i*9+0][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
-              n_n_dB0 += dB0[i*9+1][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j];
-              n_n_dB0 += dB0[i*9+2][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j];
+              n_n_dB0 -= dB0[i*9+0][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
+              n_n_dB0 -= dB0[i*9+1][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j];
+              n_n_dB0 -= dB0[i*9+2][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j];
 
-              n_n_dB0 += dB0[i*9+3][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j];
-              n_n_dB0 += dB0[i*9+4][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
-              n_n_dB0 += dB0[i*9+5][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j];
+              n_n_dB0 -= dB0[i*9+3][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j];
+              n_n_dB0 -= dB0[i*9+4][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
+              n_n_dB0 -= dB0[i*9+5][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j];
 
-              n_n_dB0 += dB0[i*9+6][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j];
-              n_n_dB0 += dB0[i*9+7][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j];
-              n_n_dB0 += dB0[i*9+8][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
+              n_n_dB0 -= dB0[i*9+6][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j];
+              n_n_dB0 -= dB0[i*9+7][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j];
+              n_n_dB0 -= dB0[i*9+8][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
 
               dAdnu[i][j] = u[i*Nnodes+j] * n_n_dB0;
             }
@@ -3844,17 +3850,17 @@ template <class Real, Integer ORDER=10> class Stellarator {
             for (Long i = 0; i < Nelem; i++) {
               for (Long j = 0; j < Nnodes; j++) {
                 Real dphi_dnu_ = 0;
-                dphi_dnu_ += nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
-                dphi_dnu_ += nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j];
-                dphi_dnu_ += nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j];
 
-                dphi_dnu_ += nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j];
-                dphi_dnu_ += nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
-                dphi_dnu_ += nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j];
 
-                dphi_dnu_ += nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j];
-                dphi_dnu_ += nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j];
-                dphi_dnu_ += nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j];
+                dphi_dnu_ += -nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
                 dphi_dnu[i][j] = dphi_dnu_;
               }
             }
@@ -4294,7 +4300,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           dg_dnu = f;
         }
 
-        auto compute_g = [&compute_gvec,&sigma,&alpha,&S,&area_elem,&normal,&compute_norm_area_elem,&compute_invA,&compute_half_n_plus_dG,&compute_B0,&compute_inner_prod,&comm] (const Vector<ElemBasis>& nu, Real eps) {
+        auto compute_g = [&compute_gvec,&sigma,&alpha,&S,&area_elem,&normal,&Jt,&Jp,&compute_harmonic_vector_potentials,&compute_norm_area_elem,&compute_invA,&compute_half_n_plus_dG,&compute_B0,&compute_inner_prod,&comm] (const Vector<ElemBasis>& nu, Real eps) {
           const Long Nelem = S.GetElemList().NElem();
           const Long Nnodes = ElemBasis::Size();
 
@@ -4309,7 +4315,9 @@ template <class Real, Integer ORDER=10> class Stellarator {
               S.Elem(i,2)[j] += eps*nu[i][j] * normal[i*COORD_DIM+2][j];
             }
           }
+          compute_harmonic_vector_potentials(Jt, Jp);
           compute_norm_area_elem(S.GetElemList(), normal, area_elem);
+          S.quadrature_BS  .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavart  , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
           S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
           S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
           S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
@@ -4328,7 +4336,9 @@ template <class Real, Integer ORDER=10> class Stellarator {
               S.Elem(i,2)[j] = X_orig[i*COORD_DIM+2][j];
             }
           }
+          compute_harmonic_vector_potentials(Jt, Jp);
           compute_norm_area_elem(S.GetElemList(), normal, area_elem);
+          S.quadrature_BS  .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.BiotSavart  , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
           S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
           S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
           S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
@@ -4451,20 +4461,21 @@ template <class Real, Integer ORDER=10> class Stellarator {
       Y = R * sctl::sin(theta);
     }
 
-    GenericKernel<BiotSavart3D  > BiotSavart  ;
-    GenericKernel<Laplace3D_FxU > Laplace_FxU ;
-    GenericKernel<Laplace3D_DxU > Laplace_DxU ;
-    GenericKernel<Laplace3D_FxdU> Laplace_FxdU;
-    GenericKernel<Laplace3D_dUxF> Laplace_dUxF;
-    GenericKernel<Laplace3D_Fxd2U> Laplace_Fxd2U;
-    GenericKernel<Laplace3D_dUxD> Laplace_dUxD;
-    GenericKernel<Laplace3D_DxdU> Laplace_DxdU;
+    GenericKernel<BiotSavart3D>     BiotSavart ;
+    GenericKernel<BiotSavartGrad3D> BiotSavartGrad;
+    GenericKernel<Laplace3D_FxU >   Laplace_FxU ;
+    GenericKernel<Laplace3D_DxU >   Laplace_DxU ;
+    GenericKernel<Laplace3D_FxdU>   Laplace_FxdU;
+    GenericKernel<Laplace3D_dUxF>   Laplace_dUxF;
+    GenericKernel<Laplace3D_dUxD>   Laplace_dUxD;
+    GenericKernel<Laplace3D_Fxd2U>  Laplace_Fxd2U;
+
+    Quadrature<Real> quadrature_BS  ;
+    Quadrature<Real> quadrature_dBS ;
     Quadrature<Real> quadrature_FxU ;
     Quadrature<Real> quadrature_DxU ;
     Quadrature<Real> quadrature_FxdU;
     Quadrature<Real> quadrature_dUxF;
-    Quadrature<Real> quadrature_Fxd2U;
-    Quadrature<Real> quadrature_dUxD;
 
     ElemLst elements;
     Vector<Long> NtNp_;