.

include/sctl/boundary_quadrature.hpp

@@ -2623,56 +2623,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
         return B;
       };
 
-      auto compute_A21adj = [&S,&area_elem,&normal](Vector<Real>& A21adj_flux, Real flux) {
-        const Long Nelem = S.GetElemList().NElem();
-        const Long Nnodes = ElemBasis::Size();
-
-        Vector<ElemBasis> density(Nelem * COORD_DIM);
-        { // Set density
-          Vector<ElemBasis> dX;
-          ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              Tensor<Real,true,COORD_DIM,2> dx;
-              dx(0,0) = dX[i*COORD_DIM*2+0][j];
-              dx(0,1) = dX[i*COORD_DIM*2+1][j];
-              dx(1,0) = dX[i*COORD_DIM*2+2][j];
-              dx(1,1) = dX[i*COORD_DIM*2+3][j];
-              dx(2,0) = dX[i*COORD_DIM*2+4][j];
-              dx(2,1) = dX[i*COORD_DIM*2+5][j];
-
-              Real s = 1 / (area_elem[i][j] * S.NtNp_[0]);
-              for (Long k = 0; k < COORD_DIM; k++) {
-                density[i*COORD_DIM+k][j] = dx(k,1) * s;
-              }
-            }
-          }
-        }
-
-        Vector<ElemBasis> Gdensity;
-        S.quadrature_FxU.Eval(Gdensity, S.GetElemList(), density, S.Laplace_FxU);
-
-        Vector<ElemBasis> nxGdensity(Nelem * COORD_DIM);
-        for (Long i = 0; i < Nelem; i++) { // Set nxGdensity
-          for (Long j = 0; j < Nnodes; j++) {
-            Tensor<Real,true,COORD_DIM> Gdensity_, n;
-            Gdensity_(0) = Gdensity[i*COORD_DIM+0][j];
-            Gdensity_(1) = Gdensity[i*COORD_DIM+1][j];
-            Gdensity_(2) = Gdensity[i*COORD_DIM+2][j];
-
-            n(0) = normal[i*COORD_DIM+0][j];
-            n(1) = normal[i*COORD_DIM+1][j];
-            n(2) = normal[i*COORD_DIM+2][j];
-
-            nxGdensity[i*COORD_DIM+0][j] = n(1) * Gdensity_(2) - n(2) * Gdensity_(1);
-            nxGdensity[i*COORD_DIM+1][j] = n(2) * Gdensity_(0) - n(0) * Gdensity_(2);
-            nxGdensity[i*COORD_DIM+2][j] = n(0) * Gdensity_(1) - n(1) * Gdensity_(0);
-          }
-        }
-        S.quadrature_dUxF.Eval(A21adj_flux, S.GetElemList(), nxGdensity, S.Laplace_dUxF);
-        A21adj_flux *= flux;
-      };
-
       auto compute_A11 = [&S,&normal,&compute_half_n_plus_dG,&compute_dot_prod](Vector<Real>& B_dot_n, const Vector<Real>& sigma) {
         const Long Nelem = S.GetElemList().NElem();
         const Long Nnodes = ElemBasis::Size();
@@ -2701,7 +2651,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
         }
       };
-      auto compute_A21 = [&S,&normal,&compute_half_n_plus_dG,&compute_poloidal_circulation_,&compute_A21adj,&compute_inner_prod](const Vector<Real>& sigma) {
+      auto compute_A21 = [&S,&normal,&compute_half_n_plus_dG,&compute_poloidal_circulation_](const Vector<Real>& sigma) {
         const Long Nelem = S.GetElemList().NElem();
         const Long Nnodes = ElemBasis::Size();
 
@@ -2712,41 +2662,40 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
         }
 
-        if (0) {
-          Vector<ElemBasis> A21_(Nelem);
-          Vector<Real> A21(Nelem*Nnodes);
-          compute_A21adj(A21, 1);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              A21_[i][j] = A21[i*Nnodes+j];
-            }
-          }
-          return compute_inner_prod(A21_, sigma_);
-        } else {
-          Vector<ElemBasis> B = compute_half_n_plus_dG(sigma_);
+        if (0) { // alternate implementation
+          //Vector<ElemBasis> A21_(Nelem);
+          //Vector<Real> A21(Nelem*Nnodes);
+          //compute_A21adj(A21, 1);
+          //for (Long i = 0; i < Nelem; i++) {
+          //  for (Long j = 0; j < Nnodes; j++) {
+          //    A21_[i][j] = A21[i*Nnodes+j];
+          //  }
+          //}
+          //return compute_inner_prod(A21_, sigma_);
+        }
+        Vector<ElemBasis> B = compute_half_n_plus_dG(sigma_);
 
-          Vector<ElemBasis> J(Nelem * COORD_DIM);
-          for (Long i = 0; i < Nelem; i++) { // Set J
-            for (Long j = 0; j < Nnodes; j++) {
-              Tensor<Real,true,COORD_DIM> b, n;
-              b(0) = B[i*COORD_DIM+0][j];
-              b(1) = B[i*COORD_DIM+1][j];
-              b(2) = B[i*COORD_DIM+2][j];
+        Vector<ElemBasis> J(Nelem * COORD_DIM);
+        for (Long i = 0; i < Nelem; i++) { // Set J
+          for (Long j = 0; j < Nnodes; j++) {
+            Tensor<Real,true,COORD_DIM> b, n;
+            b(0) = B[i*COORD_DIM+0][j];
+            b(1) = B[i*COORD_DIM+1][j];
+            b(2) = B[i*COORD_DIM+2][j];
 
-              n(0) = normal[i*COORD_DIM+0][j];
-              n(1) = normal[i*COORD_DIM+1][j];
-              n(2) = normal[i*COORD_DIM+2][j];
+            n(0) = normal[i*COORD_DIM+0][j];
+            n(1) = normal[i*COORD_DIM+1][j];
+            n(2) = normal[i*COORD_DIM+2][j];
 
-              J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1);
-              J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2);
-              J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0);
-            }
+            J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1);
+            J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2);
+            J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0);
           }
-
-          Vector<ElemBasis> A;
-          S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
-          return compute_poloidal_circulation_(A)/S.NtNp_[0];
         }
+
+        Vector<ElemBasis> A;
+        S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
+        return compute_poloidal_circulation_(A)/S.NtNp_[0];
       };
       auto compute_A22 = [&S,&compute_B0,&normal,&compute_poloidal_circulation_](const Real alpha) {
         const Long Nelem = S.GetElemList().NElem();
@@ -2860,6 +2809,55 @@ template <class Real, Integer ORDER=10> class Stellarator {
         }
         return compute_inner_prod(A12_sigma, sigma);
       };
+      auto compute_A21adj = [&S,&area_elem,&normal](Vector<Real>& A21adj_flux, Real flux) {
+        const Long Nelem = S.GetElemList().NElem();
+        const Long Nnodes = ElemBasis::Size();
+
+        Vector<ElemBasis> density(Nelem * COORD_DIM);
+        { // Set density
+          Vector<ElemBasis> dX;
+          ElemBasis::Grad(dX, S.GetElemList().ElemVector());
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              Tensor<Real,true,COORD_DIM,2> dx;
+              dx(0,0) = dX[i*COORD_DIM*2+0][j];
+              dx(0,1) = dX[i*COORD_DIM*2+1][j];
+              dx(1,0) = dX[i*COORD_DIM*2+2][j];
+              dx(1,1) = dX[i*COORD_DIM*2+3][j];
+              dx(2,0) = dX[i*COORD_DIM*2+4][j];
+              dx(2,1) = dX[i*COORD_DIM*2+5][j];
+
+              Real s = 1 / (area_elem[i][j] * S.NtNp_[0]);
+              for (Long k = 0; k < COORD_DIM; k++) {
+                density[i*COORD_DIM+k][j] = dx(k,1) * s;
+              }
+            }
+          }
+        }
+
+        Vector<ElemBasis> Gdensity;
+        S.quadrature_FxU.Eval(Gdensity, S.GetElemList(), density, S.Laplace_FxU);
+
+        Vector<ElemBasis> nxGdensity(Nelem * COORD_DIM);
+        for (Long i = 0; i < Nelem; i++) { // Set nxGdensity
+          for (Long j = 0; j < Nnodes; j++) {
+            Tensor<Real,true,COORD_DIM> Gdensity_, n;
+            Gdensity_(0) = Gdensity[i*COORD_DIM+0][j];
+            Gdensity_(1) = Gdensity[i*COORD_DIM+1][j];
+            Gdensity_(2) = Gdensity[i*COORD_DIM+2][j];
+
+            n(0) = normal[i*COORD_DIM+0][j];
+            n(1) = normal[i*COORD_DIM+1][j];
+            n(2) = normal[i*COORD_DIM+2][j];
+
+            nxGdensity[i*COORD_DIM+0][j] = n(1) * Gdensity_(2) - n(2) * Gdensity_(1);
+            nxGdensity[i*COORD_DIM+1][j] = n(2) * Gdensity_(0) - n(0) * Gdensity_(2);
+            nxGdensity[i*COORD_DIM+2][j] = n(0) * Gdensity_(1) - n(1) * Gdensity_(0);
+          }
+        }
+        S.quadrature_dUxF.Eval(A21adj_flux, S.GetElemList(), nxGdensity, S.Laplace_dUxF);
+        A21adj_flux *= flux;
+      };
       auto compute_A22adj = [&compute_A22] (const Real alpha) {
         return compute_A22(alpha);
       };
@@ -3038,7 +3036,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
       Real flux = 1.0, alpha;
-      Vector<ElemBasis> sigma;
+      Vector<ElemBasis> sigma(S.GetElemList().NElem());
       compute_invA(sigma, alpha, flux);
       Vector<ElemBasis> B = compute_half_n_plus_dG(sigma) + compute_B0(alpha);
       Real g = compute_inner_prod(B, B);
@@ -3462,19 +3460,509 @@ template <class Real, Integer ORDER=10> class Stellarator {
         const Long Nelem = S.GetElemList().NElem();
         const Long Nnodes = ElemBasis::Size();
 
-        Vector<ElemBasis> dAdnu(Nelem*Nnodes);
+        Vector<ElemBasis> dAdnu(Nelem);
         dAdnu = 0;
         return dAdnu;
       };
-      auto compute_u_dAdnu_v_11 = [&S,&comm] (const Vector<Real>& u, const Vector<Real>& v) {
+      auto compute_u_dAdnu_v_11 = [&S,&comm,&area_elem,&normal,&compute_dot_prod,&compute_grad_adj,&compute_B0,&compute_dB0] (const Vector<Real>& u, const Vector<Real>& v) {
         const Long Nelem = S.GetElemList().NElem();
         const Long Nnodes = ElemBasis::Size();
 
-        Vector<ElemBasis> dAdnu(Nelem*Nnodes);
-        dAdnu = 0;
-        return dAdnu;
+        auto compute_v = [&S,&area_elem] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> v(Nelem * COORD_DIM);
+          Vector<ElemBasis> dX;
+          ElemBasis::Grad(dX, S.GetElemList().ElemVector());
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              Tensor<Real,true,COORD_DIM,2> dx;
+              dx(0,0) = dX[i*COORD_DIM*2+0][j];
+              dx(0,1) = dX[i*COORD_DIM*2+1][j];
+              dx(1,0) = dX[i*COORD_DIM*2+2][j];
+              dx(1,1) = dX[i*COORD_DIM*2+3][j];
+              dx(2,0) = dX[i*COORD_DIM*2+4][j];
+              dx(2,1) = dX[i*COORD_DIM*2+5][j];
+
+              Real s = 1 / (area_elem[i][j] * S.NtNp_[0]);
+              for (Long k = 0; k < COORD_DIM; k++) {
+                v[i*COORD_DIM+k][j] = dx(k,1) * s;
+              }
+            }
+          }
+          return v;
+        };
+        auto compute_AxB = [&S] (const Vector<ElemBasis>& A, const Vector<ElemBasis>& B) {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> J(Nelem * COORD_DIM);
+          for (Long i = 0; i < Nelem; i++) { // Set J
+            for (Long j = 0; j < Nnodes; j++) {
+              Tensor<Real,true,COORD_DIM> a, b;
+              a(0) = A[i*COORD_DIM+0][j];
+              a(1) = A[i*COORD_DIM+1][j];
+              a(2) = A[i*COORD_DIM+2][j];
+
+              b(0) = B[i*COORD_DIM+0][j];
+              b(1) = B[i*COORD_DIM+1][j];
+              b(2) = B[i*COORD_DIM+2][j];
+
+              J[i*COORD_DIM+0][j] = a(1) * b(2) - a(2) * b(1);
+              J[i*COORD_DIM+1][j] = a(2) * b(0) - a(0) * b(2);
+              J[i*COORD_DIM+2][j] = a(0) * b(1) - a(1) * b(0);
+            }
+          }
+          return J;
+        };
+
+        auto compute_dphi_dnu0 = [&S,&normal,&compute_AxB,&compute_v,&compute_dB0] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> Gv;
+          Vector<ElemBasis> v = compute_v();
+          S.quadrature_FxU.Eval(Gv, S.GetElemList(), v, S.Laplace_FxU);
+          Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
+
+          Vector<ElemBasis> gradB = compute_dB0(1.0);
+          Vector<ElemBasis> dphi_dnu(Nelem);
+          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+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[i][j] = dphi_dnu_;
+            }
+          }
+          return dphi_dnu;
+        };
+        auto compute_dphi_dnu1 = [&S,&normal,&compute_AxB,&compute_v,&compute_B0,compute_grad_adj] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> Gv;
+          Vector<ElemBasis> v = compute_v();
+          S.quadrature_FxU.Eval(Gv, S.GetElemList(), v, S.Laplace_FxU);
+
+          Vector<ElemBasis> B = compute_B0(1.0);
+          Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
+
+          return compute_grad_adj(BxGv);
+        };
+        auto compute_dphi_dnu2 = [&S,&normal,&compute_AxB,&compute_v,&compute_B0,&compute_dot_prod] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> H(Nelem);
+          { // Set mean curvature H
+            const Vector<ElemBasis> X = S.GetElemList().ElemVector();
+            Vector<ElemBasis> dX, d2X;
+            ElemBasis::Grad(dX, X);
+            ElemBasis::Grad(d2X, dX);
+            for (Long i = 0; i < Nelem; i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                Tensor<Real,true,2,2> I, invI, II;
+                for (Long k0 = 0; k0 < 2; k0++) {
+                  for (Long k1 = 0; k1 < 2; k1++) {
+                    I(k0,k1) = 0;
+                    I(k0,k1) += dX[(i*COORD_DIM+0)*2+k0][j] * dX[(i*COORD_DIM+0)*2+k1][j];
+                    I(k0,k1) += dX[(i*COORD_DIM+1)*2+k0][j] * dX[(i*COORD_DIM+1)*2+k1][j];
+                    I(k0,k1) += dX[(i*COORD_DIM+2)*2+k0][j] * dX[(i*COORD_DIM+2)*2+k1][j];
+
+                    II(k0,k1) = 0;
+                    II(k0,k1) += d2X[(i*COORD_DIM+0)*4+k0*2+k1][j] * normal[i*COORD_DIM+0][j];
+                    II(k0,k1) += d2X[(i*COORD_DIM+1)*4+k0*2+k1][j] * normal[i*COORD_DIM+1][j];
+                    II(k0,k1) += d2X[(i*COORD_DIM+2)*4+k0*2+k1][j] * normal[i*COORD_DIM+2][j];
+                  }
+                }
+                { // Set invI
+                  Real detI = I(0,0)*I(1,1)-I(0,1)*I(1,0);
+                  invI(0,0) = I(1,1) / detI;
+                  invI(0,1) = -I(0,1) / detI;
+                  invI(1,0) = -I(1,0) / detI;
+                  invI(1,1) = I(0,0) / detI;
+                }
+                { // Set H
+                  H[i][j] = 0;
+                  H[i][j] += -0.5 * II(0,0)*invI(0,0);
+                  H[i][j] += -0.5 * II(0,1)*invI(0,1);
+                  H[i][j] += -0.5 * II(1,0)*invI(1,0);
+                  H[i][j] += -0.5 * II(1,1)*invI(1,1);
+                }
+              }
+            }
+          }
+
+          Vector<ElemBasis> Gv;
+          Vector<ElemBasis> v = compute_v();
+          S.quadrature_FxU.Eval(Gv, S.GetElemList(), v, S.Laplace_FxU);
+
+          Vector<ElemBasis> B = compute_B0(1.0);
+          Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
+          Vector<ElemBasis> n_dot_BxGv = compute_dot_prod(normal,BxGv);
+
+          Vector<ElemBasis> dphi_dnu(Nelem);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              dphi_dnu[i][j] = n_dot_BxGv[i][j] * 2*H[i][j];
+            }
+          }
+          return dphi_dnu;
+        };
+        auto compute_dphi_dnu3 = [&S,&normal,&compute_AxB,&compute_v,&compute_B0,&compute_dot_prod] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> H(Nelem);
+          { // Set mean curvature H
+            const Vector<ElemBasis> X = S.GetElemList().ElemVector();
+            Vector<ElemBasis> dX, d2X;
+            ElemBasis::Grad(dX, X);
+            ElemBasis::Grad(d2X, dX);
+            for (Long i = 0; i < Nelem; i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                Tensor<Real,true,2,2> I, invI, II;
+                for (Long k0 = 0; k0 < 2; k0++) {
+                  for (Long k1 = 0; k1 < 2; k1++) {
+                    I(k0,k1) = 0;
+                    I(k0,k1) += dX[(i*COORD_DIM+0)*2+k0][j] * dX[(i*COORD_DIM+0)*2+k1][j];
+                    I(k0,k1) += dX[(i*COORD_DIM+1)*2+k0][j] * dX[(i*COORD_DIM+1)*2+k1][j];
+                    I(k0,k1) += dX[(i*COORD_DIM+2)*2+k0][j] * dX[(i*COORD_DIM+2)*2+k1][j];
+
+                    II(k0,k1) = 0;
+                    II(k0,k1) += d2X[(i*COORD_DIM+0)*4+k0*2+k1][j] * normal[i*COORD_DIM+0][j];
+                    II(k0,k1) += d2X[(i*COORD_DIM+1)*4+k0*2+k1][j] * normal[i*COORD_DIM+1][j];
+                    II(k0,k1) += d2X[(i*COORD_DIM+2)*4+k0*2+k1][j] * normal[i*COORD_DIM+2][j];
+                  }
+                }
+                { // Set invI
+                  Real detI = I(0,0)*I(1,1)-I(0,1)*I(1,0);
+                  invI(0,0) = I(1,1) / detI;
+                  invI(0,1) = -I(0,1) / detI;
+                  invI(1,0) = -I(1,0) / detI;
+                  invI(1,1) = I(0,0) / detI;
+                }
+                { // Set H
+                  H[i][j] = 0;
+                  H[i][j] += -0.5 * II(0,0)*invI(0,0);
+                  H[i][j] += -0.5 * II(0,1)*invI(0,1);
+                  H[i][j] += -0.5 * II(1,0)*invI(1,0);
+                  H[i][j] += -0.5 * II(1,1)*invI(1,1);
+                }
+              }
+            }
+          }
+
+          Vector<ElemBasis> GnxB;
+          Vector<ElemBasis> B = compute_B0(1.0);
+          Vector<ElemBasis> nxB = compute_AxB(normal,B);
+          S.quadrature_FxU.Eval(GnxB, S.GetElemList(), nxB, S.Laplace_FxU);
+
+          Vector<ElemBasis> v = compute_v();
+          Vector<ElemBasis> v_dot_GnxB = compute_dot_prod(v,GnxB);
+
+          Vector<ElemBasis> dphi_dnu(Nelem);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              dphi_dnu[i][j] = v_dot_GnxB[i][j] * 2*H[i][j];
+            }
+          }
+          return dphi_dnu;
+        };
+        auto compute_dphi_dnu4 = [&S,&normal,&area_elem,&compute_AxB,&compute_B0] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> GnxB;
+          Vector<ElemBasis> B = compute_B0(1.0);
+          Vector<ElemBasis> nxB = compute_AxB(normal,B);
+          S.quadrature_FxU.Eval(GnxB, S.GetElemList(), nxB, S.Laplace_FxU);
+
+          Vector<ElemBasis> H(Nelem);
+          { // Set mean curvature H
+            const Vector<ElemBasis> X = S.GetElemList().ElemVector();
+            Vector<ElemBasis> dX, d2X;
+            ElemBasis::Grad(dX, X);
+            ElemBasis::Grad(d2X, dX);
+            for (Long i = 0; i < Nelem; i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                Tensor<Real,true,2,2> I, invI, II;
+                for (Long k0 = 0; k0 < 2; k0++) {
+                  for (Long k1 = 0; k1 < 2; k1++) {
+                    I(k0,k1) = 0;
+                    I(k0,k1) += dX[(i*COORD_DIM+0)*2+k0][j] * dX[(i*COORD_DIM+0)*2+k1][j];
+                    I(k0,k1) += dX[(i*COORD_DIM+1)*2+k0][j] * dX[(i*COORD_DIM+1)*2+k1][j];
+                    I(k0,k1) += dX[(i*COORD_DIM+2)*2+k0][j] * dX[(i*COORD_DIM+2)*2+k1][j];
+
+                    II(k0,k1) = 0;
+                    II(k0,k1) += d2X[(i*COORD_DIM+0)*4+k0*2+k1][j] * normal[i*COORD_DIM+0][j];
+                    II(k0,k1) += d2X[(i*COORD_DIM+1)*4+k0*2+k1][j] * normal[i*COORD_DIM+1][j];
+                    II(k0,k1) += d2X[(i*COORD_DIM+2)*4+k0*2+k1][j] * normal[i*COORD_DIM+2][j];
+                  }
+                }
+                { // Set invI
+                  Real detI = I(0,0)*I(1,1)-I(0,1)*I(1,0);
+                  invI(0,0) = I(1,1) / detI;
+                  invI(0,1) = -I(0,1) / detI;
+                  invI(1,0) = -I(1,0) / detI;
+                  invI(1,1) = I(0,0) / detI;
+                }
+                { // Set H
+                  H[i][j] = 0;
+                  H[i][j] += -0.5 * II(0,0)*invI(0,0);
+                  H[i][j] += -0.5 * II(0,1)*invI(0,1);
+                  H[i][j] += -0.5 * II(1,0)*invI(1,0);
+                  H[i][j] += -0.5 * II(1,1)*invI(1,1);
+                }
+              }
+            }
+          }
+
+          Vector<ElemBasis> dv_dnu1(Nelem), dv_dnu2(Nelem), dv_dnu3(Nelem);
+          { // Set dv_dnu1, dv_dnu2, dv_dnu3
+            Vector<ElemBasis> dX, dn, V_tmp(Nelem);
+            ElemBasis::Grad(dn, normal);
+            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
+            for (Long i = 0; i < Nelem; i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                dv_dnu1[i][j] = 0;
+                dv_dnu1[i][j] += -GnxB[i*COORD_DIM+0][j] * dX[(i*COORD_DIM+0)*2+1][j] * 2 * H[i][j] / (area_elem[i][j] * S.NtNp_[0]);
+                dv_dnu1[i][j] += -GnxB[i*COORD_DIM+1][j] * dX[(i*COORD_DIM+1)*2+1][j] * 2 * H[i][j] / (area_elem[i][j] * S.NtNp_[0]);
+                dv_dnu1[i][j] += -GnxB[i*COORD_DIM+2][j] * dX[(i*COORD_DIM+2)*2+1][j] * 2 * H[i][j] / (area_elem[i][j] * S.NtNp_[0]);
+
+                dv_dnu2[i][j] = 0;
+                dv_dnu2[i][j] += GnxB[i*COORD_DIM+0][j] * dn[(i*COORD_DIM+0)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
+                dv_dnu2[i][j] += GnxB[i*COORD_DIM+1][j] * dn[(i*COORD_DIM+1)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
+                dv_dnu2[i][j] += GnxB[i*COORD_DIM+2][j] * dn[(i*COORD_DIM+2)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
+
+                V_tmp[i][j] = 0;
+                V_tmp[i][j] += GnxB[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j] / (area_elem[i][j] * S.NtNp_[0]); //dnu[i*2+1][j] *
+                V_tmp[i][j] += GnxB[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j] / (area_elem[i][j] * S.NtNp_[0]); //dnu[i*2+1][j] *
+                V_tmp[i][j] += GnxB[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j] / (area_elem[i][j] * S.NtNp_[0]); //dnu[i*2+1][j] *
+              }
+            }
+
+            { // dv_dnu3 <-- grad_adj V_tmp
+              Vector<ElemBasis> eye(Nnodes), Mgrad;
+              eye = 0;
+              for (Long i = 0; i < Nnodes; i++) eye[i][i] = 1;
+              ElemBasis::Grad(Mgrad, eye);
+
+              Vector<ElemBasis> grad_adj_V(Nelem);
+              const auto& quad_wts = ElemBasis::QuadWts();
+              for (Long i = 0; i < Nelem; i++) {
+                for (Long j = 0; j < Nnodes; j++) {
+                  Real sum = 0;
+                  for (Long k = 0; k < Nnodes; k++) {
+                    sum += Mgrad[j*2+1][k] * V_tmp[i][k] * (area_elem[i][k] * quad_wts[k]) / (quad_wts[j] * area_elem[i][j]);
+                  }
+                  dv_dnu3[i][j] = sum;
+                }
+              }
+            }
+          }
+
+          return dv_dnu1+dv_dnu2+dv_dnu3;
+        };
+        auto compute_dphi_dnu5 = [&S,&normal,&compute_AxB,&compute_v,&compute_B0] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> dGnxB;
+          Vector<ElemBasis> B = compute_B0(1.0);
+          Vector<ElemBasis> nxB = compute_AxB(normal,B);
+          S.quadrature_FxdU.Eval(dGnxB, S.GetElemList(), nxB, S.Laplace_FxdU);
+
+          Vector<ElemBasis> v = compute_v();
+          Vector<ElemBasis> dphi_dnu(Nelem);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              Real dphi_dnu_ = 0;
+              dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGnxB[(i*COORD_DIM+0)*COORD_DIM+0][j] * v[i*COORD_DIM+0][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGnxB[(i*COORD_DIM+0)*COORD_DIM+1][j] * v[i*COORD_DIM+0][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGnxB[(i*COORD_DIM+0)*COORD_DIM+2][j] * v[i*COORD_DIM+0][j];
+
+              dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGnxB[(i*COORD_DIM+1)*COORD_DIM+0][j] * v[i*COORD_DIM+1][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGnxB[(i*COORD_DIM+1)*COORD_DIM+1][j] * v[i*COORD_DIM+1][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGnxB[(i*COORD_DIM+1)*COORD_DIM+2][j] * v[i*COORD_DIM+1][j];
+
+              dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGnxB[(i*COORD_DIM+2)*COORD_DIM+0][j] * v[i*COORD_DIM+2][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGnxB[(i*COORD_DIM+2)*COORD_DIM+1][j] * v[i*COORD_DIM+2][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGnxB[(i*COORD_DIM+2)*COORD_DIM+2][j] * v[i*COORD_DIM+2][j];
+              dphi_dnu[i][j] = dphi_dnu_;
+            }
+          }
+          return dphi_dnu;
+        };
+        auto compute_dphi_dnu6 = [&S,&normal,&compute_AxB,&compute_v,&compute_B0] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> B = compute_B0(1.0);
+          Vector<ElemBasis> nxB = compute_AxB(normal,B);
+
+          Vector<ElemBasis> dGv;
+          Vector<ElemBasis> v = compute_v();
+          S.quadrature_FxdU.Eval(dGv, S.GetElemList(), v, S.Laplace_FxdU);
+
+          Vector<ElemBasis> dphi_dnu(Nelem);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              Real dphi_dnu_ = 0;
+              dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGv[(i*COORD_DIM+0)*COORD_DIM+0][j] * nxB[i*COORD_DIM+0][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGv[(i*COORD_DIM+0)*COORD_DIM+1][j] * nxB[i*COORD_DIM+0][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGv[(i*COORD_DIM+0)*COORD_DIM+2][j] * nxB[i*COORD_DIM+0][j];
+
+              dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGv[(i*COORD_DIM+1)*COORD_DIM+0][j] * nxB[i*COORD_DIM+1][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGv[(i*COORD_DIM+1)*COORD_DIM+1][j] * nxB[i*COORD_DIM+1][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGv[(i*COORD_DIM+1)*COORD_DIM+2][j] * nxB[i*COORD_DIM+1][j];
+
+              dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGv[(i*COORD_DIM+2)*COORD_DIM+0][j] * nxB[i*COORD_DIM+2][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGv[(i*COORD_DIM+2)*COORD_DIM+1][j] * nxB[i*COORD_DIM+2][j];
+              dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGv[(i*COORD_DIM+2)*COORD_DIM+2][j] * nxB[i*COORD_DIM+2][j];
+              dphi_dnu[i][j] = dphi_dnu_;
+            }
+          }
+          return dphi_dnu;
+        };
+
+        auto dphi_dnu0 = compute_dphi_dnu0();
+        auto dphi_dnu1 = compute_dphi_dnu1();
+        auto dphi_dnu2 = compute_dphi_dnu2();
+        auto dphi_dnu3 = compute_dphi_dnu3();
+        auto dphi_dnu4 = compute_dphi_dnu4();
+        auto dphi_dnu5 = compute_dphi_dnu5();
+        auto dphi_dnu6 = compute_dphi_dnu6();
+
+        return (dphi_dnu0+dphi_dnu1+dphi_dnu2+dphi_dnu3+dphi_dnu4+dphi_dnu5+dphi_dnu6) * (u[Nelem*Nnodes] * v[Nelem*Nnodes]);
       };
 
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+      if (0) {
+        const Long Nelem = S.GetElemList().NElem();
+        const Long Nnodes = ElemBasis::Size();
+
+        auto compute_v = [&S,&area_elem] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> v(Nelem * COORD_DIM);
+          Vector<ElemBasis> dX;
+          ElemBasis::Grad(dX, S.GetElemList().ElemVector());
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              Tensor<Real,true,COORD_DIM,2> dx;
+              dx(0,0) = dX[i*COORD_DIM*2+0][j];
+              dx(0,1) = dX[i*COORD_DIM*2+1][j];
+              dx(1,0) = dX[i*COORD_DIM*2+2][j];
+              dx(1,1) = dX[i*COORD_DIM*2+3][j];
+              dx(2,0) = dX[i*COORD_DIM*2+4][j];
+              dx(2,1) = dX[i*COORD_DIM*2+5][j];
+
+              Real s = 1 / (area_elem[i][j] * S.NtNp_[0]);
+              for (Long k = 0; k < COORD_DIM; k++) {
+                v[i*COORD_DIM+k][j] = dx(k,1) * s;
+              }
+            }
+          }
+          return v;
+        };
+        auto compute_AxB = [&S] (const Vector<ElemBasis>& A, const Vector<ElemBasis>& B) {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> J(Nelem * COORD_DIM);
+          for (Long i = 0; i < Nelem; i++) { // Set J
+            for (Long j = 0; j < Nnodes; j++) {
+              Tensor<Real,true,COORD_DIM> a, b;
+              a(0) = A[i*COORD_DIM+0][j];
+              a(1) = A[i*COORD_DIM+1][j];
+              a(2) = A[i*COORD_DIM+2][j];
+
+              b(0) = B[i*COORD_DIM+0][j];
+              b(1) = B[i*COORD_DIM+1][j];
+              b(2) = B[i*COORD_DIM+2][j];
+
+              J[i*COORD_DIM+0][j] = a(1) * b(2) - a(2) * b(1);
+              J[i*COORD_DIM+1][j] = a(2) * b(0) - a(0) * b(2);
+              J[i*COORD_DIM+2][j] = a(0) * b(1) - a(1) * b(0);
+            }
+          }
+          return J;
+        };
+        auto compute_dphi_dnu0 = [&S,&normal,&compute_AxB,&compute_v,&compute_dB0] () {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+          return 0;
+        };
+
+        Vector<ElemBasis> nu(Nelem);
+        nu = area_elem;
+
+        Real dphi_dnu0 = compute_inner_prod(nu, compute_dphi_dnu0());
+        std::cout<<dphi_dnu0<<' ';
+        std::cout<<'\n';
+
+        auto compute_flux = [&S,&comm,&normal,&area_elem,&compute_norm_area_elem,&compute_B0,&compute_AxB,&compute_v,&compute_inner_prod] (const Vector<ElemBasis>& nu, Real eps) {
+          const Long Nelem = S.GetElemList().NElem();
+          const Long Nnodes = ElemBasis::Size();
+
+          Vector<ElemBasis> X_orig(Nelem*COORD_DIM);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              X_orig[i*COORD_DIM+0][j] = S.Elem(i,0)[j];
+              X_orig[i*COORD_DIM+1][j] = S.Elem(i,1)[j];
+              X_orig[i*COORD_DIM+2][j] = S.Elem(i,2)[j];
+              S.Elem(i,0)[j] += eps*nu[i][j] * normal[i*COORD_DIM+0][j];
+              S.Elem(i,1)[j] += eps*nu[i][j] * normal[i*COORD_DIM+1][j];
+              S.Elem(i,2)[j] += eps*nu[i][j] * normal[i*COORD_DIM+2][j];
+            }
+          }
+          compute_norm_area_elem(normal, area_elem);
+          S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
+
+          Vector<ElemBasis> v = compute_v();
+          Vector<ElemBasis> B = compute_B0(1.0);
+          Vector<ElemBasis> J = compute_AxB(normal,B);
+          Vector<ElemBasis> A;
+          S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
+          Real flux = compute_inner_prod(v, A);
+
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              S.Elem(i,0)[j] = X_orig[i*COORD_DIM+0][j];
+              S.Elem(i,1)[j] = X_orig[i*COORD_DIM+1][j];
+              S.Elem(i,2)[j] = X_orig[i*COORD_DIM+2][j];
+            }
+          }
+          compute_norm_area_elem(normal, area_elem);
+          S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
+
+          return flux;
+        };
+        Real dphi_dnu = (compute_flux(nu,1e-3)-compute_flux(nu,-1e-3)) / 2e-3;
+        std::cout<<"dphi_dnu = "<<dphi_dnu<<'\n';
+
+        Real phi = compute_flux(nu,0);
+        std::cout<<"phi = "<<phi<<'\n';
+      }
+
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
       auto compute_Av = [&S,&area_elem,&normal,&compute_norm_area_elem,&compute_A,&comm] (const Vector<Real>& v, const Vector<ElemBasis>& nu, Real eps) {
         const Long Nelem = S.GetElemList().NElem();
         const Long Nnodes = ElemBasis::Size();
@@ -3665,7 +4153,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
         const Long Nnodes = ElemBasis::Size();
 
         Vector<ElemBasis> nu(Nelem);
-        nu = 1; //area_elem;
+        nu = area_elem;
 
         if (1) {
           Real dg_dnu0 = compute_inner_prod(nu, compute_dg_dnu(sigma, alpha, B));