.

include/sctl/boundary_quadrature.hpp

@@ -2998,31 +2998,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
         ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
         ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-        auto compute_B0 = [&S, &Jp](const Real alpha) {
-          Vector<ElemBasis> B0;
-          EvalQuadrature(B0, S.quadrature_BS, S, Jp, S.BiotSavart);
-          return B0*alpha;
-        };
-        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_grad_adj = [&S,&area_elem] (const Vector<ElemBasis>& V) {
           const Long Nelem = S.NElem();
           const Long Nnodes = ElemBasis::Size();
@@ -3096,10 +3071,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
           return grad_adj_V;
         };
 
-        auto compute_u_dAdnu_v_00 = [&S,&normal,&H,&compute_B,&compute_dB,&compute_half_n_plus_dG,&compute_grad_adj] (const Vector<Real>& u_, const Vector<ElemBasis>& v, Real alpha, Real beta) {
-          const Long Nelem = S.NElem();
+        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) {
           const Long Nnodes = ElemBasis::Size();
+          const Long Nelem = S.NElem();
 
+          Vector<ElemBasis> dAdnu0(Nelem), dAdnu1(Nelem), dAdnu2(Nelem), dAdnu3(Nelem);
           Vector<ElemBasis> u(Nelem), u_n(Nelem*COORD_DIM);
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
@@ -3110,12 +3086,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
           }
 
-          Vector<ElemBasis> dAdnu0(Nelem), dAdnu1(Nelem), dAdnu2(Nelem), dAdnu3(Nelem);
-          dAdnu0 = 0;
-          dAdnu1 = 0;
-          dAdnu2 = 0;
-          dAdnu3 = 0;
-
           // dAdnu0 = u B \cdot grad_nu
           Vector<ElemBasis> B = compute_B(v, alpha, beta);
           Vector<ElemBasis> u_B(Nelem*COORD_DIM);
@@ -3128,33 +3098,33 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
           dAdnu0 = compute_grad_adj(u_B)*(-1.0);
 
-          // dAdnu1 = (2H) v (I/2 + \nabla G)^T [u n]
-          EvalQuadrature(dAdnu1, S.quadrature_dUxF, S, u_n, S.Laplace_dUxF);
+          // dAdnu1 = (u n) \cdot (n \cdnot \nabla) B
+          Vector<ElemBasis> dB = compute_dB(v, alpha, beta);
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
-              dAdnu1[i][j] += 0.5 * u_n[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
-              dAdnu1[i][j] += 0.5 * u_n[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
-              dAdnu1[i][j] += 0.5 * u_n[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
-              dAdnu1[i][j] *= -2*H[i][j] * v[i][j];
+              dAdnu1[i][j] = 0;
+              dAdnu1[i][j] -= dB[i*9+0][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+0][j];
+              dAdnu1[i][j] -= dB[i*9+1][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+1][j];
+              dAdnu1[i][j] -= dB[i*9+2][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+2][j];
+
+              dAdnu1[i][j] -= dB[i*9+3][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+0][j];
+              dAdnu1[i][j] -= dB[i*9+4][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+1][j];
+              dAdnu1[i][j] -= dB[i*9+5][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+2][j];
+
+              dAdnu1[i][j] -= dB[i*9+6][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+0][j];
+              dAdnu1[i][j] -= dB[i*9+7][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+1][j];
+              dAdnu1[i][j] -= dB[i*9+8][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+2][j];
             }
           }
 
-          // dAdnu2 = (u n) \cdot (n \cdnot \nabla) B
-          Vector<ElemBasis> d2G_v = compute_dB(v, alpha, beta);
+          // dAdnu2 = (2H) v (I/2 + \nabla G)^T [u n]
+          EvalQuadrature(dAdnu2, S.quadrature_dUxF, S, u_n, S.Laplace_dUxF);
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
-              dAdnu2[i][j] = 0;
-              dAdnu2[i][j] -= d2G_v[i*9+0][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+0][j];
-              dAdnu2[i][j] -= d2G_v[i*9+1][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+1][j];
-              dAdnu2[i][j] -= d2G_v[i*9+2][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+2][j];
-
-              dAdnu2[i][j] -= d2G_v[i*9+3][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+0][j];
-              dAdnu2[i][j] -= d2G_v[i*9+4][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+1][j];
-              dAdnu2[i][j] -= d2G_v[i*9+5][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+2][j];
-
-              dAdnu2[i][j] -= d2G_v[i*9+6][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+0][j];
-              dAdnu2[i][j] -= d2G_v[i*9+7][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+1][j];
-              dAdnu2[i][j] -= d2G_v[i*9+8][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+2][j];
+              dAdnu2[i][j] += 0.5 * u_n[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
+              dAdnu2[i][j] += 0.5 * u_n[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
+              dAdnu2[i][j] += 0.5 * u_n[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
+              dAdnu2[i][j] *= -2*H[i][j] * v[i][j];
             }
           }
 
@@ -3172,16 +3142,31 @@ template <class Real, Integer ORDER=10> class Stellarator {
 
           return dAdnu0 + dAdnu1 + dAdnu2 + dAdnu3;
         };
-        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<Real>& v) {
+
+        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> sigma(Nelem);
+          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++) {
-              sigma[i][j] = v[i*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) {
+          const Long Nnodes = ElemBasis::Size();
+          const Long Nelem = S.NElem();
 
           auto compute_v = [&S,&area_elem] () {
             const Long Nelem = S.NElem();
@@ -3460,7 +3445,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, const Vector<Real>& v) {
+        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) {
           const Long Nelem = S.NElem();
           const Long Nnodes = ElemBasis::Size();
 
@@ -3696,26 +3681,16 @@ template <class Real, Integer ORDER=10> class Stellarator {
           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]);
+          return (dphi_dnu0+dphi_dnu1+dphi_dnu2+dphi_dnu3+dphi_dnu4+dphi_dnu5+dphi_dnu6) * (u[Nelem*Nnodes] * alpha);
         };
+
         { // Set dg_dnu -= dg_dsigma invA dA_dnu sigma
-          Vector<Real> sigma_(Nelem*Nnodes+S.Nsurf());
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              sigma_[i*Nnodes+j] = sigma[i][j];
-            }
-          }
-          if (S.Nsurf() >= 1) sigma_[Nelem*Nnodes+0] = alpha;
-          if (S.Nsurf() >= 2) sigma_[Nelem*Nnodes+1] = beta;
-
-          auto dg_dnu1 = compute_u_dAdnu_v_00(dg_dsigma_invA, sigma, alpha, beta)*(-1);
-          // auto dg_dnu2 = compute_u_dAdnu_v_01(dg_dsigma_invA, sigma_)*(-1);
-          auto dg_dnu3 = compute_u_dAdnu_v_10(dg_dsigma_invA, sigma_)*(-1);
-          auto dg_dnu4 = compute_u_dAdnu_v_11(dg_dsigma_invA, sigma_)*(-1);
-          dg_dnu += dg_dnu1;
-          // dg_dnu += dg_dnu2;
-          dg_dnu += dg_dnu3;
-          dg_dnu += dg_dnu4;
+          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);
+          dg_dnu -= dg_dnu1;
+          dg_dnu -= dg_dnu3;
+          dg_dnu -= dg_dnu4;
         }
         return dg_dnu;
       };