.

include/sctl/boundary_quadrature.hpp

@@ -2476,13 +2476,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
       S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU   , order_singular, order_direct, -1.0, comm);
       S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU  , order_singular, order_direct, -1.0, comm);
       S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF  , order_singular, order_direct, -1.0, comm);
-
-      //Quadrature<Real> quadrature_Fxd2U;
-      //Quadrature<Real> quadrature_dUxF;
-      //Quadrature<Real> quadrature_dUxD;
-      //quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
-      //quadrature_dUxF .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
-      //quadrature_dUxD .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
+      s.quadrature_dUxD .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
+      s.quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
 
       auto compute_B0_deprecated = [&S](const Real alpha) { // alpha/|r| \hat{\theta}
         const Vector<ElemBasis> X = S.GetElemList().ElemVector();
@@ -3008,11 +3003,12 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
 
           // dg_dnu4 = (2H) sigma (\nabla G)^T [2 B]
-          Quadrature<Real> quadrature_dUxF;
-          quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
-          quadrature_dUxF.Eval(dg_dnu4, S.GetElemList(), v, S.Laplace_dUxF);
+          S.quadrature_dUxF.Eval(dg_dnu4, S.GetElemList(), v, S.Laplace_dUxF);
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
+              dg_dnu4[i][j] += 0.5 * v[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
+              dg_dnu4[i][j] += 0.5 * v[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
+              dg_dnu4[i][j] += 0.5 * v[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
               dg_dnu4[i][j] *= 2*H[i][j] * sigma[i][j];
             }
           }
@@ -3204,11 +3200,12 @@ template <class Real, Integer ORDER=10> class Stellarator {
           dAdnu0 = compute_grad_adj(u_B)*(-1.0);
 
           // dAdnu1 = (2H) v (I/2 + \nabla G)^T [u n]
-          Quadrature<Real> quadrature_dUxF;
-          quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
-          quadrature_dUxF.Eval(dAdnu1, S.GetElemList(), u_n, S.Laplace_dUxF);
+          S.quadrature_dUxF.Eval(dAdnu1, S.GetElemList(), u_n, S.Laplace_dUxF);
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
+              dAdnu1[i][j] += 0.5 * u_n[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
+              dAdnu1[i][j] += 0.5 * u_n[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
+              dAdnu1[i][j] += 0.5 * u_n[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
               dAdnu1[i][j] *= -2*H[i][j] * v[i][j];
             }
           }
@@ -3696,9 +3693,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             Vector<ElemBasis> nxGv = compute_AxB(Gv,normal);
 
             Vector<ElemBasis> dphi_dnu(Nelem);
-            Quadrature<Real> quadrature_dUxF;
-            quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
-            quadrature_dUxF.Eval(dphi_dnu, S.GetElemList(), nxGv, S.Laplace_dUxF);
+            S.quadrature_dUxF.Eval(dphi_dnu, S.GetElemList(), nxGv, S.Laplace_dUxF);
             for (Long i = 0; i < Nelem; i++) {
               for (Long j = 0; j < Nnodes; j++) {
                 dphi_dnu[i][j] *= -2*H[i][j] * sigma[i][j];
@@ -4348,7 +4343,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
         };
         {
           Vector<ElemBasis> nu(Nelem);
-          nu = area_elem;
+          nu = dg_dnu;
 
           Real eps = 1e-4;
           Real g0 = compute_g(nu,-eps);
@@ -4358,7 +4353,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
         }
         {
           Vector<ElemBasis> nu(Nelem);
-          nu = 1;
+          nu = area_elem;
 
           Real eps = 1e-4;
           Real g0 = compute_g(nu,-eps);
@@ -4368,7 +4363,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
         }
         {
           Vector<ElemBasis> nu(Nelem);
-          nu = dg_dnu;
+          nu = 1;
 
           Real eps = 1e-4;
           Real g0 = compute_g(nu,-eps);
@@ -4476,6 +4471,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
     Quadrature<Real> quadrature_DxU ;
     Quadrature<Real> quadrature_FxdU;
     Quadrature<Real> quadrature_dUxF;
+    Quadrature<Real> quadrature_dUxD;
+    Quadrature<Real> quadrature_Fxd2U;
 
     ElemLst elements;
     Vector<Long> NtNp_;