.

include/sctl/boundary_quadrature.hpp

@@ -2954,7 +2954,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
         }
 
-        // dg_dnu = (B*B) 2H
+        // dg_dnu0 = (B*B) 2H
         for (Long i = 0; i < Nelem; i++) {
           for (Long j = 0; j < Nnodes; j++) {
             dg_dnu0[i][j] = 0;
@@ -2964,7 +2964,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
         }
 
-        // dg_dnu1 = (2 B) \cdot (n \cdnot nabla) \nabla G[sigma]
+        // dg_dnu1 = (2 B) \cdot (n \cdnot \nabla) \nabla G[sigma]
         Vector<ElemBasis> d2Gsigma;
         Quadrature<Real> quadrature_Fxd2U;
         quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
@@ -2986,7 +2986,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
         }
 
-        // dg_dnu2 = (2 B) \alpha dB0_dnu \hat{\theta}
+        // dg_dnu2 = (2 B) \alpha (n \cdot \nabla) \hat{\theta} / |r|
         Vector<ElemBasis> dB0 = compute_dB0(alpha);
         for (Long i = 0; i < Nelem; i++) {
           for (Long j = 0; j < Nnodes; j++) {
@@ -3402,7 +3402,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
         }
 
-        // dAdnu3 = (v (\nabla D)^T [u n]
+        // dAdnu3 = (v n \cdot \nabla D[u]
         Vector<ElemBasis> nablaDt_u_n;
         Quadrature<Real> quadrature_dUxD;
         quadrature_dUxD.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
@@ -5220,11 +5220,26 @@ template <class Real, Integer ORDER=10> class Stellarator {
           vtu.AddElems(S.GetElemList(), dg_dnu, ORDER);
           vtu.WriteVTK("dg_dnu", comm);
         }
+        { // Save data
+          Matrix<Real> M(S.NtNp_[0]*ORDER, S.NtNp_[1]*ORDER);
+          for (Long tt = 0; tt < S.NtNp_[0]; tt++) {
+            for (Long pp = 0; pp < S.NtNp_[1]; pp++) {
+              for (Long t = 0; t < ORDER; t++) {
+                for (Long p = 0; p < ORDER; p++) {
+                  Long elem_idx = tt * S.NtNp_[1] + pp;
+                  Long node_idx = p * ORDER + t;
+                  M[tt*ORDER+t][pp*ORDER+p] = dg_dnu[elem_idx][node_idx];
+                }
+              }
+            }
+          }
+          M.Write("dg_dnu.mat");
+        }
         { // filter dg_dnu and write VTU
           const Long Nelem = S.GetElemList().NElem();
           const Long Nnodes = ElemBasis::Size();
           const Integer INTERP_ORDER = 12;
-          Long Nt = S.NtNp_[0]*ORDER/3, Np = S.NtNp_[1]*ORDER/3;
+          Long Nt = S.NtNp_[0]*ORDER/5, Np = S.NtNp_[1]*ORDER/5;
 
           Matrix<Real> M(Nt, Np); M = 0;
           const auto& quad_wts = ElemBasis::QuadWts();
@@ -5327,6 +5342,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
             vtu.AddElems(S.GetElemList(), f, ORDER);
             vtu.WriteVTK("dg_dnu_filtered", comm);
           }
+          dg_dnu = f;
         }
 
         auto compute_g = [&sigma,&alpha,&S,&area_elem,&normal,&compute_norm_area_elem,&compute_invA,&compute_half_n_plus_dG,&compute_B0,&compute_inner_prod,&comm] (const Vector<ElemBasis>& nu, Real eps) {