.

include/sctl/boundary_quadrature.hpp

@@ -3141,33 +3141,45 @@ 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_B,&compute_dB] (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, bool toroidal_flux) {
           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] () {
+          auto compute_v = [&S,&area_elem,&toroidal_flux] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
+            Real scal[2];
+            if (S.Nsurf() == 1) {
+              SCTL_ASSERT(toroidal_flux == true);
+              scal[0] = 1.0 / S.NTor(0);
+              scal[1] = 0;
+            } else if (S.Nsurf() == 2) {
+              if (toroidal_flux == true) {
+                scal[0] = -1.0 / S.NTor(0);
+                scal[1] = 1.0 / S.NTor(1);
+              } else {
+                scal[0] = 1.0 / S.NPol(0);
+                scal[1] = -1.0 / S.NPol(1);
+              }
+            } else {
+              SCTL_ASSERT(false);
+            }
+
             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;
+            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_;
+                for (Long j = 0; j < Nnodes; j++) {
+                  Real s = scal[k] / area_elem[i][j];
+                  v[i*COORD_DIM+0][j] = dX[i*COORD_DIM*2+0+(toroidal_flux?1:0)][j] * s;
+                  v[i*COORD_DIM+1][j] = dX[i*COORD_DIM*2+2+(toroidal_flux?1:0)][j] * s;
+                  v[i*COORD_DIM+2][j] = dX[i*COORD_DIM*2+4+(toroidal_flux?1:0)][j] * s;
                 }
               }
             }
@@ -3245,7 +3257,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             return dphi_dnu;
           };
-          auto compute_dphi_dnu3 = [&S,&normal,&area_elem,&H,&compute_AxB,&B] () {
+          auto compute_dphi_dnu3 = [&S,&normal,&area_elem,&H,&compute_AxB,&compute_v,&B] () {
             const Long Nelem = S.NElem();
             const Long Nnodes = ElemBasis::Size();
 
@@ -3253,6 +3265,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> dv_dnu1(Nelem), dv_dnu2(Nelem);
             { // Set dv_dnu1, dv_dnu2
               Vector<ElemBasis> dX, dGnxB;
@@ -3261,9 +3274,9 @@ template <class Real, Integer ORDER=10> class Stellarator {
               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_dnu1[i][j] += -GnxB[i*COORD_DIM+0][j] * v[i*COORD_DIM+0][j] * 2 * H[i][j];
+                  dv_dnu1[i][j] += -GnxB[i*COORD_DIM+1][j] * v[i*COORD_DIM+1][j] * 2 * H[i][j];
+                  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]);
@@ -3418,8 +3431,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
 
         { // Set dg_dnu -= dg_dsigma invA dA_dnu sigma
           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];
+          if (S.Nsurf() >= 1) dg_dnu -= compute_u_dAdnu_v_1(sigma, alpha, beta,  true) * dg_dsigma_invA[Nelem*Nnodes+0];
+          if (S.Nsurf() >= 2) dg_dnu -= compute_u_dAdnu_v_1(sigma, alpha, beta, false) * dg_dsigma_invA[Nelem*Nnodes+1];
         }
         return dg_dnu;
       };