.

include/sctl/boundary_quadrature.hpp

@@ -3148,7 +3148,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           Vector<ElemBasis> B = compute_B(sigma, alpha, beta);
           Vector<ElemBasis> gradB = compute_dB(sigma, alpha, beta);
 
-          auto compute_v = [&S,&area_elem,&toroidal_flux] () {
+          auto compute_v = [&S,&area_elem,&toroidal_flux] (const Vector<ElemBasis>& X) {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
@@ -3171,7 +3171,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
 
             Vector<ElemBasis> v(Nelem * COORD_DIM);
             Vector<ElemBasis> dX;
-            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
+            ElemBasis::Grad(dX, X);
             for (Long k = 0; k < S.Nsurf(); k++) {
               for (Long i_ = 0; i_ < S.NTor(k)*S.NPol(k); i_++) {
                 Long i = S.ElemDsp(k) + i_;
@@ -3214,7 +3214,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
             Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
             return compute_grad_adj(BxGv)*(-1.0);
@@ -3224,7 +3224,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
 
             Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
@@ -3246,7 +3246,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             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 = compute_v(S.GetElemList().ElemVector());
             Vector<ElemBasis> v_dot_GnxB = compute_dot_prod(v,GnxB);
 
             Vector<ElemBasis> dphi_dnu(Nelem);
@@ -3265,12 +3265,10 @@ template <class Real, Integer ORDER=10> class Stellarator {
             Vector<ElemBasis> nxB = compute_AxB(normal,B);
             EvalQuadrature(GnxB, S.quadrature_FxU, S, nxB, S.Laplace_FxU);
 
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> dGnxB = compute_v(GnxB);
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             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;
@@ -3279,9 +3277,9 @@ template <class Real, Integer ORDER=10> class Stellarator {
                   dv_dnu1[i][j] += -GnxB[i*COORD_DIM+2][j] * v[i*COORD_DIM+2][j] * 2 * H[i][j];
 
                   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]);
+                  dv_dnu2[i][j] += -dGnxB[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
+                  dv_dnu2[i][j] += -dGnxB[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
+                  dv_dnu2[i][j] += -dGnxB[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
                 }
               }
             }
@@ -3295,7 +3293,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             Vector<ElemBasis> nxB = compute_AxB(normal,B);
             EvalQuadrature(dGnxB, S.quadrature_FxdU, S, nxB, S.Laplace_FxdU);
 
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             Vector<ElemBasis> dphi_dnu(Nelem);
             for (Long i = 0; i < Nelem; i++) {
               for (Long j = 0; j < Nnodes; j++) {
@@ -3323,7 +3321,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             Vector<ElemBasis> nxB = compute_AxB(normal,B);
 
             Vector<ElemBasis> dGv;
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             EvalQuadrature(dGv, S.quadrature_FxdU, S, v, S.Laplace_FxdU);
 
             Vector<ElemBasis> dphi_dnu(Nelem);
@@ -3351,7 +3349,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
             Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
 
@@ -3380,7 +3378,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
             Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
 
@@ -3398,7 +3396,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             const Long Nnodes = ElemBasis::Size();
 
             Vector<ElemBasis> Gv;
-            Vector<ElemBasis> v = compute_v();
+            Vector<ElemBasis> v = compute_v(S.GetElemList().ElemVector());
             EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU);
             Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
 
@@ -3452,6 +3450,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
         Vector<Long> NtNp;
         NtNp.PushBack(20);
         NtNp.PushBack(4);
+        NtNp.PushBack(20);
+        NtNp.PushBack(4);
         S = Stellarator<Real,ORDER>(NtNp);
       }
       if (S.Nsurf() == 1) flux_pol = 0.0;