.

include/sctl/boundary_quadrature.hpp

@@ -2526,31 +2526,23 @@ template <class Real, Integer ORDER=10> class Stellarator {
             dudX_V[i*2+0][j] = 0;
             dudX_V[i*2+1][j] = 0;
 
-            dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+0)*2+0][j] * V[i*COORD_DIM+0][j];
-            dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+1)*2+0][j] * V[i*COORD_DIM+1][j];
-            dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+2)*2+0][j] * V[i*COORD_DIM+2][j];
+            dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+0)*2+0][j] * V[i*COORD_DIM+0][j] * area_elem[i][j];
+            dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+1)*2+0][j] * V[i*COORD_DIM+1][j] * area_elem[i][j];
+            dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+2)*2+0][j] * V[i*COORD_DIM+2][j] * area_elem[i][j];
 
-            dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+0)*2+1][j] * V[i*COORD_DIM+0][j];
-            dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+1)*2+1][j] * V[i*COORD_DIM+1][j];
-            dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+2)*2+1][j] * V[i*COORD_DIM+2][j];
+            dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+0)*2+1][j] * V[i*COORD_DIM+0][j] * area_elem[i][j];
+            dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+1)*2+1][j] * V[i*COORD_DIM+1][j] * area_elem[i][j];
+            dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+2)*2+1][j] * V[i*COORD_DIM+2][j] * area_elem[i][j];
           }
         }
 
-        Vector<ElemBasis> eye(Nnodes), Mgrad;
-        eye = 0;
-        for (Long i = 0; i < Nnodes; i++) eye[i][i] = 1;
-        ElemBasis::Grad(Mgrad, eye);
+        Vector<ElemBasis> grad_dudX_V;
+        ElemBasis::Grad(grad_dudX_V, dudX_V);
 
         Vector<ElemBasis> grad_adj_V(Nelem);
-        const auto& quad_wts = ElemBasis::QuadWts();
         for (Long i = 0; i < Nelem; i++) {
           for (Long j = 0; j < Nnodes; j++) {
-            Real sum = 0;
-            for (Long k = 0; k < Nnodes; k++) {
-              sum += Mgrad[j*2+0][k] * dudX_V[i*2+0][k] * (area_elem[i][k] * quad_wts[k]) / (quad_wts[j] * area_elem[i][j]);
-              sum += Mgrad[j*2+1][k] * dudX_V[i*2+1][k] * (area_elem[i][k] * quad_wts[k]) / (quad_wts[j] * area_elem[i][j]);
-            }
-            grad_adj_V[i][j] = sum;
+            grad_adj_V[i][j] = -(grad_dudX_V[(i*2+0)*2+0][j] + grad_dudX_V[(i*2+1)*2+1][j]) / area_elem[i][j];
           }
         }
         return grad_adj_V;
@@ -3725,27 +3717,18 @@ template <class Real, Integer ORDER=10> class Stellarator {
                 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] / (area_elem[i][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] / (area_elem[i][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] / (area_elem[i][j] * S.NtNp_[0]); //dnu[i*2+1][j] *
+                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> eye(Nnodes), Mgrad;
-              eye = 0;
-              for (Long i = 0; i < Nnodes; i++) eye[i][i] = 1;
-              ElemBasis::Grad(Mgrad, eye);
-
-              Vector<ElemBasis> grad_adj_V(Nelem);
-              const auto& quad_wts = ElemBasis::QuadWts();
+              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++) {
-                  Real sum = 0;
-                  for (Long k = 0; k < Nnodes; k++) {
-                    sum += Mgrad[j*2+1][k] * V_tmp[i][k] * (area_elem[i][k] * quad_wts[k]) / (quad_wts[j] * area_elem[i][j]);
-                  }
-                  dv_dnu3[i][j] = sum;
+                  dv_dnu3[i][j] = -grad_adj_V_tmp[i*2+1][j] / area_elem[i][j];
                 }
               }
             }
@@ -4267,27 +4250,18 @@ template <class Real, Integer ORDER=10> class Stellarator {
                 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] / (area_elem[i][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] / (area_elem[i][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] / (area_elem[i][j] * S.NtNp_[0]); //dnu[i*2+1][j] *
+                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> eye(Nnodes), Mgrad;
-              eye = 0;
-              for (Long i = 0; i < Nnodes; i++) eye[i][i] = 1;
-              ElemBasis::Grad(Mgrad, eye);
-
-              Vector<ElemBasis> grad_adj_V(Nelem);
-              const auto& quad_wts = ElemBasis::QuadWts();
+              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++) {
-                  Real sum = 0;
-                  for (Long k = 0; k < Nnodes; k++) {
-                    sum += Mgrad[j*2+1][k] * V_tmp[i][k] * (area_elem[i][k] * quad_wts[k]) / (quad_wts[j] * area_elem[i][j]);
-                  }
-                  dv_dnu3[i][j] = sum;
+                  dv_dnu3[i][j] = -grad_adj_V_tmp[i*2+1][j] / area_elem[i][j];
                 }
               }
             }
@@ -5149,7 +5123,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
 
           Vector<Real> dg_dsigma_invA = compute_invAadj(dg_dsigma);
-
           { // Set dg_dnu = - dg_dsigma invA dA_dnu sigma
             Vector<Real> sigma_(Nelem*Nnodes+1);
             for (Long i = 0; i < Nelem; i++) {