.

include/sctl/boundary_quadrature.hpp

@@ -3755,6 +3755,113 @@ template <class Real, Integer ORDER=10> class Stellarator {
               }
             }
           }
+          for (Long i = 0; i < S.Nsurf(); i++) { // filter dXdt
+            const Long Nt0 = S.NTor(i);
+            const Long Np0 = S.NPol(i);
+            const Long Nelem = Nt0 * Np0;
+            const Long offset = S.ElemDsp(i);
+            const Long Nnodes = ElemBasis::Size();
+            const Long INTERP_ORDER = 12;
+            const Long Nt = Nt0*ORDER/5;
+            const Long Np = Np0*ORDER/5;
+            for (Long k = 0; k < COORD_DIM; k++) {
+              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 < Nt0; tt++) { // Set M
+                for (Long pp = 0; pp < Np0; pp++) {
+                  for (Long t = 0; t < ORDER; t++) {
+                    for (Long p = 0; p < ORDER; p++) {
+                      Real theta = (tt + Mnodes[0][t]) / Nt0;
+                      Real phi   = (pp + Mnodes[0][p]) / Np0;
+                      Long i = (Long)(theta * Nt);
+                      Long j = (Long)(phi   * Np);
+                      Real x = theta * Nt - i;
+                      Real y = phi   * Np - j;
+
+                      Long elem_idx = tt * Np0 + 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] += dXdt[(offset+elem_idx)*COORD_DIM+k][node_idx] * quad_wts[node_idx] * Interp0[ii] * Interp1[jj] / (Nt0 * Np0) * (Nt * Np);
+                        }
+                      }
+                    }
+                  }
+                }
+              }
+
+              for (Long tt = 0; tt < Nt0; tt++) {
+                for (Long pp = 0; pp < Np0; 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]) / Nt0;
+                      Real phi   = (pp + Mnodes[0][p]) / Np0;
+                      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 * Np0 + pp;
+                      Long node_idx = p * ORDER + t;
+                      dXdt[(offset+elem_idx)*COORD_DIM+k][node_idx] = f0;
+                    }
+                  }
+                }
+              }
+            }
+          }
 
           { // Update S, dt
             Stellarator<Real,ORDER> S0 = S, S1 = S, S2 = S;