.

include/sctl/boundary_quadrature.hpp

@@ -3070,22 +3070,6 @@ 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();
-
-          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();
@@ -3157,10 +3141,13 @@ template <class Real, Integer ORDER=10> class Stellarator {
 
           return dAdnu0 + dAdnu1 + dAdnu2 + dAdnu3;
         };
-        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) {
+        auto compute_u_dAdnu_v_1 = [&S,&area_elem,&normal,&H,&compute_grad_adj,&compute_B,&compute_dB] (const Vector<ElemBasis>& sigma, Real alpha, Real beta) {
           const Long Nnodes = ElemBasis::Size();
           const Long Nelem = S.NElem();
 
+          Vector<ElemBasis> B = compute_B(sigma, alpha, beta);
+          Vector<ElemBasis> gradB = compute_dB(sigma, alpha, beta);
+
           auto compute_v = [&S,&area_elem] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
@@ -3210,20 +3197,17 @@ template <class Real, Integer ORDER=10> class Stellarator {
             return J;
           };
 
-          auto compute_dphi_dnu0 = [&S,&normal,&compute_AxB,&compute_v,&compute_half_n_plus_dG,compute_grad_adj,sigma] () {
+          auto compute_dphi_dnu0 = [&S,&normal,&compute_AxB,&compute_v,&B,compute_grad_adj] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> Gv;
             Vector<ElemBasis> v = compute_v();
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
-
-            Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
             Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
-
             return compute_grad_adj(BxGv)*(-1.0);
           };
-          auto compute_dphi_dnu1 = [&S,&normal,&H,&compute_AxB,&compute_v,&compute_half_n_plus_dG,sigma] () {
+          auto compute_dphi_dnu1 = [&S,&normal,&H,&compute_AxB,&compute_v,&B] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
@@ -3231,7 +3215,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
             Vector<ElemBasis> v = compute_v();
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
 
-            Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
             Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
             Vector<ElemBasis> n_dot_BxGv = compute_dot_prod(normal,BxGv);
 
@@ -3243,12 +3226,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dphi_dnu;
           };
-          auto compute_dphi_dnu2 = [&S,&normal,&H,&compute_AxB,&compute_v,&compute_half_n_plus_dG,sigma] () {
+          auto compute_dphi_dnu2 = [&S,&normal,&H,&compute_AxB,&compute_v,&B] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> GnxB;
-            Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
             Vector<ElemBasis> nxB = compute_AxB(normal,B);
             EvalQuadrature(GnxB, S.quadrature_FxU, S, nxB, S.Laplace_FxU);
 
@@ -3263,12 +3245,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dphi_dnu;
           };
-          auto compute_dphi_dnu3 = [&S,&normal,&area_elem,&H,&compute_AxB,&compute_half_n_plus_dG,sigma] () {
+          auto compute_dphi_dnu3 = [&S,&normal,&area_elem,&H,&compute_AxB,&B] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> GnxB;
-            Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
             Vector<ElemBasis> nxB = compute_AxB(normal,B);
             EvalQuadrature(GnxB, S.quadrature_FxU, S, nxB, S.Laplace_FxU);
 
@@ -3293,12 +3274,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dv_dnu1 + dv_dnu2;
           };
-          auto compute_dphi_dnu4 = [&S,&normal,&compute_AxB,&compute_v,&compute_half_n_plus_dG,sigma] () {
+          auto compute_dphi_dnu4 = [&S,&normal,&compute_AxB,&compute_v,&B] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> dGnxB;
-            Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
             Vector<ElemBasis> nxB = compute_AxB(normal,B);
             EvalQuadrature(dGnxB, S.quadrature_FxdU, S, nxB, S.Laplace_FxdU);
 
@@ -3323,11 +3303,10 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dphi_dnu;
           };
-          auto compute_dphi_dnu5 = [&S,&normal,&compute_AxB,&compute_v,&compute_half_n_plus_dG,sigma] () {
+          auto compute_dphi_dnu5 = [&S,&normal,&compute_AxB,&compute_v,&B] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
-            Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
             Vector<ElemBasis> nxB = compute_AxB(normal,B);
 
             Vector<ElemBasis> dGv;
@@ -3354,7 +3333,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dphi_dnu;
           };
-          auto compute_dphi_dnu6 = [&S,&normal,&compute_AxB,&compute_v,sigma] () {
+          auto compute_dphi_dnu6 = [&S,&normal,&compute_AxB,&compute_v,&gradB] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
@@ -3363,8 +3342,6 @@ 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;
-            EvalQuadrature(gradB, S.quadrature_Fxd2U, S, sigma, S.Laplace_Fxd2U);
             Vector<ElemBasis> dphi_dnu(Nelem);
             for (Long i = 0; i < Nelem; i++) {
               for (Long j = 0; j < Nnodes; j++) {
@@ -3385,7 +3362,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dphi_dnu;
           };
-          auto compute_dphi_dnu7 = [&S,&normal,&H,&compute_AxB,&compute_v,sigma] () {
+          auto compute_dphi_dnu7 = [&S,&normal,&H,&compute_AxB,&compute_v,&sigma] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
@@ -3403,7 +3380,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dphi_dnu;
           };
-          auto compute_dphi_dnu8 = [&S,&normal,&H,&compute_AxB,&compute_v,sigma] () {
+          auto compute_dphi_dnu8 = [&S,&normal,&H,&compute_AxB,&compute_v,&sigma] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
@@ -3436,254 +3413,13 @@ template <class Real, Integer ORDER=10> class Stellarator {
           auto dphi_dnu7 = compute_dphi_dnu7();
           auto dphi_dnu8 = compute_dphi_dnu8();
 
-          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_dB] (const Vector<Real>& u, Real alpha, Real beta) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-
-          auto compute_v = [&S,&area_elem] () {
-            const Long Nelem = S.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.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_dB] () {
-            const Long Nelem = S.NElem();
-            const Long Nnodes = ElemBasis::Size();
-
-            Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
-            EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
-            Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
-
-            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++) {
-                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_B,compute_grad_adj] () {
-            const Long Nelem = S.NElem();
-            const Long Nnodes = ElemBasis::Size();
-
-            Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
-            EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
-
-            Vector<ElemBasis> B = compute_B(Vector<ElemBasis>(), 1.0, 0);
-            Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
-
-            return compute_grad_adj(BxGv)*(-1.0);
-          };
-          auto compute_dphi_dnu2 = [&S,&normal,&H,&compute_AxB,&compute_v,&compute_B] () {
-            const Long Nelem = S.NElem();
-            const Long Nnodes = ElemBasis::Size();
-
-            Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
-            EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
-
-            Vector<ElemBasis> B = compute_B(Vector<ElemBasis>(), 1.0, 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,&H,&compute_AxB,&compute_v,&compute_B] () {
-            const Long Nelem = S.NElem();
-            const Long Nnodes = ElemBasis::Size();
-
-            Vector<ElemBasis> GnxB;
-            Vector<ElemBasis> B = compute_B(Vector<ElemBasis>(), 1.0, 0);
-            Vector<ElemBasis> nxB = compute_AxB(normal,B);
-            EvalQuadrature(GnxB, S.quadrature_FxU, S, 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,&H,&compute_AxB,&compute_B] () {
-            const Long Nelem = S.NElem();
-            const Long Nnodes = ElemBasis::Size();
-
-            Vector<ElemBasis> GnxB;
-            Vector<ElemBasis> B = compute_B(Vector<ElemBasis>(), 1.0, 0);
-            Vector<ElemBasis> nxB = compute_AxB(normal,B);
-            EvalQuadrature(GnxB, S.quadrature_FxU, S, nxB, S.Laplace_FxU);
-
-            Vector<ElemBasis> dv_dnu1(Nelem), dv_dnu2(Nelem);
-            { // Set dv_dnu1, dv_dnu2
-              Vector<ElemBasis> dX, dGnxB;
-              ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-              ElemBasis::Grad(dGnxB, GnxB);
-              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] += -dGnxB[(i*COORD_DIM+0)*2+1][j] * normal[i*COORD_DIM+0][j] / (area_elem[i][j] * S.NtNp_[0]);
-                  dv_dnu2[i][j] += -dGnxB[(i*COORD_DIM+1)*2+1][j] * normal[i*COORD_DIM+1][j] / (area_elem[i][j] * S.NtNp_[0]);
-                  dv_dnu2[i][j] += -dGnxB[(i*COORD_DIM+2)*2+1][j] * normal[i*COORD_DIM+2][j] / (area_elem[i][j] * S.NtNp_[0]);
-                }
-              }
-            }
-            return dv_dnu1 + dv_dnu2;
-          };
-          auto compute_dphi_dnu5 = [&S,&normal,&compute_AxB,&compute_v,&compute_B] () {
-            const Long Nelem = S.NElem();
-            const Long Nnodes = ElemBasis::Size();
-
-            Vector<ElemBasis> dGnxB;
-            Vector<ElemBasis> B = compute_B(Vector<ElemBasis>(), 1.0, 0);
-            Vector<ElemBasis> nxB = compute_AxB(normal,B);
-            EvalQuadrature(dGnxB, S.quadrature_FxdU, S, 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_B] () {
-            const Long Nelem = S.NElem();
-            const Long Nnodes = ElemBasis::Size();
-
-            Vector<ElemBasis> B = compute_B(Vector<ElemBasis>(), 1.0, 0);
-            Vector<ElemBasis> nxB = compute_AxB(normal,B);
-
-            Vector<ElemBasis> dGv;
-            Vector<ElemBasis> v = compute_v();
-            EvalQuadrature(dGv, S.quadrature_FxdU, S, 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] * alpha);
+          return (dphi_dnu0+dphi_dnu1+dphi_dnu2+dphi_dnu3+dphi_dnu4+dphi_dnu5+dphi_dnu6+dphi_dnu7+dphi_dnu8);
         };
 
         { // Set dg_dnu -= dg_dsigma invA dA_dnu sigma
-          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_dnu2;
+          dg_dnu -= compute_u_dAdnu_v_0(dg_dsigma_invA, sigma, alpha, beta);
+          if (S.Nsurf() >= 1) dg_dnu -= compute_u_dAdnu_v_1(sigma, alpha, beta) * dg_dsigma_invA[Nelem*Nnodes+0];
+          // if (S.Nsurf() >= 2) dg_dnu -= compute_u_dAdnu_v_2(sigma, alpha, beta) * dg_dsigma_invA[Nelem*Nnodes+1];
         }
         return dg_dnu;
       };