.

include/sctl/boundary_quadrature.hpp

@@ -5120,18 +5120,13 @@ template <class Real, Integer ORDER=10> class Stellarator {
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-      if (0) { // test grad_g
+      if (1) { // test grad_g
         const Long Nelem = S.GetElemList().NElem();
         const Long Nnodes = ElemBasis::Size();
 
         Vector<ElemBasis> dg_dnu;
         { // Compute dg_dnu
           dg_dnu = compute_dg_dnu(sigma, alpha, B);
-          { // Write VTU
-            VTUData vtu;
-            vtu.AddElems(S.GetElemList(), dg_dnu, ORDER);
-            vtu.WriteVTK("dg_dnu_", comm);
-          }
 
           Vector<Real> dg_dsigma(Nelem*Nnodes+1);
           { // Set dg_dsigma
@@ -5143,30 +5138,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             dg_dsigma[Nelem*Nnodes] = compute_dg_dalpha(B);
           }
-          { // Write VTU
-            VTUData vtu;
-            Vector<ElemBasis> dg_dsigma_(Nelem);
-            for (Long i = 0; i < Nelem; i++) {
-              for (Long j = 0; j < Nnodes; j++) {
-                dg_dsigma_[i][j] = dg_dsigma[i*Nnodes+j];
-              }
-            }
-            vtu.AddElems(S.GetElemList(), dg_dsigma_, ORDER);
-            vtu.WriteVTK("dg_dsigma", comm);
-          }
 
           Vector<Real> dg_dsigma_invA = compute_invAadj(dg_dsigma);
-          { // Write VTU
-            VTUData vtu;
-            Vector<ElemBasis> dg_dsigma_invA_(Nelem);
-            for (Long i = 0; i < Nelem; i++) {
-              for (Long j = 0; j < Nnodes; j++) {
-                dg_dsigma_invA_[i][j] = dg_dsigma_invA[i*Nnodes+j];
-              }
-            }
-            vtu.AddElems(S.GetElemList(), dg_dsigma_invA_, ORDER);
-            vtu.WriteVTK("dg_dsigma_invA", comm);
-          }
 
           { // Set dg_dnu = - dg_dsigma invA dA_dnu sigma
             Vector<Real> sigma_(Nelem*Nnodes+1);
@@ -5181,26 +5154,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
             auto dg_dnu2 = compute_u_dAdnu_v_01(dg_dsigma_invA, sigma_)*(-1);
             auto dg_dnu3 = compute_u_dAdnu_v_10(dg_dsigma_invA, sigma_)*(-1);
             auto dg_dnu4 = compute_u_dAdnu_v_11(dg_dsigma_invA, sigma_)*(-1);
-            { // Write VTU
-              VTUData vtu;
-              vtu.AddElems(S.GetElemList(), dg_dnu1, ORDER);
-              vtu.WriteVTK("dg_dnu1", comm);
-            }
-            { // Write VTU
-              VTUData vtu;
-              vtu.AddElems(S.GetElemList(), dg_dnu2, ORDER);
-              vtu.WriteVTK("dg_dnu2", comm);
-            }
-            { // Write VTU
-              VTUData vtu;
-              vtu.AddElems(S.GetElemList(), dg_dnu3, ORDER);
-              vtu.WriteVTK("dg_dnu3", comm);
-            }
-            { // Write VTU
-              VTUData vtu;
-              vtu.AddElems(S.GetElemList(), dg_dnu4, ORDER);
-              vtu.WriteVTK("dg_dnu4", comm);
-            }
 
             {
               //Vector<ElemBasis> nu(Nelem);
@@ -5267,6 +5220,114 @@ template <class Real, Integer ORDER=10> class Stellarator {
           vtu.AddElems(S.GetElemList(), dg_dnu, ORDER);
           vtu.WriteVTK("dg_dnu", comm);
         }
+        { // 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/2, Np = S.NtNp_[1]*ORDER/2;
+
+          Matrix<Real> M(Nt, Np); M = 0;
+          const auto& quad_wts = ElemBasis::QuadWts();
+          const Matrix<Real>& Mnodes = Basis<Real,1,ORDER>::Nodes();
+          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++) {
+                  Real theta = (tt + Mnodes[0][t]) / S.NtNp_[0];
+                  Real phi   = (pp + Mnodes[0][p]) / S.NtNp_[1];
+                  Long i = (Long)(theta * Nt);
+                  Long j = (Long)(phi   * Np);
+                  Real x = theta * Nt - i;
+                  Real y = phi   * Np - j;
+
+                  Long elem_idx = tt * S.NtNp_[1] + pp;
+                  Long node_idx = p * ORDER + t;
+
+                  Vector<Real> Interp0(INTERP_ORDER);
+                  Vector<Real> Interp1(INTERP_ORDER);
+                  { // Set Interp0, Interp1
+                    auto node = [] (Long i) {
+                      return (Real)i - (INTERP_ORDER-1)/2;
+                    };
+                    for (Long i = 0; i < INTERP_ORDER; i++) {
+                      Real wt_x = 1, wt_y = 1;
+                      for (Long j = 0; j < INTERP_ORDER; j++) {
+                        if (j != i) {
+                          wt_x *= (x - node(j)) / (node(i) - node(j));
+                          wt_y *= (y - node(j)) / (node(i) - node(j));
+                        }
+                        Interp0[i] = wt_x;
+                        Interp1[i] = wt_y;
+                      }
+                    }
+                  }
+
+                  for (Long ii = 0; ii < INTERP_ORDER; ii++) {
+                    for (Long jj = 0; jj < INTERP_ORDER; jj++) {
+                      Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt;
+                      Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np;
+                      M[idx_i][idx_j] += dg_dnu[elem_idx][node_idx] * area_elem[elem_idx][node_idx] * quad_wts[node_idx] * Interp0[ii] * Interp1[jj];
+                    }
+                  }
+                }
+              }
+            }
+          }
+
+          Vector<ElemBasis> f(Nelem);
+          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++) {
+                  Matrix<Real> Mnodes = Basis<Real,1,ORDER>::Nodes();
+                  Real theta = (tt + Mnodes[0][t]) / S.NtNp_[0];
+                  Real phi   = (pp + Mnodes[0][p]) / S.NtNp_[1];
+                  Long i = (Long)(theta * Nt);
+                  Long j = (Long)(phi   * Np);
+                  Real x = theta * Nt - i;
+                  Real y = phi   * Np - j;
+
+                  Vector<Real> Interp0(INTERP_ORDER);
+                  Vector<Real> Interp1(INTERP_ORDER);
+                  { // Set Interp0, Interp1
+                    auto node = [] (Long i) {
+                      return (Real)i - (INTERP_ORDER-1)/2;
+                    };
+                    for (Long i = 0; i < INTERP_ORDER; i++) {
+                      Real wt_x = 1, wt_y = 1;
+                      for (Long j = 0; j < INTERP_ORDER; j++) {
+                        if (j != i) {
+                          wt_x *= (x - node(j)) / (node(i) - node(j));
+                          wt_y *= (y - node(j)) / (node(i) - node(j));
+                        }
+                        Interp0[i] = wt_x;
+                        Interp1[i] = wt_y;
+                      }
+                    }
+                  }
+
+                  Real f0 = 0;
+                  for (Long ii = 0; ii < INTERP_ORDER; ii++) {
+                    for (Long jj = 0; jj < INTERP_ORDER; jj++) {
+                      Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt;
+                      Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np;
+                      f0 += Interp0[ii] * Interp1[jj] * M[idx_i][idx_j];
+                    }
+                  }
+
+                  Long elem_idx = tt * S.NtNp_[1] + pp;
+                  Long node_idx = p * ORDER + t;
+                  f[elem_idx][node_idx] = f0;
+                }
+              }
+            }
+          }
+          { // Write VTU
+            VTUData vtu;
+            vtu.AddElems(S.GetElemList(), f, ORDER);
+            vtu.WriteVTK("dg_dnu_filtered", comm);
+          }
+        }
 
         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) {
           const Long Nelem = S.GetElemList().NElem();