.

include/sctl/boundary_quadrature.hpp

@@ -243,6 +243,11 @@ template <class Real, Integer DIM, Integer ORDER> class Basis {
       for (Long i = 0; i < Size(); i++) X[i] *= a;
       return X;
     }
+    Basis operator+(Real a) const {
+      Basis X = (*this);
+      for (Long i = 0; i < Size(); i++) X[i] += a;
+      return X;
+    }
     Basis& operator+=(const Basis& X) {
       for (Long i = 0; i < Size(); i++) (*this)[i] += X[i];
       return *this;
@@ -259,6 +264,10 @@ template <class Real, Integer DIM, Integer ORDER> class Basis {
       for (Long i = 0; i < Size(); i++) (*this)[i] *= a;
       return *this;
     }
+    Basis& operator+=(Real a) {
+      for (Long i = 0; i < Size(); i++) (*this)[i] += a;
+      return *this;
+    }
     Basis& operator=(Real a) {
       for (Long i = 0; i < Size(); i++) (*this)[i] = a;
       return *this;
@@ -2541,7 +2550,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
               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] = sum;
           }
         }
         return grad_adj_V;
@@ -3368,7 +3377,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             u_B[i*COORD_DIM+2][j] = u[i][j] * B[i*COORD_DIM+2][j];
           }
         }
-        dAdnu0 = compute_grad_adj(u_B);
+        dAdnu0 = compute_grad_adj(u_B)*(-1.0);
 
         // dAdnu1 = (2H) v (I/2 + \nabla G)^T [u n]
         Quadrature<Real> quadrature_dUxF;
@@ -3452,7 +3461,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             u_B0[i*COORD_DIM+2][j] = u[i*Nnodes+j] * B0[i*COORD_DIM+2][j];
           }
         }
-        dAdnu += compute_grad_adj(u_B0);
+        dAdnu -= compute_grad_adj(u_B0);
 
         return dAdnu;
       };
@@ -3527,7 +3536,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
           Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
 
-          return compute_grad_adj(BxGv);
+          return compute_grad_adj(BxGv)*(-1.0);
         };
         auto compute_dphi_dnu1 = [&S,&normal,&compute_AxB,&compute_v,&compute_half_n_plus_dG,&compute_dot_prod,sigma] () {
           const Long Nelem = S.GetElemList().NElem();
@@ -4069,7 +4078,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           Vector<ElemBasis> B = compute_B0(1.0);
           Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
 
-          return compute_grad_adj(BxGv);
+          return compute_grad_adj(BxGv)*(-1.0);
         };
         auto compute_dphi_dnu2 = [&S,&normal,&compute_AxB,&compute_v,&compute_B0,&compute_dot_prod] () {
           const Long Nelem = S.GetElemList().NElem();
@@ -4426,7 +4435,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           Vector<ElemBasis> B = compute_half_n_plus_dG(sigma);
           Vector<ElemBasis> BxGv = compute_AxB(B,Gv);
 
-          return compute_grad_adj(BxGv);
+          return compute_grad_adj(BxGv)*(-1.0);
         };
         auto compute_dphi_dnu1 = [&S,&normal,&compute_AxB,&compute_v,&compute_half_n_plus_dG,&compute_dot_prod,sigma] () {
           const Long Nelem = S.GetElemList().NElem();