.

include/sctl/boundary_quadrature.hpp

@@ -3070,8 +3070,23 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
           return grad_adj_V;
         };
+        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();
 
-        auto compute_u_dAdnu_v_00 = [&S,&normal,&H,&compute_B,&compute_dB,&compute_grad_adj] (const Vector<Real>& u_, const Vector<ElemBasis>& v, Real alpha, Real beta) {
+          Vector<ElemBasis> B;
+          EvalQuadrature(B, S.quadrature_FxdU, S, 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;
+        };
+
+        auto compute_u_dAdnu_v_0 = [&S,&normal,&H,&compute_B,&compute_dB,&compute_grad_adj] (const Vector<Real>& u_, const Vector<ElemBasis>& v, Real alpha, Real beta) {
           const Long Nnodes = ElemBasis::Size();
           const Long Nelem = S.NElem();
 
@@ -3142,29 +3157,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 
           return dAdnu0 + dAdnu1 + dAdnu2 + dAdnu3;
         };
-
-        auto compute_dB0 = [&S, &Jp](const Real alpha) {
-          Vector<ElemBasis> dB0;
-          EvalQuadrature(dB0, S.quadrature_dBS, S, Jp, S.BiotSavartGrad);
-          return dB0*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;
-          EvalQuadrature(B, S.quadrature_FxdU, S, 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;
-        };
-
-        auto compute_u_dAdnu_v_10 = [&S,&area_elem,&normal,&H,&compute_grad_adj,&compute_half_n_plus_dG] (const Vector<Real>& u, const Vector<ElemBasis>& sigma) {
+        auto compute_u_dAdnu_v_1 = [&S,&area_elem,&normal,&H,&compute_grad_adj,&compute_half_n_plus_dG] (const Vector<Real>& u, const Vector<ElemBasis>& sigma, Real alpha, Real beta) {
           const Long Nnodes = ElemBasis::Size();
           const Long Nelem = S.NElem();
 
@@ -3445,7 +3438,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 
           return (dphi_dnu0+dphi_dnu1+dphi_dnu2+dphi_dnu3+dphi_dnu4+dphi_dnu5+dphi_dnu6+dphi_dnu7+dphi_dnu8) * u[Nelem*Nnodes];
         };
-        auto compute_u_dAdnu_v_11 = [&S,&area_elem,&normal,&H,&compute_grad_adj,&compute_B,&compute_dB0] (const Vector<Real>& u, Real alpha, Real beta) {
+        auto compute_u_dAdnu_v_11 = [&S,&area_elem,&normal,&H,&compute_grad_adj,&compute_B,&compute_dB] (const Vector<Real>& u, Real alpha, Real beta) {
           const Long Nelem = S.NElem();
           const Long Nnodes = ElemBasis::Size();
 
@@ -3498,7 +3491,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             return J;
           };
 
-          auto compute_dphi_dnu0 = [&S,&normal,&compute_AxB,&compute_v,&compute_dB0] () {
+          auto compute_dphi_dnu0 = [&S,&normal,&compute_AxB,&compute_v,&compute_dB] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
@@ -3507,7 +3500,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
             Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
 
-            Vector<ElemBasis> gradB = compute_dB0(1.0);
+            Vector<ElemBasis> gradB = compute_dB(Vector<ElemBasis>(), 1.0, 0);
             Vector<ElemBasis> dphi_dnu(Nelem);
             for (Long i = 0; i < Nelem; i++) {
               for (Long j = 0; j < Nnodes; j++) {
@@ -3685,12 +3678,12 @@ template <class Real, Integer ORDER=10> class Stellarator {
         };
 
         { // Set dg_dnu -= dg_dsigma invA dA_dnu sigma
-          auto dg_dnu1 = compute_u_dAdnu_v_00(dg_dsigma_invA, sigma, alpha, beta);
-          auto dg_dnu3 = compute_u_dAdnu_v_10(dg_dsigma_invA, sigma);
-          auto dg_dnu4 = compute_u_dAdnu_v_11(dg_dsigma_invA, alpha, beta);
+          auto dg_dnu0 = compute_u_dAdnu_v_0(dg_dsigma_invA, sigma, alpha, beta);
+          auto dg_dnu1 = compute_u_dAdnu_v_1(dg_dsigma_invA, sigma, alpha, beta);
+          auto dg_dnu2 = compute_u_dAdnu_v_11(dg_dsigma_invA, alpha, beta);
+          dg_dnu -= dg_dnu0;
           dg_dnu -= dg_dnu1;
-          dg_dnu -= dg_dnu3;
-          dg_dnu -= dg_dnu4;
+          dg_dnu -= dg_dnu2;
         }
         return dg_dnu;
       };