.

include/sctl/boundary_quadrature.hpp

@@ -2247,8 +2247,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
       Stellarator<Real,ORDER> S;
       { // Init S
         Vector<Long> NtNp;
-        NtNp.PushBack(40);
-        NtNp.PushBack(8);
+        NtNp.PushBack(20);
+        NtNp.PushBack(4);
         S = Stellarator<Real,ORDER>(NtNp);
       }
 
@@ -3699,11 +3699,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
           }
 
-          Vector<ElemBasis> dv_dnu1(Nelem), dv_dnu2(Nelem), dv_dnu3(Nelem);
-          { // Set dv_dnu1, dv_dnu2, dv_dnu3
-            Vector<ElemBasis> dX, dn, V_tmp(Nelem);
-            ElemBasis::Grad(dn, normal);
+          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;
@@ -3712,29 +3712,13 @@ template <class Real, Integer ORDER=10> class Stellarator {
                 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] += GnxB[i*COORD_DIM+0][j] * dn[(i*COORD_DIM+0)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
-                dv_dnu2[i][j] += GnxB[i*COORD_DIM+1][j] * dn[(i*COORD_DIM+1)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
-                dv_dnu2[i][j] += GnxB[i*COORD_DIM+2][j] * dn[(i*COORD_DIM+2)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
-
-                V_tmp[i][j] = 0;
-                V_tmp[i][j] += GnxB[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j] / S.NtNp_[0]; //dnu[i*2+1][j] *
-                V_tmp[i][j] += GnxB[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j] / S.NtNp_[0]; //dnu[i*2+1][j] *
-                V_tmp[i][j] += GnxB[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j] / S.NtNp_[0]; //dnu[i*2+1][j] *
-              }
-            }
-
-            { // dv_dnu3 <-- grad_adj V_tmp
-              Vector<ElemBasis> grad_adj_V_tmp(Nelem);
-              ElemBasis::Grad(grad_adj_V_tmp, V_tmp);
-              for (Long i = 0; i < Nelem; i++) {
-                for (Long j = 0; j < Nnodes; j++) {
-                  dv_dnu3[i][j] = -grad_adj_V_tmp[i*2+1][j] / area_elem[i][j];
-                }
+                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+dv_dnu3;
+          return dv_dnu1 + dv_dnu2;
         };
         auto compute_dphi_dnu4 = [&S,&normal,&compute_AxB,&compute_v,&compute_half_n_plus_dG,sigma] () {
           const Long Nelem = S.GetElemList().NElem();
@@ -4232,11 +4216,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
           }
 
-          Vector<ElemBasis> dv_dnu1(Nelem), dv_dnu2(Nelem), dv_dnu3(Nelem);
-          { // Set dv_dnu1, dv_dnu2, dv_dnu3
-            Vector<ElemBasis> dX, dn, V_tmp(Nelem);
-            ElemBasis::Grad(dn, normal);
+          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;
@@ -4245,29 +4229,13 @@ template <class Real, Integer ORDER=10> class Stellarator {
                 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] += GnxB[i*COORD_DIM+0][j] * dn[(i*COORD_DIM+0)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
-                dv_dnu2[i][j] += GnxB[i*COORD_DIM+1][j] * dn[(i*COORD_DIM+1)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
-                dv_dnu2[i][j] += GnxB[i*COORD_DIM+2][j] * dn[(i*COORD_DIM+2)*2+1][j] / (area_elem[i][j] * S.NtNp_[0]);
-
-                V_tmp[i][j] = 0;
-                V_tmp[i][j] += GnxB[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j] / S.NtNp_[0]; //dnu[i*2+1][j] *
-                V_tmp[i][j] += GnxB[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j] / S.NtNp_[0]; //dnu[i*2+1][j] *
-                V_tmp[i][j] += GnxB[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j] / S.NtNp_[0]; //dnu[i*2+1][j] *
-              }
-            }
-
-            { // dv_dnu3 <-- grad_adj V_tmp
-              Vector<ElemBasis> grad_adj_V_tmp(Nelem);
-              ElemBasis::Grad(grad_adj_V_tmp, V_tmp);
-              for (Long i = 0; i < Nelem; i++) {
-                for (Long j = 0; j < Nnodes; j++) {
-                  dv_dnu3[i][j] = -grad_adj_V_tmp[i*2+1][j] / area_elem[i][j];
-                }
+                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+dv_dnu3;
+          return dv_dnu1 + dv_dnu2;
         };
         auto compute_dphi_dnu5 = [&S,&normal,&compute_AxB,&compute_v,&compute_B0] () {
           const Long Nelem = S.GetElemList().NElem();
@@ -5217,7 +5185,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
           M.Write("dg_dnu.mat");
         }
-        { // filter dg_dnu and write VTU
+        if (0) { // filter dg_dnu and write VTU
           const Long Nelem = S.GetElemList().NElem();
           const Long Nnodes = ElemBasis::Size();
           const Integer INTERP_ORDER = 12;