update SCTL

include/sctl.hpp

@@ -20,6 +20,9 @@
 #define SCTL_QUOTEME_1(x) #x
 #define SCTL_INCLUDE(x) SCTL_QUOTEME(SCTL_NAMESPACE/x)
 
+// Tensor
+#include SCTL_INCLUDE(tensor.hpp)
+
 // Tree
 #include SCTL_INCLUDE(tree.hpp)
 
@@ -67,4 +70,8 @@
 #include SCTL_INCLUDE(stacktrace.h)
 const int sgh = SCTL_NAMESPACE::SetSigHandler(); // Set signal handler
 
+// Boundary quadrature, Kernel functions
+#include SCTL_INCLUDE(kernel_functions.hpp)
+#include SCTL_INCLUDE(boundary_quadrature.hpp)
+
 #endif //_SCTL_HPP_

include/sctl/boundary_quadrature.hpp

@@ -0,0 +1,1789 @@
+#ifndef _SCTL_BOUNDARY_QUADRATURE_HPP_
+#define _SCTL_BOUNDARY_QUADRATURE_HPP_
+
+#include <mutex>
+#include <atomic>
+#include <tuple>
+
+namespace SCTL_NAMESPACE {
+
+template <class Real, Integer DIM, Integer ORDER> class Basis {
+  public:
+    using ValueType = Real;
+
+    // class EvalOperator {
+    //   public:
+    // };
+    using EvalOpType = Matrix<ValueType>;
+
+    static constexpr Long Dim() {
+      return DIM;
+    }
+    static constexpr Long Size() {
+      return pow<DIM,Long>(ORDER);
+    }
+    static const Matrix<ValueType>& Nodes() {
+      static Matrix<ValueType> nodes_(DIM,Size());
+      auto nodes_1d = [](Integer i) {
+        return 0.5 - 0.5 * sctl::cos<ValueType>((2*i+1) * const_pi<ValueType>() / (2*ORDER));
+      };
+      { // Set nodes_
+        static std::mutex mutex;
+        static std::atomic<Integer> first_time(true);
+        if (first_time.load(std::memory_order_relaxed)) {
+          std::lock_guard<std::mutex> guard(mutex);
+          if (first_time.load(std::memory_order_relaxed)) {
+            Integer N = 1;
+            for (Integer d = 0; d < DIM; d++) {
+              for (Integer j = 0; j < ORDER; j++) {
+                for (Integer i = 0; i < N; i++) {
+                  for (Integer k = 0; k < d; k++) {
+                    nodes_[k][j*N+i] = nodes_[k][i];
+                  }
+                  nodes_[d][j*N+i] = nodes_1d(j);
+                }
+              }
+              N *= ORDER;
+            }
+
+            std::atomic_thread_fence(std::memory_order_seq_cst);
+            first_time.store(false);
+          }
+        }
+      }
+      return nodes_;
+    }
+
+    static void Grad(Vector<Basis>& dX, const Vector<Basis>& X) {
+      static Matrix<ValueType> GradOp[DIM];
+      static std::mutex mutex;
+      static std::atomic<Integer> first_time(true);
+      if (first_time.load(std::memory_order_relaxed)) {
+        std::lock_guard<std::mutex> guard(mutex);
+        if (first_time.load(std::memory_order_relaxed)) {
+          { // Set GradOp
+            auto nodes = Basis<ValueType,1,ORDER>::Nodes();
+            SCTL_ASSERT(nodes.Dim(1) == ORDER);
+            Matrix<ValueType> M(ORDER, ORDER);
+            for (Integer i = 0; i < ORDER; i++) { // Set M
+              Real x = nodes[0][i];
+              for (Integer j = 0; j < ORDER; j++) {
+                M[j][i] = 0;
+                for (Integer l = 0; l < ORDER; l++) {
+                  if (l != j) {
+                    Real M_ = 1;
+                    for (Integer k = 0; k < ORDER; k++) {
+                      if (k != j && k != l) M_ *= (x - nodes[0][k]);
+                      if (k != j) M_ /= (nodes[0][j] - nodes[0][k]);
+                    }
+                    M[j][i] += M_;
+                  }
+                }
+              }
+            }
+            for (Integer d = 0; d < DIM; d++) {
+              GradOp[d].ReInit(Size(), Size());
+              GradOp[d] = 0;
+              Integer stride0 = sctl::pow<Integer>(ORDER, d);
+              Integer repeat0 = sctl::pow<Integer>(ORDER, d);
+              Integer stride1 = sctl::pow<Integer>(ORDER, d+1);
+              Integer repeat1 = sctl::pow<Integer>(ORDER, DIM-d-1);
+              for (Integer k1 = 0; k1 < repeat1; k1++) {
+                for (Integer i = 0; i < ORDER; i++) {
+                  for (Integer j = 0; j < ORDER; j++) {
+                    for (Integer k0 = 0; k0 < repeat0; k0++) {
+                      GradOp[d][k1*stride1 + i*stride0 + k0][k1*stride1 + j*stride0 + k0] = M[i][j];
+                    }
+                  }
+                }
+              }
+            }
+          }
+          std::atomic_thread_fence(std::memory_order_seq_cst);
+          first_time.store(false);
+        }
+      }
+
+      if (dX.Dim() != X.Dim()*DIM) dX.ReInit(X.Dim()*DIM);
+      for (Long i = 0; i < X.Dim(); i++) {
+        const Matrix<ValueType> Vi(1, Size(), (Iterator<ValueType>)(ConstIterator<ValueType>)X[i].NodeValues_, false);
+        for (Integer k = 0; k < DIM; k++) {
+          Matrix<ValueType> Vo(1, Size(), dX[i*DIM+k].NodeValues_, false);
+          Matrix<ValueType>::GEMM(Vo, Vi, GradOp[k]);
+        }
+      }
+    }
+    static EvalOpType SetupEval(const Matrix<ValueType>& X) {
+      Long N = X.Dim(1);
+      SCTL_ASSERT(X.Dim(0) == DIM);
+      Matrix<ValueType> M(Size(), N);
+      { // Set M
+        auto nodes = Basis<ValueType,1,ORDER>::Nodes();
+        Integer NN = nodes.Dim(1);
+        Matrix<ValueType> M_(NN, DIM*N);
+        for (Long i = 0; i < DIM*N; i++) {
+          ValueType x = X[0][i];
+          for (Integer j = 0; j < NN; j++) {
+            ValueType y = 1;
+            for (Integer k = 0; k < NN; k++) {
+              y *= (j==k ? 1 : (nodes[0][k] - x) / (nodes[0][k] - nodes[0][j]));
+            }
+            M_[j][i] = y;
+          }
+        }
+        if (DIM == 1) {
+          SCTL_ASSERT(M.Dim(0) == M_.Dim(0));
+          SCTL_ASSERT(M.Dim(1) == M_.Dim(1));
+          M = M_;
+        } else {
+          Integer NNN = 1;
+          M = 1;
+          for (Integer d = 0; d < DIM; d++) {
+            for (Integer k = 1; k < NN; k++) {
+              for (Integer j = 0; j < NNN; j++) {
+                for (Long i = 0; i < N; i++) {
+                  M[k*NNN+j][i] = M[j][i] * M_[k][d*N+i];
+                }
+              }
+            }
+            { // k = 0
+              for (Integer j = 0; j < NNN; j++) {
+                for (Long i = 0; i < N; i++) {
+                  M[j][i] *= M_[0][d*N+i];
+                }
+              }
+            }
+            NNN *= NN;
+          }
+        }
+      }
+      return M;
+    }
+    static void Eval(Matrix<ValueType>& Y, const Vector<Basis>& X, const EvalOpType& M) {
+      Long N0 = X.Dim();
+      Long N1 = M.Dim(1);
+      SCTL_ASSERT(M.Dim(0) == Size());
+      if (Y.Dim(0) != N0 || Y.Dim(1) != N1) Y.ReInit(N0, N1);
+      for (Long i = 0; i < N0; i++) {
+        const Matrix<ValueType> X_(1,Size(),(Iterator<ValueType>)(ConstIterator<ValueType>)X[i].NodeValues_,false);
+        Matrix<ValueType> Y_(1,N1,Y[i],false);
+        Matrix<ValueType>::GEMM(Y_,X_,M);
+      }
+    }
+
+    const ValueType& operator[](Long i) const {
+      SCTL_ASSERT(i < Size());
+      return NodeValues_[i];
+    }
+    ValueType& operator[](Long i) {
+      SCTL_ASSERT(i < Size());
+      return NodeValues_[i];
+    }
+
+  private:
+    StaticArray<ValueType,Size()> NodeValues_;
+};
+
+template <Integer COORD_DIM, class Basis> class ElemList {
+  public:
+
+    using CoordBasis = Basis;
+    using CoordType = typename CoordBasis::ValueType;
+
+    static constexpr Integer CoordDim() {
+      return COORD_DIM;
+    }
+    static constexpr Integer ElemDim() {
+      return CoordBasis::Dim();
+    }
+
+    ElemList(Long Nelem = 0) {
+      ReInit(Nelem);
+    }
+    void ReInit(Long Nelem = 0) {
+      Nelem_ = Nelem;
+      X_.ReInit(Nelem_ * COORD_DIM);
+    }
+    void ReInit(const Vector<CoordBasis>& X) {
+      Nelem_ = X.Dim() / COORD_DIM;
+      SCTL_ASSERT(X.Dim() == Nelem_ * COORD_DIM);
+      X_ = X;
+    }
+    Long NElem() const {
+      return Nelem_;
+    }
+
+    CoordBasis& operator()(Long elem, Integer dim) {
+      SCTL_ASSERT(elem >= 0 && elem < Nelem_);
+      SCTL_ASSERT(dim >= 0 && dim < COORD_DIM);
+      return X_[elem*COORD_DIM+dim];
+    }
+    const CoordBasis& operator()(Long elem, Integer dim) const {
+      SCTL_ASSERT(elem >= 0 && elem < Nelem_);
+      SCTL_ASSERT(dim >= 0 && dim < COORD_DIM);
+      return X_[elem*COORD_DIM+dim];
+    }
+    const Vector<CoordBasis>& ElemVector() const {
+      return X_;
+    }
+
+  private:
+    static_assert(CoordBasis::Dim() <= CoordDim(), "Basis dimension can not be greater than COORD_DIM.");
+    Vector<CoordBasis> X_;
+    Long Nelem_;
+
+    mutable Vector<CoordBasis> dX_;
+};
+
+template <class Real> class Quadrature {
+
+    static Real machine_epsilon() {
+      Real eps=1;
+      while(eps*(Real)0.5+(Real)1.0>1.0) eps*=0.5;
+      return eps;
+    }
+
+    template <Integer DIM> static void DuffyQuad(Matrix<Real>& nodes, Vector<Real>& weights, const Vector<Real>& coord, Integer order, Real adapt = -1.0) {
+      SCTL_ASSERT(coord.Dim() == DIM);
+      static Real eps = machine_epsilon()*16;
+
+      Matrix<Real> qx;
+      Vector<Real> qw;
+      { // Set qx, qw
+        Vector<Real> qx0, qw0;
+        ChebBasis<Real>::quad_rule(order, qx0, qw0);
+
+        Integer N = sctl::pow<DIM,Integer>(order);
+        qx.ReInit(DIM,N);
+        qw.ReInit(N);
+        qw[0] = 1;
+
+        Integer N_ = 1;
+        for (Integer d = 0; d < DIM; d++) {
+          for (Integer j = 0; j < order; j++) {
+            for (Integer i = 0; i < N_; i++) {
+              for (Integer k = 0; k < d; k++) {
+                qx[k][j*N_+i] = qx[k][i];
+              }
+              qx[d][j*N_+i] = qx0[j];
+              qw[j*N_+i] = qw[i];
+            }
+          }
+          for (Integer j = 0; j < order; j++) {
+            for (Integer i = 0; i < N_; i++) {
+              qw[j*N_+i] *= qw0[j];
+            }
+          }
+          N_ *= order;
+        }
+      }
+
+      Vector<Real> X;
+      { // Set X
+        StaticArray<Real,2*DIM+2> X_;
+        X_[0] = 0;
+        X_[1] = adapt;
+        for (Integer i = 0; i < DIM; i++) {
+          X_[2*i+2] = sctl::fabs<Real>(coord[i]);
+          X_[2*i+3] = sctl::fabs<Real>(coord[i]-1);
+        }
+        std::sort((Iterator<Real>)X_, (Iterator<Real>)X_+2*DIM+2);
+
+        X.PushBack(std::max<Real>(0, X_[2*DIM]-1));
+        for (Integer i = 0; i < 2*DIM+2; i++) {
+          if (X[X.Dim()-1] < X_[i]) {
+            if (X.Dim())
+            X.PushBack(X_[i]);
+          }
+        }
+
+        /////////////////////////////////////////////////////////////////////////////////////////////////
+        Vector<Real> r(1);
+        r[0] = X[0];
+        for (Integer i = 1; i < X.Dim(); i++) {
+          while (r[r.Dim() - 1] > 0.0 && (order*0.5) * r[r.Dim() - 1] < X[i]) r.PushBack((order*0.5) * r[r.Dim() - 1]); // TODO
+          r.PushBack(X[i]);
+        }
+        X = r;
+        /////////////////////////////////////////////////////////////////////////////////////////////////
+      }
+
+      Vector<Real> nds, wts;
+      for (Integer k = 0; k < X.Dim()-1; k++) {
+        for (Integer dd = 0; dd < 2*DIM; dd++) {
+          Integer d0 = (dd>>1);
+          StaticArray<Real,2*DIM> range0, range1;
+          { // Set range0, range1
+            Integer d1 = (dd%2?1:-1);
+            for (Integer d = 0; d < DIM; d++) {
+              range0[d*2+0] = std::max<Real>(0,std::min<Real>(1,coord[d] - X[k]  ));
+              range0[d*2+1] = std::max<Real>(0,std::min<Real>(1,coord[d] + X[k]  ));
+              range1[d*2+0] = std::max<Real>(0,std::min<Real>(1,coord[d] - X[k+1]));
+              range1[d*2+1] = std::max<Real>(0,std::min<Real>(1,coord[d] + X[k+1]));
+            }
+            range0[d0*2+0] = std::max<Real>(0,std::min<Real>(1,coord[d0] + d1*X[k+0]));
+            range0[d0*2+1] = std::max<Real>(0,std::min<Real>(1,coord[d0] + d1*X[k+0]));
+            range1[d0*2+0] = std::max<Real>(0,std::min<Real>(1,coord[d0] + d1*X[k+1]));
+            range1[d0*2+1] = std::max<Real>(0,std::min<Real>(1,coord[d0] + d1*X[k+1]));
+          }
+          { // if volume(range0, range1) == 0 then continue
+            Real v0 = 1, v1 = 1;
+            for (Integer d = 0; d < DIM; d++) {
+              if (d == d0) {
+                v0 *= sctl::fabs<Real>(range0[d*2+0]-range1[d*2+0]);
+                v1 *= sctl::fabs<Real>(range0[d*2+0]-range1[d*2+0]);
+              } else {
+                v0 *= range0[d*2+1]-range0[d*2+0];
+                v1 *= range1[d*2+1]-range1[d*2+0];
+              }
+            }
+            if (v0 < eps && v1 < eps) continue;
+          }
+          for (Integer i = 0; i < qx.Dim(1); i++) { // Set nds, wts
+            Real w = qw[i];
+            Real z = qx[d0][i];
+            for (Integer d = 0; d < DIM; d++) {
+              Real y = qx[d][i];
+              nds.PushBack((range0[d*2+0]*(1-y) + range0[d*2+1]*y)*(1-z) + (range1[d*2+0]*(1-y) + range1[d*2+1]*y)*z);
+              if (d == d0) {
+                w *= abs(range1[d*2+0] - range0[d*2+0]);
+              } else {
+                w *= (range0[d*2+1] - range0[d*2+0])*(1-z) + (range1[d*2+1] - range1[d*2+0])*z;
+              }
+            }
+            wts.PushBack(w);
+          }
+        }
+      }
+      nodes = Matrix<Real>(nds.Dim()/DIM,DIM,nds.begin()).Transpose();
+      weights = wts;
+    }
+
+    template <Integer DIM> static void TensorProductGaussQuad(Matrix<Real>& nodes, Vector<Real>& weights, Integer order) {
+      Vector<Real> coord(DIM);
+      coord = 0;
+      coord[0] = -10;
+      DuffyQuad<DIM>(nodes, weights, coord, order);
+    }
+
+
+
+    template <class DensityBasis, class ElemList, class Kernel> static void SetupSingular(Matrix<Real>& M_singular, const Matrix<Real>& trg_nds, const ElemList& elem_lst, const Kernel& kernel, Integer order_singular = 10, Integer order_direct = 10) {
+      using CoordBasis = typename ElemList::CoordBasis;
+      using CoordEvalOpType = typename CoordBasis::EvalOpType;
+      using DensityEvalOpType = typename DensityBasis::EvalOpType;
+
+      constexpr Integer CoordDim = ElemList::CoordDim();
+      constexpr Integer ElemDim = ElemList::ElemDim();
+      constexpr Integer KDIM0 = Kernel::SrcDim();
+      constexpr Integer KDIM1 = Kernel::TrgDim();
+
+      const Long Nelem = elem_lst.NElem();
+      const Integer Ntrg = trg_nds.Dim(1);
+      SCTL_ASSERT(trg_nds.Dim(0) == ElemDim);
+
+      Vector<Real> Xt;
+      { // Set Xt
+        auto Meval = CoordBasis::SetupEval(trg_nds);
+        eval_basis(Xt, elem_lst.ElemVector(), ElemList::CoordDim(), trg_nds.Dim(1), Meval);
+      }
+      SCTL_ASSERT(Xt.Dim() == Nelem * Ntrg * CoordDim);
+
+      const Vector<CoordBasis>& X = elem_lst.ElemVector();
+      Vector<CoordBasis> dX;
+      CoordBasis::Grad(dX, X);
+
+      auto& M = M_singular;
+      M.ReInit(Nelem * KDIM0 * DensityBasis::Size(), KDIM1 * Ntrg);
+      #pragma omp parallel for schedule(static)
+      for (Long i = 0; i < Ntrg; i++) { // Set M (singular)
+        Matrix<Real> quad_nds;
+        Vector<Real> quad_wts;
+        { // Set quad_nds, quad_wts
+          StaticArray<Real,ElemDim> trg_node_;
+          for (Integer k = 0; k < ElemDim; k++) {
+            trg_node_[k] = trg_nds[k][i];
+          }
+          Vector<Real> trg_node(ElemDim, trg_node_, false);
+          DuffyQuad<ElemDim>(quad_nds, quad_wts, trg_node, order_singular);
+        }
+        const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds);
+        Integer Nnds = quad_wts.Dim();
+
+        Vector<Real> X_, dX_, Xa_, Xn_;
+        { // Set X_, dX_
+          eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp);
+          eval_basis(dX_, dX, CoordDim * ElemDim, Nnds, CoordEvalOp);
+        }
+        if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_
+          Long N = Nelem*Nnds;
+          Xa_.ReInit(N);
+          Xn_.ReInit(N*CoordDim);
+          for (Long j = 0; j < N; j++) {
+            StaticArray<Real,CoordDim> normal;
+            normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3];
+            normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5];
+            normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1];
+            Xa_[j] = sctl::sqrt<Real>(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+            Real invXa = 1/Xa_[j];
+            Xn_[j*3+0] = normal[0] * invXa;
+            Xn_[j*3+1] = normal[1] * invXa;
+            Xn_[j*3+2] = normal[2] * invXa;
+          }
+        }
+
+        DensityEvalOpType DensityEvalOp;
+        if (std::is_same<CoordBasis,DensityBasis>::value) {
+          DensityEvalOp = CoordEvalOp;
+        } else {
+          DensityEvalOp = DensityBasis::SetupEval(quad_nds);
+        }
+
+        for (Long j = 0; j < Nelem; j++) {
+          Matrix<Real> M__(Nnds * KDIM0, KDIM1);
+          { // Set kernel matrix M__
+            const Vector<Real> X0_(CoordDim, (Iterator<Real>)Xt.begin() + (j * Ntrg + i) * CoordDim, false);
+            const Vector<Real> X__(Nnds * CoordDim, X_.begin() + j * Nnds * CoordDim, false);
+            const Vector<Real> Xn__(Nnds * CoordDim, Xn_.begin() + j * Nnds * CoordDim, false);
+            kernel.template KernelMatrix<Real>(M__, X0_, X__, Xn__);
+          }
+          for (Long k0 = 0; k0 < KDIM0; k0++) {
+            for (Long k1 = 0; k1 < KDIM1; k1++) {
+              for (Long l = 0; l < DensityBasis::Size(); l++) {
+                Real M_lk = 0;
+                for (Long n = 0; n < Nnds; n++) {
+                  Real quad_wt = Xa_[j * Nnds + n] * quad_wts[n];
+                  M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1];
+                }
+                M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1 * Ntrg + i] = M_lk;
+              }
+            }
+          }
+        }
+      }
+      { // Set M (subtract direct)
+        Matrix<Real> quad_nds;
+        Vector<Real> quad_wts;
+        TensorProductGaussQuad<ElemDim>(quad_nds, quad_wts, order_direct);
+        const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds);
+        Integer Nnds = quad_wts.Dim();
+
+        Vector<Real> X_, dX_, Xa_, Xn_;
+        { // Set X_, dX_
+          eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp);
+          eval_basis(dX_, dX, CoordDim * ElemDim, Nnds, CoordEvalOp);
+        }
+        if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_
+          Long N = Nelem*Nnds;
+          Xa_.ReInit(N);
+          Xn_.ReInit(N*CoordDim);
+          for (Long j = 0; j < N; j++) {
+            StaticArray<Real,CoordDim> normal;
+            normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3];
+            normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5];
+            normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1];
+            Xa_[j] = sctl::sqrt<Real>(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+            Real invXa = 1/Xa_[j];
+            Xn_[j*3+0] = normal[0] * invXa;
+            Xn_[j*3+1] = normal[1] * invXa;
+            Xn_[j*3+2] = normal[2] * invXa;
+          }
+        }
+
+        DensityEvalOpType DensityEvalOp;
+        if (std::is_same<CoordBasis,DensityBasis>::value) {
+          DensityEvalOp = CoordEvalOp;
+        } else {
+          DensityEvalOp = DensityBasis::SetupEval(quad_nds);
+        }
+
+        #pragma omp parallel for schedule(static)
+        for (Long i = 0; i < Ntrg; i++) { // Subtract direct contribution
+          for (Long j = 0; j < Nelem; j++) {
+            Matrix<Real> M__(Nnds * KDIM0, KDIM1);
+            { // Set kernel matrix M__
+              const Vector<Real> X0_(CoordDim, (Iterator<Real>)Xt.begin() + (j * Ntrg + i) * CoordDim, false);
+              const Vector<Real> X__(Nnds * CoordDim, X_.begin() + j * Nnds * CoordDim, false);
+              const Vector<Real> Xn__(Nnds * CoordDim, Xn_.begin() + j * Nnds * CoordDim, false);
+              kernel.template KernelMatrix<Real>(M__, X0_, X__, Xn__);
+            }
+            for (Long k0 = 0; k0 < KDIM0; k0++) {
+              for (Long k1 = 0; k1 < KDIM1; k1++) {
+                for (Long l = 0; l < DensityBasis::Size(); l++) {
+                  Real M_lk = 0;
+                  for (Long n = 0; n < Nnds; n++) {
+                    Real quad_wt = Xa_[j * Nnds + n] * quad_wts[n];
+                    M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1];
+                  }
+                  M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1 * Ntrg + i] -= M_lk;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+
+    template <class DensityBasis> static void EvalSingular(Matrix<Real>& U, const Vector<DensityBasis>& density, const Matrix<Real>& M, Integer KDIM0_, Integer KDIM1_) {
+      if (M.Dim(0) == 0 || M.Dim(1) == 0) {
+        U.ReInit(0,0);
+        return;
+      }
+
+      const Long Ntrg = M.Dim(1) / KDIM1_;
+      SCTL_ASSERT(M.Dim(1) == KDIM1_ * Ntrg);
+
+      const Long Nelem = M.Dim(0) / (KDIM0_ * DensityBasis::Size());
+      SCTL_ASSERT(M.Dim(0) == Nelem * KDIM0_ * DensityBasis::Size());
+
+      const Integer dof = density.Dim() / (Nelem * KDIM0_);
+      SCTL_ASSERT(density.Dim() == Nelem * dof * KDIM0_);
+
+      if (U.Dim(0) != Nelem * dof * KDIM1_ || U.Dim(1) != Ntrg) {
+        U.ReInit(Nelem * dof * KDIM1_, Ntrg);
+        U = 0;
+      }
+      for (Long j = 0; j < Nelem; j++) {
+        const Matrix<Real> M_(KDIM0_ * DensityBasis::Size(), KDIM1_ * Ntrg, (Iterator<Real>)M[j * KDIM0_ * DensityBasis::Size()], false);
+        Matrix<Real> U_(dof, KDIM1_ * Ntrg, U[j*dof*KDIM1_], false);
+        Matrix<Real> F_(dof, KDIM0_ * DensityBasis::Size());
+        for (Long i = 0; i < dof; i++) {
+          for (Long k = 0; k < KDIM0_; k++) {
+            for (Long l = 0; l < DensityBasis::Size(); l++) {
+              F_[i][k * DensityBasis::Size() + l] = density[(j * dof + i) * KDIM0_ + k][l];
+            }
+          }
+        }
+        Matrix<Real>::GEMM(U_, F_, M_);
+      }
+    }
+
+
+
+    template <Integer DIM> struct PointData {
+      bool operator<(const PointData& p) const {
+        return mid < p.mid;
+      }
+
+      Long rank;
+      Long surf_rank;
+      Morton<DIM> mid;
+      StaticArray<Real,DIM> coord;
+      Real radius2;
+    };
+
+    template <class T1, class T2> struct Pair {
+      Pair() {}
+
+      Pair(T1 x, T2 y) : first(x), second(y) {}
+
+      bool operator<(const Pair& p) const {
+        return (first < p.first) || (((first == p.first) && (second < p.second)));
+      }
+
+      T1 first;
+      T2 second;
+    };
+
+    template <class ElemList> static void BuildNbrList(Vector<Pair<Long,Long>>& pair_lst, const Vector<Real>& Xt, const Vector<Long>& trg_surf, const ElemList& elem_lst, Real distance_factor, Real period_length, const Comm& comm) {
+      using CoordBasis = typename ElemList::CoordBasis;
+      constexpr Integer CoordDim = ElemList::CoordDim();
+      constexpr Integer ElemDim = ElemList::ElemDim();
+      using PtData = PointData<CoordDim>;
+      const Integer rank = comm.Rank();
+
+      Real R0 = 0;
+      StaticArray<Real,CoordDim> X0;
+      { // Find bounding box
+        Long N = Xt.Dim() / CoordDim;
+        SCTL_ASSERT(Xt.Dim() == N * CoordDim);
+        SCTL_ASSERT(N);
+
+        StaticArray<Real,CoordDim*2> Xloc;
+        StaticArray<Real,CoordDim*2> Xglb;
+        for (Integer k = 0; k < CoordDim; k++) {
+          Xloc[0*CoordDim+k] = Xt[k];
+          Xloc[1*CoordDim+k] = Xt[k];
+        }
+        for (Long i = 0; i < N; i++) {
+          for (Integer k = 0; k < CoordDim; k++) {
+            Xloc[0*CoordDim+k] = std::min<Real>(Xloc[0*CoordDim+k], Xt[i*CoordDim+k]);
+            Xloc[1*CoordDim+k] = std::max<Real>(Xloc[1*CoordDim+k], Xt[i*CoordDim+k]);
+          }
+        }
+        comm.Allreduce((ConstIterator<Real>)Xloc+0*CoordDim, (Iterator<Real>)Xglb+0*CoordDim, CoordDim, Comm::CommOp::MIN);
+        comm.Allreduce((ConstIterator<Real>)Xloc+1*CoordDim, (Iterator<Real>)Xglb+1*CoordDim, CoordDim, Comm::CommOp::MAX);
+        for (Integer k = 0; k < CoordDim; k++) {
+          R0 = std::max(R0, Xglb[1*CoordDim+k]-Xglb[0*CoordDim+k]);
+        }
+
+        R0 = R0 * 2.0;
+        for (Integer k = 0; k < CoordDim; k++) {
+          X0[k] = Xglb[k] - R0*0.25;
+        }
+      }
+      if (period_length > 0) {
+        R0 = period_length;
+      }
+
+      Vector<PtData> PtSrc, PtTrg;
+      Integer order_upsample = (Integer)(const_pi<Real>() / distance_factor + 0.5);
+      { // Set PtSrc
+        const Vector<CoordBasis>& X_elem_lst = elem_lst.ElemVector();
+        Vector<CoordBasis> dX_elem_lst;
+        CoordBasis::Grad(dX_elem_lst, X_elem_lst);
+
+        Matrix<Real> nds;
+        Vector<Real> wts;
+        TensorProductGaussQuad<ElemDim>(nds, wts, order_upsample);
+        const Long Nnds = nds.Dim(1);
+
+        Vector<Real> X, dX;
+        const auto CoordEvalOp = CoordBasis::SetupEval(nds);
+        eval_basis(X, X_elem_lst, CoordDim, Nnds, CoordEvalOp);
+        eval_basis(dX, dX_elem_lst, CoordDim * ElemDim, Nnds, CoordEvalOp);
+
+        const Long N = X.Dim() / CoordDim;
+        const Long Nelem = elem_lst.NElem();
+        SCTL_ASSERT(X.Dim() == N * CoordDim);
+        SCTL_ASSERT(N == Nelem * Nnds);
+
+        Long rank_offset, surf_rank_offset;
+        { // Set rank_offset, surf_rank_offset
+          comm.Scan(Ptr2ConstItr<Long>(&N,1), Ptr2Itr<Long>(&rank_offset,1), 1, Comm::CommOp::SUM);
+          comm.Scan(Ptr2ConstItr<Long>(&Nelem,1), Ptr2Itr<Long>(&surf_rank_offset,1), 1, Comm::CommOp::SUM);
+          surf_rank_offset -= Nelem;
+          rank_offset -= N;
+        }
+
+        PtSrc.ReInit(N);
+        const Real R0inv = 1.0 / R0;
+        for (Long i = 0; i < N; i++) { // Set coord
+          for (Integer k = 0; k < CoordDim; k++) {
+            PtSrc[i].coord[k] = (X[i*CoordDim+k] - X0[k]) * R0inv;
+          }
+        }
+        if (period_length > 0) { // Wrap-around coord
+          for (Long i = 0; i < N; i++) {
+            auto& x = PtSrc[i].coord;
+            for (Integer k = 0; k < CoordDim; k++) {
+              x[k] -= (Long)(x[k]);
+            }
+          }
+        }
+        for (Long i = 0; i < N; i++) { // Set radius2, mid, rank
+          Integer depth = 0;
+          { // Set radius2, depth
+            Real radius2 = 0;
+            for (Integer k0 = 0; k0 < ElemDim; k0++) {
+              Real R2 = 0;
+              for (Integer k1 = 0; k1 < CoordDim; k1++) {
+                Real dX_ = dX[(i*CoordDim+k1)*ElemDim+k0];
+                R2 += dX_*dX_;
+              }
+              radius2 = std::max(radius2, R2);
+            }
+            radius2 *= R0inv*R0inv * distance_factor*distance_factor;
+            PtSrc[i].radius2 = radius2;
+
+            Long Rinv = (Long)(1.0/radius2);
+            while (Rinv > 0) {
+              Rinv = (Rinv>>2);
+              depth++;
+            }
+          }
+          PtSrc[i].mid = Morton<CoordDim>((Iterator<Real>)PtSrc[i].coord, std::min(Morton<CoordDim>::MaxDepth(),depth));
+          PtSrc[i].rank = rank_offset + i;
+        }
+        for (Long i = 0 ; i < Nelem; i++) { // Set surf_rank
+          for (Long j = 0; j < Nnds; j++) {
+            PtSrc[i*Nnds+j].surf_rank = surf_rank_offset + i;
+          }
+        }
+
+        Vector<PtData> PtSrcSorted;
+        comm.HyperQuickSort(PtSrc, PtSrcSorted);
+        PtSrc.Swap(PtSrcSorted);
+      }
+      { // Set PtTrg
+        const Long N = Xt.Dim() / CoordDim;
+        SCTL_ASSERT(Xt.Dim() == N * CoordDim);
+
+        Long rank_offset;
+        { // Set rank_offset
+          comm.Scan(Ptr2ConstItr<Long>(&N,1), Ptr2Itr<Long>(&rank_offset,1), 1, Comm::CommOp::SUM);
+          rank_offset -= N;
+        }
+
+        PtTrg.ReInit(N);
+        const Real R0inv = 1.0 / R0;
+        for (Long i = 0; i < N; i++) { // Set coord
+          for (Integer k = 0; k < CoordDim; k++) {
+            PtTrg[i].coord[k] = (Xt[i*CoordDim+k] - X0[k]) * R0inv;
+          }
+        }
+        if (period_length > 0) { // Wrap-around coord
+          for (Long i = 0; i < N; i++) {
+            auto& x = PtTrg[i].coord;
+            for (Integer k = 0; k < CoordDim; k++) {
+              x[k] -= (Long)(x[k]);
+            }
+          }
+        }
+        for (Long i = 0; i < N; i++) { // Set radius2, mid, rank
+          PtTrg[i].radius2 = 0;
+          PtTrg[i].mid = Morton<CoordDim>((Iterator<Real>)PtTrg[i].coord);
+          PtTrg[i].rank = rank_offset + i;
+        }
+        if (trg_surf.Dim()) { // Set surf_rank
+          SCTL_ASSERT(trg_surf.Dim() == N);
+          for (Long i = 0; i < N; i++) {
+            PtTrg[i].surf_rank = trg_surf[i];
+          }
+        } else {
+          for (Long i = 0; i < N; i++) {
+            PtTrg[i].surf_rank = -1;
+          }
+        }
+
+        Vector<PtData> PtTrgSorted;
+        comm.HyperQuickSort(PtTrg, PtTrgSorted);
+        PtTrg.Swap(PtTrgSorted);
+      }
+
+      Tree<CoordDim> tree(comm);
+      { // Init tree
+        Vector<Real> Xall(PtSrc.Dim()+PtTrg.Dim());
+        { // Set Xall
+          Xall.ReInit((PtSrc.Dim()+PtTrg.Dim())*CoordDim);
+          Long Nsrc = PtSrc.Dim();
+          Long Ntrg = PtTrg.Dim();
+          for (Long i = 0; i < Nsrc; i++) {
+            for (Integer k = 0; k < CoordDim; k++) {
+              Xall[i*CoordDim+k] = PtSrc[i].coord[k];
+            }
+          }
+          for (Long i = 0; i < Ntrg; i++) {
+            for (Integer k = 0; k < CoordDim; k++) {
+              Xall[(Nsrc+i)*CoordDim+k] = PtTrg[i].coord[k];
+            }
+          }
+        }
+        tree.UpdateRefinement(Xall, 1000, true, period_length>0);
+      }
+      { // Repartition PtSrc, PtTrg
+        PtData splitter;
+        splitter.mid = tree.GetPartitionMID()[rank];
+        comm.PartitionS(PtSrc, splitter);
+        comm.PartitionS(PtTrg, splitter);
+      }
+      { // Add tree data PtSrc
+        const auto& node_mid = tree.GetNodeMID();
+        const Long N = node_mid.Dim();
+        SCTL_ASSERT(N);
+
+        Vector<Long> dsp(N), cnt(N);
+        for (Long i = 0; i < N; i++) {
+          PtData m0;
+          m0.mid = node_mid[i];
+          dsp[i] = std::lower_bound(PtSrc.begin(), PtSrc.end(), m0) - PtSrc.begin();
+        }
+        for (Long i = 0; i < N-1; i++) {
+          cnt[i] = dsp[i+1] - dsp[i];
+        }
+        cnt[N-1] = PtSrc.Dim() - dsp[N-1];
+        tree.AddData("PtSrc", PtSrc, cnt);
+      }
+      tree.template Broadcast<PtData>("PtSrc");
+
+      { // Build pair_lst
+        Vector<Long> cnt;
+        Vector<PtData> PtSrc;
+        tree.GetData(PtSrc, cnt, "PtSrc");
+        const auto& node_mid = tree.GetNodeMID();
+        const auto& node_attr = tree.GetNodeAttr();
+
+        Vector<Morton<CoordDim>> nbr_mid_tmp;
+        for (Long i = 0; i < node_mid.Dim(); i++) {
+          if (node_attr[i].Leaf && !node_attr[i].Ghost) {
+            Vector<Morton<CoordDim>> child_mid;
+            node_mid[i].Children(child_mid);
+            for (const auto& trg_mid : child_mid) {
+              Integer d0 = trg_mid.Depth();
+              Vector<PtData> Src, Trg;
+              { // Set Trg
+                PtData m0, m1;
+                m0.mid = trg_mid;
+                m1.mid = trg_mid.Next();
+                Long a = std::lower_bound(PtTrg.begin(), PtTrg.end(), m0) - PtTrg.begin();
+                Long b = std::lower_bound(PtTrg.begin(), PtTrg.end(), m1) - PtTrg.begin();
+                Trg.ReInit(b-a, PtTrg.begin()+a, false);
+                if (!Trg.Dim()) continue;
+              }
+
+              Vector<std::set<Long>> near_elem(Trg.Dim());
+              for (Integer d = 0; d <= d0; d++) {
+                trg_mid.NbrList(nbr_mid_tmp, d, period_length>0);
+                for (const auto& src_mid : nbr_mid_tmp) { // Set Src
+                  PtData m0, m1;
+                  m0.mid = src_mid;
+                  m1.mid = (d==d0 ? src_mid.Next() : src_mid.Ancestor(d+1));
+                  Long a = std::lower_bound(PtSrc.begin(), PtSrc.end(), m0) - PtSrc.begin();
+                  Long b = std::lower_bound(PtSrc.begin(), PtSrc.end(), m1) - PtSrc.begin();
+                  Src.ReInit(b-a, PtSrc.begin()+a, false);
+                  if (!Src.Dim()) continue;
+
+                  for (Long t = 0; t < Trg.Dim(); t++) { // set near_elem[t] <-- {s : dist(s,t) < radius(s)}
+                    for (Long s = 0; s < Src.Dim(); s++) {
+                      if (Trg[t].surf_rank != Src[s].surf_rank) {
+                        Real R2 = 0;
+                        for (Integer k = 0; k < CoordDim; k++) {
+                          Real dx = (Src[s].coord[k] - Trg[t].coord[k]);
+                          R2 += dx * dx;
+                        }
+                        if (R2 < Src[s].radius2) {
+                          near_elem[t].insert(Src[s].surf_rank);
+                        }
+                      }
+                    }
+                  }
+                }
+              }
+
+              for (Long t = 0; t < Trg.Dim(); t++) { // Set pair_lst
+                for (Long elem_idx : near_elem[t]) {
+                  pair_lst.PushBack(Pair<Long,Long>(elem_idx,Trg[t].rank));
+                }
+              }
+            }
+          }
+        }
+      }
+      { // Sort and repartition pair_lst
+        Vector<Pair<Long,Long>> pair_lst_sorted;
+        comm.HyperQuickSort(pair_lst, pair_lst_sorted);
+
+        Long surf_rank_offset;
+        const Long Nelem = elem_lst.NElem();
+        comm.Scan(Ptr2ConstItr<Long>(&Nelem,1), Ptr2Itr<Long>(&surf_rank_offset,1), 1, Comm::CommOp::SUM);
+        surf_rank_offset -= Nelem;
+
+        comm.PartitionS(pair_lst_sorted, Pair<Long,Long>(surf_rank_offset,0));
+        pair_lst.Swap(pair_lst_sorted);
+      }
+    }
+
+    template <class ElemList> static void BuildNbrListDeprecated(Vector<Pair<Long,Long>>& pair_lst, const Vector<Real>& Xt, const ElemList& elem_lst, const Matrix<Real>& surf_nds, Real distance_factor) {
+      using CoordBasis = typename ElemList::CoordBasis;
+      constexpr Integer CoordDim = ElemList::CoordDim();
+      constexpr Integer ElemDim = ElemList::ElemDim();
+      const Long Nelem = elem_lst.NElem();
+
+      const Long Ntrg = Xt.Dim() / CoordDim;
+      SCTL_ASSERT(Xt.Dim() == Ntrg * CoordDim);
+
+      Long Nnds, Nsurf_nds;
+      Vector<Real> X_surf, X, dX;
+      Integer order_upsample = (Integer)(const_pi<Real>() / distance_factor + 0.5);
+      { // Set X, dX
+        const Vector<CoordBasis>& X_elem_lst = elem_lst.ElemVector();
+        Vector<CoordBasis> dX_elem_lst;
+        CoordBasis::Grad(dX_elem_lst, X_elem_lst);
+
+        Matrix<Real> nds_upsample;
+        Vector<Real> wts_upsample;
+        TensorProductGaussQuad<ElemDim>(nds_upsample, wts_upsample, order_upsample);
+
+        Nnds = nds_upsample.Dim(1);
+        const auto CoordEvalOp = CoordBasis::SetupEval(nds_upsample);
+        eval_basis(X, X_elem_lst, CoordDim, nds_upsample.Dim(1), CoordEvalOp);
+        eval_basis(dX, dX_elem_lst, CoordDim * ElemDim, nds_upsample.Dim(1), CoordEvalOp);
+
+        Nsurf_nds = surf_nds.Dim(1);
+        const auto CoordEvalOp_surf = CoordBasis::SetupEval(surf_nds);
+        eval_basis(X_surf, X_elem_lst, CoordDim, Nsurf_nds, CoordEvalOp_surf);
+      }
+
+      Real d2 = distance_factor * distance_factor;
+      for (Long i = 0; i < Nelem; i++) {
+        std::set<Long> near_pts;
+        std::set<Long> self_pts;
+        for (Long j = 0; j < Nnds; j++) {
+          Real R2_max = 0;
+          StaticArray<Real, CoordDim> X0;
+          for (Integer k = 0; k < CoordDim; k++) {
+            X0[k] = X[(i*Nnds+j)*CoordDim+k];
+          }
+          for (Integer k0 = 0; k0 < ElemDim; k0++) {
+            Real R2 = 0;
+            for (Integer k1 = 0; k1 < CoordDim; k1++) {
+              Real dX_ = dX[((i*Nnds+j)*CoordDim+k1)*ElemDim+k0];
+              R2 += dX_*dX_;
+            }
+            R2_max = std::max(R2_max, R2*d2);
+          }
+
+          for (Long k = 0; k < Ntrg; k++) {
+            Real R2 = 0;
+            for (Integer l = 0; l < CoordDim; l++) {
+              Real dX = Xt[k*CoordDim+l]- X0[l];
+              R2 += dX * dX;
+            }
+            if (R2 < R2_max) near_pts.insert(k);
+          }
+        }
+        for (Long j = 0; j < Nsurf_nds; j++) {
+          StaticArray<Real, CoordDim> X0;
+          for (Integer k = 0; k < CoordDim; k++) {
+            X0[k] = X_surf[(i*Nsurf_nds+j)*CoordDim+k];
+          }
+          for (Long k = 0; k < Ntrg; k++) {
+            Real R2 = 0;
+            for (Integer l = 0; l < CoordDim; l++) {
+              Real dX = Xt[k*CoordDim+l]- X0[l];
+              R2 += dX * dX;
+            }
+            if (R2 == 0) self_pts.insert(k);
+          }
+        }
+        for (Long trg_idx : self_pts) {
+          near_pts.erase(trg_idx);
+        }
+        for (Long trg_idx : near_pts) {
+          pair_lst.PushBack(Pair<Long,Long>(i,trg_idx));
+        }
+      }
+    }
+
+    template <class DensityBasis, class ElemList, class Kernel> static void SetupNearSingular(Matrix<Real>& M_near_singular, Vector<Pair<Long,Long>>& pair_lst, const Vector<Real>& Xt_, const Vector<Long>& trg_surf, const ElemList& elem_lst, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm) {
+      static_assert(std::is_same<Real,typename DensityBasis::ValueType>::value);
+      static_assert(std::is_same<Real,typename ElemList::CoordType>::value);
+      static_assert(DensityBasis::Dim() == ElemList::ElemDim());
+      using CoordBasis = typename ElemList::CoordBasis;
+      using CoordEvalOpType = typename CoordBasis::EvalOpType;
+      using DensityEvalOpType = typename DensityBasis::EvalOpType;
+
+      constexpr Integer CoordDim = ElemList::CoordDim();
+      constexpr Integer ElemDim = ElemList::ElemDim();
+      constexpr Integer KDIM0 = Kernel::SrcDim();
+      constexpr Integer KDIM1 = Kernel::TrgDim();
+      const Long Nelem = elem_lst.NElem();
+
+      BuildNbrList(pair_lst, Xt_, trg_surf, elem_lst, 2.5/order_direct, period_length, comm);
+      const Long Ninterac = pair_lst.Dim();
+
+      Vector<Real> Xt;
+      { // Set Xt
+        Integer rank = comm.Rank();
+        Integer np = comm.Size();
+
+        Vector<Long> splitter_ranks;
+        { // Set splitter_ranks
+          Vector<Long> cnt(np);
+          const Long N = Xt_.Dim() / CoordDim;
+          comm.Allgather(Ptr2ConstItr<Long>(&N,1), 1, cnt.begin(), 1);
+          scan(splitter_ranks, cnt);
+        }
+
+        Vector<Long> scatter_index, recv_index, recv_cnt(np), recv_dsp(np);
+        { // Set scatter_index, recv_index, recv_cnt, recv_dsp
+          { // Set scatter_index, recv_index
+            Vector<Pair<Long,Long>> scatter_pair(pair_lst.Dim());
+            for (Long i = 0; i < pair_lst.Dim(); i++) {
+              scatter_pair[i] = Pair<Long,Long>(pair_lst[i].second,i);
+            }
+            omp_par::merge_sort(scatter_pair.begin(), scatter_pair.end());
+
+            recv_index.ReInit(scatter_pair.Dim());
+            scatter_index.ReInit(scatter_pair.Dim());
+            for (Long i = 0; i < scatter_index.Dim(); i++) {
+              recv_index[i] = scatter_pair[i].first;
+              scatter_index[i] = scatter_pair[i].second;
+            }
+          }
+          for (Integer i = 0; i < np; i++) {
+            recv_dsp[i] = std::lower_bound(recv_index.begin(), recv_index.end(), splitter_ranks[i]) - recv_index.begin();
+          }
+          for (Integer i = 0; i < np-1; i++) {
+            recv_cnt[i] = recv_dsp[i+1] - recv_dsp[i];
+          }
+          recv_cnt[np-1] = recv_index.Dim() - recv_dsp[np-1];
+        }
+
+        Vector<Long> send_index, send_cnt(np), send_dsp(np);
+        { // Set send_index, send_cnt, send_dsp
+          comm.Alltoall(recv_cnt.begin(), 1, send_cnt.begin(), 1);
+          scan(send_dsp, send_cnt);
+          send_index.ReInit(send_cnt[np-1] + send_dsp[np-1]);
+          comm.Alltoallv(recv_index.begin(), recv_cnt.begin(), recv_dsp.begin(), send_index.begin(), send_cnt.begin(), send_dsp.begin());
+        }
+
+        Vector<Real> Xt_send(send_index.Dim() * CoordDim);
+        for (Long i = 0; i < send_index.Dim(); i++) { // Set Xt_send
+          Long idx = send_index[i] - splitter_ranks[rank];
+          for (Integer k = 0; k < CoordDim; k++) {
+            Xt_send[i*CoordDim+k] = Xt_[idx*CoordDim+k];
+          }
+        }
+
+        Vector<Real> Xt_recv(recv_index.Dim() * CoordDim);
+        { // Set Xt_recv
+          for (Long i = 0; i < np; i++) {
+            send_cnt[i] *= CoordDim;
+            send_dsp[i] *= CoordDim;
+            recv_cnt[i] *= CoordDim;
+            recv_dsp[i] *= CoordDim;
+          }
+          comm.Alltoallv(Xt_send.begin(), send_cnt.begin(), send_dsp.begin(), Xt_recv.begin(), recv_cnt.begin(), recv_dsp.begin());
+        }
+
+        Xt.ReInit(scatter_index.Dim() * CoordDim);
+        for (Long i = 0; i < scatter_index.Dim(); i++) { // Set Xt
+          Long idx = scatter_index[i];
+          for (Integer k = 0; k < CoordDim; k++) {
+            Xt[idx*CoordDim+k] = Xt_recv[i*CoordDim+k];
+          }
+        }
+      }
+
+      const Vector<CoordBasis>& X = elem_lst.ElemVector();
+      Vector<CoordBasis> dX;
+      CoordBasis::Grad(dX, X);
+
+      Long elem_rank_offset;
+      { // Set elem_rank_offset
+        comm.Scan(Ptr2ConstItr<Long>(&Nelem,1), Ptr2Itr<Long>(&elem_rank_offset,1), 1, Comm::CommOp::SUM);
+        elem_rank_offset -= Nelem;
+      }
+
+      auto& M = M_near_singular;
+      M.ReInit(Ninterac * KDIM0 * DensityBasis::Size(), KDIM1);
+      #pragma omp parallel for schedule(static)
+      for (Long j = 0; j < Ninterac; j++) { // Set M (near-singular)
+        const Long src_idx = pair_lst[j].first - elem_rank_offset;
+
+        Real adapt = -1.0;
+        Tensor<Real,true,ElemDim,1> u0;
+        { // Set u0 (project target point to the surface patch in parameter space)
+          ConstIterator<Real> Xt_ = Xt.begin() + j * CoordDim;
+          const auto& nodes = CoordBasis::Nodes();
+
+          Long min_idx = -1;
+          Real min_R2 = 1e10;
+          for (Long i = 0; i < CoordBasis::Size(); i++) {
+            Real R2 = 0;
+            for (Integer k = 0; k < CoordDim; k++) {
+              Real dX = X[src_idx * CoordDim + k][i] - Xt_[k];
+              R2 += dX * dX;
+            }
+            if (R2 < min_R2) {
+              min_R2 = R2;
+              min_idx = i;
+            }
+          }
+          SCTL_ASSERT(min_idx >= 0);
+          for (Integer k = 0; k < ElemDim; k++) {
+            u0(k,0) = nodes[k][min_idx];
+          }
+
+          for (Integer i = 0; i < 2; i++) { // iterate
+            Matrix<Real> X_, dX_;
+            for (Integer k = 0; k < ElemDim; k++) {
+              u0(k,0) = std::min(1.0, u0(k,0));
+              u0(k,0) = std::max(0.0, u0(k,0));
+            }
+            const auto eval_op = CoordBasis::SetupEval(Matrix<Real>(ElemDim,1,u0.begin(),false));
+            CoordBasis::Eval(X_, Vector<CoordBasis>(CoordDim,(Iterator<CoordBasis>)X.begin()+src_idx*CoordDim,false),eval_op);
+            CoordBasis::Eval(dX_, Vector<CoordBasis>(CoordDim*ElemDim,dX.begin()+src_idx*CoordDim*ElemDim,false),eval_op);
+
+            const Tensor<Real,false,CoordDim,1> x0((Iterator<Real>)Xt_);
+            const Tensor<Real,false,CoordDim,1> x(X_.begin());
+            const Tensor<Real,false,CoordDim,ElemDim> x_u(dX_.begin());
+            auto inv = [](const Tensor<Real,true,2,2>& M) {
+              Tensor<Real,true,2,2> Minv;
+              Real det_inv = 1.0 / (M(0,0)*M(1,1) - M(1,0)*M(0,1));
+              Minv(0,0) = M(1,1) * det_inv;
+              Minv(0,1) =-M(0,1) * det_inv;
+              Minv(1,0) =-M(1,0) * det_inv;
+              Minv(1,1) = M(0,0) * det_inv;
+              return Minv;
+            };
+            auto du = inv(x_u.RotateRight()*x_u) * x_u.RotateRight()*(x0-x);
+            u0 = u0 + du;
+
+            auto x_u_squared = x_u.RotateRight() * x_u;
+            adapt = sctl::sqrt<Real>( ((x0-x).RotateRight()*(x0-x))(0,0) / std::max<Real>(x_u_squared(0,0),x_u_squared(1,1)) );
+          }
+        }
+
+        Matrix<Real> quad_nds;
+        Vector<Real> quad_wts;
+        DuffyQuad<ElemDim>(quad_nds, quad_wts, Vector<Real>(ElemDim,u0.begin(),false), order_singular, adapt);
+        const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds);
+        Integer Nnds = quad_wts.Dim();
+
+        Vector<Real> X_, dX_, Xa_, Xn_;
+        { // Set X_, dX_
+          const Vector<CoordBasis> X__(CoordDim, (Iterator<CoordBasis>)X.begin() + src_idx * CoordDim, false);
+          const Vector<CoordBasis> dX__(CoordDim * ElemDim, (Iterator<CoordBasis>)dX.begin() + src_idx * CoordDim * ElemDim, false);
+          eval_basis(X_, X__, CoordDim, Nnds, CoordEvalOp);
+          eval_basis(dX_, dX__, CoordDim * ElemDim, Nnds, CoordEvalOp);
+        }
+        if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_
+          Xa_.ReInit(Nnds);
+          Xn_.ReInit(Nnds*CoordDim);
+          for (Long j = 0; j < Nnds; j++) {
+            StaticArray<Real,CoordDim> normal;
+            normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3];
+            normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5];
+            normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1];
+            Xa_[j] = sctl::sqrt<Real>(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+            Real invXa = 1/Xa_[j];
+            Xn_[j*3+0] = normal[0] * invXa;
+            Xn_[j*3+1] = normal[1] * invXa;
+            Xn_[j*3+2] = normal[2] * invXa;
+          }
+        }
+
+        DensityEvalOpType DensityEvalOp;
+        if (std::is_same<CoordBasis,DensityBasis>::value) {
+          DensityEvalOp = CoordEvalOp;
+        } else {
+          DensityEvalOp = DensityBasis::SetupEval(quad_nds);
+        }
+
+        Matrix<Real> M__(Nnds * KDIM0, KDIM1);
+        { // Set kernel matrix M__
+          const Vector<Real> X0_(CoordDim, (Iterator<Real>)Xt.begin() + j * CoordDim, false);
+          kernel.template KernelMatrix<Real>(M__, X0_, X_, Xn_);
+        }
+        for (Long k0 = 0; k0 < KDIM0; k0++) {
+          for (Long k1 = 0; k1 < KDIM1; k1++) {
+            for (Long l = 0; l < DensityBasis::Size(); l++) {
+              Real M_lk = 0;
+              for (Long n = 0; n < Nnds; n++) {
+                Real quad_wt = Xa_[n] * quad_wts[n];
+                M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1];
+              }
+              M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1] = M_lk;
+            }
+          }
+        }
+      }
+      { // Set M (subtract direct)
+        Matrix<Real> quad_nds;
+        Vector<Real> quad_wts;
+        TensorProductGaussQuad<ElemDim>(quad_nds, quad_wts, order_direct);
+        const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds);
+        Integer Nnds = quad_wts.Dim();
+
+        Vector<Real> X_, dX_, Xa_, Xn_;
+        { // Set X_, dX_
+          eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp);
+          eval_basis(dX_, dX, CoordDim * ElemDim, Nnds, CoordEvalOp);
+        }
+        if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_
+          Long N = Nelem*Nnds;
+          Xa_.ReInit(N);
+          Xn_.ReInit(N*CoordDim);
+          for (Long j = 0; j < N; j++) {
+            StaticArray<Real,CoordDim> normal;
+            normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3];
+            normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5];
+            normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1];
+            Xa_[j] = sctl::sqrt<Real>(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+            Real invXa = 1/Xa_[j];
+            Xn_[j*3+0] = normal[0] * invXa;
+            Xn_[j*3+1] = normal[1] * invXa;
+            Xn_[j*3+2] = normal[2] * invXa;
+          }
+        }
+
+        DensityEvalOpType DensityEvalOp;
+        if (std::is_same<CoordBasis,DensityBasis>::value) {
+          DensityEvalOp = CoordEvalOp;
+        } else {
+          DensityEvalOp = DensityBasis::SetupEval(quad_nds);
+        }
+
+        #pragma omp parallel for schedule(static)
+        for (Long j = 0; j < Ninterac; j++) { // Subtract direct contribution
+          const Long src_idx = pair_lst[j].first - elem_rank_offset;
+
+          Matrix<Real> M__(Nnds * KDIM0, KDIM1);
+          { // Set kernel matrix M__
+            const Vector<Real> X0_(CoordDim, (Iterator<Real>)Xt.begin() + j * CoordDim, false);
+            Vector<Real> X__(Nnds * CoordDim, X_.begin() + src_idx * Nnds * CoordDim, false);
+            Vector<Real> Xn__(Nnds * CoordDim, Xn_.begin() + src_idx * Nnds * CoordDim, false);
+            kernel.template KernelMatrix<Real>(M__, X0_, X__, Xn__);
+          }
+          for (Long k0 = 0; k0 < KDIM0; k0++) {
+            for (Long k1 = 0; k1 < KDIM1; k1++) {
+              for (Long l = 0; l < DensityBasis::Size(); l++) {
+                Real M_lk = 0;
+                for (Long n = 0; n < Nnds; n++) {
+                  Real quad_wt = Xa_[src_idx * Nnds + n] * quad_wts[n];
+                  M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1];
+                }
+                M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1] -= M_lk;
+              }
+            }
+          }
+        }
+      }
+    }
+
+    template <class DensityBasis> static void EvalNearSingular(Vector<Real>& U, const Vector<DensityBasis>& density, const Matrix<Real>& M, const Vector<Pair<Long,Long>>& pair_lst, Long Nelem_, Long Ntrg_, Integer KDIM0_, Integer KDIM1_, const Comm& comm) {
+      const Long Ninterac = pair_lst.Dim();
+      const Integer dof = density.Dim() / Nelem_ / KDIM0_;
+      SCTL_ASSERT(density.Dim() == Nelem_ * dof * KDIM0_);
+
+      Long elem_rank_offset;
+      { // Set elem_rank_offset
+        comm.Scan(Ptr2ConstItr<Long>(&Nelem_,1), Ptr2Itr<Long>(&elem_rank_offset,1), 1, Comm::CommOp::SUM);
+        elem_rank_offset -= Nelem_;
+      }
+
+      Vector<Real> U_loc(Ninterac*dof*KDIM1_);
+      for (Long j = 0; j < Ninterac; j++) {
+        const Long src_idx = pair_lst[j].first - elem_rank_offset;
+        const Matrix<Real> M_(KDIM0_ * DensityBasis::Size(), KDIM1_, (Iterator<Real>)M[j * KDIM0_ * DensityBasis::Size()], false);
+        Matrix<Real> U_(dof, KDIM1_, U_loc.begin() + j*dof*KDIM1_, false);
+        Matrix<Real> F_(dof, KDIM0_ * DensityBasis::Size());
+        for (Long i = 0; i < dof; i++) {
+          for (Long k = 0; k < KDIM0_; k++) {
+            for (Long l = 0; l < DensityBasis::Size(); l++) {
+              F_[i][k * DensityBasis::Size() + l] = density[(src_idx * dof + i) * KDIM0_ + k][l];
+            }
+          }
+        }
+        Matrix<Real>::GEMM(U_, F_, M_);
+      }
+
+      if (U.Dim() != Ntrg_ * dof * KDIM1_) {
+        U.ReInit(Ntrg_ * dof * KDIM1_);
+        U = 0;
+      }
+      { // Set U
+        Integer rank = comm.Rank();
+        Integer np = comm.Size();
+
+        Vector<Long> splitter_ranks;
+        { // Set splitter_ranks
+          Vector<Long> cnt(np);
+          comm.Allgather(Ptr2ConstItr<Long>(&Ntrg_,1), 1, cnt.begin(), 1);
+          scan(splitter_ranks, cnt);
+        }
+
+        Vector<Long> scatter_index, send_index, send_cnt(np), send_dsp(np);
+        { // Set scatter_index, send_index, send_cnt, send_dsp
+          { // Set scatter_index, send_index
+            Vector<Pair<Long,Long>> scatter_pair(pair_lst.Dim());
+            for (Long i = 0; i < pair_lst.Dim(); i++) {
+              scatter_pair[i] = Pair<Long,Long>(pair_lst[i].second,i);
+            }
+            omp_par::merge_sort(scatter_pair.begin(), scatter_pair.end());
+
+            send_index.ReInit(scatter_pair.Dim());
+            scatter_index.ReInit(scatter_pair.Dim());
+            for (Long i = 0; i < scatter_index.Dim(); i++) {
+              send_index[i] = scatter_pair[i].first;
+              scatter_index[i] = scatter_pair[i].second;
+            }
+          }
+          for (Integer i = 0; i < np; i++) {
+            send_dsp[i] = std::lower_bound(send_index.begin(), send_index.end(), splitter_ranks[i]) - send_index.begin();
+          }
+          for (Integer i = 0; i < np-1; i++) {
+            send_cnt[i] = send_dsp[i+1] - send_dsp[i];
+          }
+          send_cnt[np-1] = send_index.Dim() - send_dsp[np-1];
+        }
+
+        Vector<Long> recv_index, recv_cnt(np), recv_dsp(np);
+        { // Set recv_index, recv_cnt, recv_dsp
+          comm.Alltoall(send_cnt.begin(), 1, recv_cnt.begin(), 1);
+          scan(recv_dsp, recv_cnt);
+          recv_index.ReInit(recv_cnt[np-1] + recv_dsp[np-1]);
+          comm.Alltoallv(send_index.begin(), send_cnt.begin(), send_dsp.begin(), recv_index.begin(), recv_cnt.begin(), recv_dsp.begin());
+        }
+
+        Vector<Real> U_send(scatter_index.Dim() * dof * KDIM1_);
+        for (Long i = 0; i < scatter_index.Dim(); i++) {
+          Long idx = scatter_index[i]*dof*KDIM1_;
+          for (Long k = 0; k < dof * KDIM1_; k++) {
+            U_send[i*dof*KDIM1_ + k] = U_loc[idx + k];
+          }
+        }
+
+        Vector<Real> U_recv(recv_index.Dim() * dof * KDIM1_);
+        { // Set U_recv
+          for (Long i = 0; i < np; i++) {
+            send_cnt[i] *= dof * KDIM1_;
+            send_dsp[i] *= dof * KDIM1_;
+            recv_cnt[i] *= dof * KDIM1_;
+            recv_dsp[i] *= dof * KDIM1_;
+          }
+          comm.Alltoallv(U_send.begin(), send_cnt.begin(), send_dsp.begin(), U_recv.begin(), recv_cnt.begin(), recv_dsp.begin());
+        }
+
+        for (Long i = 0; i < recv_index.Dim(); i++) { // Set U
+          Long idx = (recv_index[i] - splitter_ranks[rank]) * dof * KDIM1_;
+          for (Integer k = 0; k < dof * KDIM1_; k++) {
+            U[idx + k] += U_recv[i*dof*KDIM1_ + k];
+          }
+        }
+      }
+    }
+
+
+
+    template <class ElemList, class DensityBasis, class Kernel> static void Direct(Vector<Real>& U, const Vector<Real>& Xt, const ElemList& elem_lst, const Vector<DensityBasis>& density, const Kernel& kernel, Integer order_direct, const Comm& comm) {
+      using CoordBasis = typename ElemList::CoordBasis;
+      using CoordEvalOpType = typename CoordBasis::EvalOpType;
+      using DensityEvalOpType = typename DensityBasis::EvalOpType;
+
+      constexpr Integer CoordDim = ElemList::CoordDim();
+      constexpr Integer ElemDim = ElemList::ElemDim();
+      constexpr Integer KDIM0 = Kernel::SrcDim();
+      constexpr Integer KDIM1 = Kernel::TrgDim();
+      const Long Nelem = elem_lst.NElem();
+      const Integer dof = density.Dim() / Nelem / KDIM0;
+      SCTL_ASSERT(density.Dim() == Nelem * dof * KDIM0);
+
+      Matrix<Real> quad_nds;
+      Vector<Real> quad_wts;
+      TensorProductGaussQuad<ElemDim>(quad_nds, quad_wts, order_direct);
+      const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds);
+      Integer Nnds = quad_wts.Dim();
+
+      const Vector<CoordBasis>& X = elem_lst.ElemVector();
+      Vector<CoordBasis> dX;
+      CoordBasis::Grad(dX, X);
+
+      Vector<Real> X_, dX_, Xa_, Xn_;
+      eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp);
+      eval_basis(dX_, dX, CoordDim*ElemDim, Nnds, CoordEvalOp);
+      if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_
+        Long N = Nelem*Nnds;
+        Xa_.ReInit(N);
+        Xn_.ReInit(N*CoordDim);
+        for (Long j = 0; j < N; j++) {
+          StaticArray<Real,CoordDim> normal;
+          normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3];
+          normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5];
+          normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1];
+          Xa_[j] = sctl::sqrt<Real>(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+          Real invXa = 1/Xa_[j];
+          Xn_[j*3+0] = normal[0] * invXa;
+          Xn_[j*3+1] = normal[1] * invXa;
+          Xn_[j*3+2] = normal[2] * invXa;
+        }
+      }
+
+      Vector<Real> Fa_;
+      { // Set Fa_
+        Vector<Real> F_;
+        if (std::is_same<CoordBasis,DensityBasis>::value) {
+          eval_basis(F_, density, dof * KDIM0, Nnds, CoordEvalOp);
+        } else {
+          const DensityEvalOpType EvalOp = DensityBasis::SetupEval(quad_nds);
+          eval_basis(F_, density, dof * KDIM0, Nnds, EvalOp);
+        }
+
+        Fa_.ReInit(F_.Dim());
+        const Integer DensityDOF = dof * KDIM0;
+        SCTL_ASSERT(F_.Dim() == Nelem * Nnds * DensityDOF);
+        for (Long j = 0; j < Nelem; j++) {
+          for (Integer k = 0; k < Nnds; k++) {
+            Long idx = j * Nnds + k;
+            Real quad_wt = Xa_[idx] * quad_wts[k];
+            for (Integer l = 0; l < DensityDOF; l++) {
+              Fa_[idx * DensityDOF + l] = F_[idx * DensityDOF + l] * quad_wt;
+            }
+          }
+        }
+      }
+
+      { // Evaluate potential
+        const Long Ntrg = Xt.Dim() / CoordDim;
+        SCTL_ASSERT(Xt.Dim() == Ntrg * CoordDim);
+        if (U.Dim() != Ntrg * dof * KDIM1) {
+          U.ReInit(Ntrg * dof * KDIM1);
+          U = 0;
+        }
+        ParticleFMM<Real,CoordDim>::Eval(U, Xt, X_, Xn_, Fa_, kernel, comm);
+      }
+    }
+
+  public:
+
+    template <class DensityBasis, class ElemList, class Kernel> void Setup(const ElemList& elem_lst, const Vector<Real>& Xt, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm) {
+      order_direct_ = order_direct;
+      period_length_ = period_length;
+      comm_ = comm;
+
+      Profile::Tic("Setup", &comm_);
+      static_assert(std::is_same<Real,typename DensityBasis::ValueType>::value);
+      static_assert(std::is_same<Real,typename ElemList::CoordType>::value);
+      static_assert(DensityBasis::Dim() == ElemList::ElemDim());
+
+      Xt_ = Xt;
+      M_singular.ReInit(0,0);
+
+      Profile::Tic("SetupNearSingular", &comm_);
+      SetupNearSingular<DensityBasis>(M_near_singular, pair_lst, Xt_, Vector<Long>(), elem_lst, kernel, order_singular, order_direct_, period_length_, comm_);
+      Profile::Toc();
+
+      Profile::Toc();
+    }
+
+    template <class DensityBasis, class PotentialBasis, class ElemList, class Kernel> void Setup(const ElemList& elem_lst, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm) {
+      order_direct_ = order_direct;
+      period_length_ = period_length;
+      comm_ = comm;
+
+      Profile::Tic("Setup", &comm_);
+      static_assert(std::is_same<Real,typename PotentialBasis::ValueType>::value);
+      static_assert(std::is_same<Real,typename DensityBasis::ValueType>::value);
+      static_assert(std::is_same<Real,typename ElemList::CoordType>::value);
+      static_assert(PotentialBasis::Dim() == ElemList::ElemDim());
+      static_assert(DensityBasis::Dim() == ElemList::ElemDim());
+
+      Vector<Long> trg_surf;
+      { // Set Xt_
+        using CoordBasis = typename ElemList::CoordBasis;
+        Matrix<Real> trg_nds = PotentialBasis::Nodes();
+        auto Meval = CoordBasis::SetupEval(trg_nds);
+        eval_basis(Xt_, elem_lst.ElemVector(), ElemList::CoordDim(), trg_nds.Dim(1), Meval);
+
+        { // Set trg_surf
+          const Long Nelem = elem_lst.NElem();
+          const Long Nnds  = trg_nds.Dim(1);
+          Long elem_offset;
+          { // Set elem_offset
+            comm.Scan(Ptr2ConstItr<Long>(&Nelem,1), Ptr2Itr<Long>(&elem_offset,1), 1, Comm::CommOp::SUM);
+            elem_offset -= Nelem;
+          }
+          trg_surf.ReInit(elem_lst.NElem() * trg_nds.Dim(1));
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnds; j++) {
+              trg_surf[i*Nnds+j] = elem_offset + i;
+            }
+          }
+        }
+      }
+
+      Profile::Tic("SetupSingular", &comm_);
+      SetupSingular<DensityBasis>(M_singular, PotentialBasis::Nodes(), elem_lst, kernel, order_singular, order_direct_);
+      Profile::Toc();
+
+      Profile::Tic("SetupNearSingular", &comm_);
+      SetupNearSingular<DensityBasis>(M_near_singular, pair_lst, Xt_, trg_surf, elem_lst, kernel, order_singular, order_direct_, period_length_, comm_);
+      Profile::Toc();
+
+      Profile::Toc();
+    }
+
+    template <class DensityBasis, class PotentialBasis, class ElemList, class Kernel> void Eval(Vector<PotentialBasis>& U, const ElemList& elements, const Vector<DensityBasis>& F, const Kernel& kernel) {
+      Profile::Tic("Eval", &comm_);
+      Matrix<Real> U_singular;
+      Vector<Real> U_direct, U_near_sing;
+
+      Profile::Tic("EvalDirect", &comm_);
+      Direct(U_direct, Xt_, elements, F, kernel, order_direct_, comm_);
+      Profile::Toc();
+
+      Profile::Tic("EvalSingular", &comm_);
+      EvalSingular(U_singular, F, M_singular, kernel.SrcDim(), kernel.TrgDim());
+      Profile::Toc();
+
+      Profile::Tic("EvalNearSingular", &comm_);
+      EvalNearSingular(U_near_sing, F, M_near_singular, pair_lst, elements.NElem(), Xt_.Dim() / ElemList::CoordDim(), kernel.SrcDim(), kernel.TrgDim(), comm_);
+      SCTL_ASSERT(U_near_sing.Dim() == U_direct.Dim());
+      Profile::Toc();
+
+      if (U.Dim() != elements.NElem() * kernel.TrgDim()) {
+        U.ReInit(elements.NElem() * kernel.TrgDim());
+      }
+      for (int i = 0; i < elements.NElem(); i++) {
+        for (int j = 0; j < PotentialBasis::Size(); j++) {
+          for (int k = 0; k < kernel.TrgDim(); k++) {
+            Real& U_ = U[i*kernel.TrgDim()+k][j];
+            U_ = 0;
+            U_ += U_direct   [(i*PotentialBasis::Size()+j)*kernel.TrgDim()+k];
+            U_ += U_near_sing[(i*PotentialBasis::Size()+j)*kernel.TrgDim()+k];
+            U_ += U_singular[i*kernel.TrgDim()+k][j];
+            U_ *= kernel.template ScaleFactor<Real>();
+          }
+        }
+      }
+      Profile::Toc();
+    }
+
+    template <class DensityBasis, class ElemList, class Kernel> void Eval(Vector<Real>& U, const ElemList& elements, const Vector<DensityBasis>& F, const Kernel& kernel) {
+      Profile::Tic("Eval", &comm_);
+      Matrix<Real> U_singular;
+      Vector<Real> U_direct, U_near_sing;
+
+      Profile::Tic("EvalDirect", &comm_);
+      Direct(U_direct, Xt_, elements, F, kernel, order_direct_, comm_);
+      Profile::Toc();
+
+      Profile::Tic("EvalSingular", &comm_);
+      EvalSingular(U_singular, F, M_singular, kernel.SrcDim(), kernel.TrgDim());
+      Profile::Toc();
+
+      Profile::Tic("EvalNearSingular", &comm_);
+      EvalNearSingular(U_near_sing, F, M_near_singular, pair_lst, elements.NElem(), Xt_.Dim() / ElemList::CoordDim(), kernel.SrcDim(), kernel.TrgDim(), comm_);
+      SCTL_ASSERT(U_near_sing.Dim() == U_direct.Dim());
+      Profile::Toc();
+
+      if (U.Dim() != U_direct.Dim()) {
+        U.ReInit(U_direct.Dim());
+      }
+      for (int i = 0; i < U.Dim(); i++) {
+        U[i] = (U_direct[i] + U_near_sing[i]) * kernel.template ScaleFactor<Real>();
+      }
+      if (U_singular.Dim(1)) {
+        for (int i = 0; i < elements.NElem(); i++) {
+          for (int j = 0; j < U_singular.Dim(1); j++) {
+            for (int k = 0; k < kernel.TrgDim(); k++) {
+              Real& U_ = U[(i*U_singular.Dim(1)+j)*kernel.TrgDim()+k];
+              U_ += U_singular[i*kernel.TrgDim()+k][j] * kernel.template ScaleFactor<Real>();
+            }
+          }
+        }
+      }
+      Profile::Toc();
+    }
+
+    template <Integer ORDER = 5> static void test(Integer order_singular = 10, Integer order_direct = 5, const Comm& comm = Comm::World()) {
+      constexpr Integer COORD_DIM = 3;
+      constexpr Integer ELEM_DIM = COORD_DIM-1;
+      using ElemList = ElemList<COORD_DIM, Basis<Real, ELEM_DIM, ORDER>>;
+      using DensityBasis = Basis<Real, ELEM_DIM, ORDER>;
+      using PotentialBasis = Basis<Real, ELEM_DIM, ORDER>;
+
+      int np = comm.Size();
+      int rank = comm.Rank();
+      auto build_torus = [rank,np](ElemList& elements, long Nt, long Np, Real Rmajor, Real Rminor){
+        auto nodes = ElemList::CoordBasis::Nodes();
+        auto torus = [](Real theta, Real phi, Real Rmajor, Real Rminor) {
+          Real R = Rmajor + Rminor * cos<Real>(phi);
+          Real X = R * cos<Real>(theta);
+          Real Y = R * sin<Real>(theta);
+          Real Z = Rminor * sin<Real>(phi);
+          return std::make_tuple(X,Y,Z);
+        };
+
+        long start = Nt*Np*(rank+0)/np;
+        long end   = Nt*Np*(rank+1)/np;
+        elements.ReInit(end - start);
+        for (long ii = start; ii < end; ii++) {
+          long i = ii / Np;
+          long j = ii % Np;
+          for (int k = 0; k < nodes.Dim(1); k++) {
+            Real X, Y, Z;
+            Real theta = 2 * const_pi<Real>() * (i + nodes[0][k]) / Nt;
+            Real phi   = 2 * const_pi<Real>() * (j + nodes[1][k]) / Np;
+            std::tie(X,Y,Z) = torus(theta, phi, Rmajor, Rminor);
+            elements(ii-start,0)[k] = X;
+            elements(ii-start,1)[k] = Y;
+            elements(ii-start,2)[k] = Z;
+          }
+        }
+      };
+      ElemList elements_src, elements_trg;
+      build_torus(elements_src, 28, 16, 2, 1.0);
+      build_torus(elements_trg, 29, 17, 2, 0.99);
+
+      Vector<Real> Xt;
+      Vector<PotentialBasis> U_onsurf, U_offsurf;
+      Vector<DensityBasis> density_sl, density_dl;
+      { // Set Xt, elements_src, elements_trg, density_sl, density_dl, U
+        Real X0[COORD_DIM] = {3,2,1};
+        std::function<void(Real*,Real*,Real*)> potential = [X0](Real* U, Real* X, Real* Xn) {
+          Real dX[COORD_DIM] = {X[0]-X0[0],X[1]-X0[1],X[2]-X0[2]};
+          Real Rinv = 1/sqrt(dX[0]*dX[0]+dX[1]*dX[1]+dX[2]*dX[2]);
+          U[0] = Rinv;
+        };
+        std::function<void(Real*,Real*,Real*)> potential_normal_derivative = [X0](Real* U, Real* X, Real* Xn) {
+          Real dX[COORD_DIM] = {X[0]-X0[0],X[1]-X0[1],X[2]-X0[2]};
+          Real Rinv = 1/sqrt(dX[0]*dX[0]+dX[1]*dX[1]+dX[2]*dX[2]);
+          Real RdotN = dX[0]*Xn[0]+dX[1]*Xn[1]+dX[2]*Xn[2];
+          U[0] = -RdotN * Rinv*Rinv*Rinv;
+        };
+
+        DiscretizeSurfaceFn<COORD_DIM,1>(density_sl, elements_src, potential_normal_derivative);
+        DiscretizeSurfaceFn<COORD_DIM,1>(density_dl, elements_src, potential);
+        DiscretizeSurfaceFn<COORD_DIM,1>(U_onsurf  , elements_src, potential);
+        DiscretizeSurfaceFn<COORD_DIM,1>(U_offsurf , elements_trg, potential);
+
+        for (long i = 0; i < elements_trg.NElem(); i++) { // Set Xt
+          for (long j = 0; j < PotentialBasis::Size(); j++) {
+            for (int k = 0; k < COORD_DIM; k++) {
+              Xt.PushBack(elements_trg(i,k)[j]);
+            }
+          }
+        }
+      }
+
+      GenericKernel<Laplace3D_DxU> Laplace_DxU;
+      GenericKernel<Laplace3D_FxU> Laplace_FxU;
+
+      Profile::Enable(true);
+      if (1) { // Greeen's identity test (Laplace, on-surface)
+        Profile::Tic("OnSurface", &comm);
+        Quadrature<Real> quadrature_DxU, quadrature_FxU;
+        quadrature_FxU.Setup<DensityBasis, PotentialBasis>(elements_src, Laplace_FxU, order_singular, order_direct, -1.0, comm);
+        quadrature_DxU.Setup<DensityBasis, PotentialBasis>(elements_src, Laplace_DxU, order_singular, order_direct, -1.0, comm);
+
+        Vector<PotentialBasis> U_sl, U_dl;
+        quadrature_FxU.Eval(U_sl, elements_src, density_sl, Laplace_FxU);
+        quadrature_DxU.Eval(U_dl, elements_src, density_dl, Laplace_DxU);
+        Profile::Toc();
+
+        Real max_err = 0;
+        Vector<PotentialBasis> err(U_onsurf.Dim());
+        for (long i = 0; i < U_sl.Dim(); i++) {
+          for (long j = 0; j < PotentialBasis::Size(); j++) {
+            err[i][j] = 0.5*U_onsurf[i][j] - (U_sl[i][j] + U_dl[i][j]);
+            max_err = std::max<Real>(max_err, fabs(err[i][j]));
+          }
+        }
+        { // Print error
+          Real glb_err;
+          comm.Allreduce(Ptr2ConstItr<Real>(&max_err,1), Ptr2Itr<Real>(&glb_err,1), 1, Comm::CommOp::MAX);
+          if (!comm.Rank()) std::cout<<"Error = "<<glb_err<<'\n';
+        }
+        { // Write VTK output
+          VTUData vtu;
+          vtu.AddElems(elements_src, err, ORDER);
+          vtu.WriteVTK("err", comm);
+        }
+        { // Write VTK output
+          VTUData vtu;
+          vtu.AddElems(elements_src, U_onsurf, ORDER);
+          vtu.WriteVTK("U", comm);
+        }
+      }
+      if (1) { // Greeen's identity test (Laplace, off-surface)
+        Profile::Tic("OffSurface", &comm);
+        Quadrature<Real> quadrature_DxU, quadrature_FxU;
+        quadrature_FxU.Setup<DensityBasis>(elements_src, Xt, Laplace_FxU, order_singular, order_direct, -1.0, comm);
+        quadrature_DxU.Setup<DensityBasis>(elements_src, Xt, Laplace_DxU, order_singular, order_direct, -1.0, comm);
+
+        Vector<Real> U_sl, U_dl;
+        quadrature_FxU.Eval(U_sl, elements_src, density_sl, Laplace_FxU);
+        quadrature_DxU.Eval(U_dl, elements_src, density_dl, Laplace_DxU);
+        Profile::Toc();
+
+        Real max_err = 0;
+        Vector<PotentialBasis> err(elements_trg.NElem());
+        for (long i = 0; i < elements_trg.NElem(); i++) {
+          for (long j = 0; j < PotentialBasis::Size(); j++) {
+            err[i][j] = U_offsurf[i][j] - (U_sl[i*PotentialBasis::Size()+j] + U_dl[i*PotentialBasis::Size()+j]);
+            max_err = std::max<Real>(max_err, fabs(err[i][j]));
+          }
+        }
+        { // Print error
+          Real glb_err;
+          comm.Allreduce(Ptr2ConstItr<Real>(&max_err,1), Ptr2Itr<Real>(&glb_err,1), 1, Comm::CommOp::MAX);
+          if (!comm.Rank()) std::cout<<"Error = "<<glb_err<<'\n';
+        }
+        { // Write VTK output
+          VTUData vtu;
+          vtu.AddElems(elements_trg, err, ORDER);
+          vtu.WriteVTK("err", comm);
+        }
+        { // Write VTK output
+          VTUData vtu;
+          vtu.AddElems(elements_trg, U_offsurf, ORDER);
+          vtu.WriteVTK("U", comm);
+        }
+      }
+      Profile::print(&comm);
+    }
+
+  private:
+
+    static void scan(Vector<Long>& dsp, const Vector<Long>& cnt) {
+      dsp.ReInit(cnt.Dim());
+      if (cnt.Dim()) dsp[0] = 0;
+      omp_par::scan(cnt.begin(), dsp.begin(), cnt.Dim());
+    }
+
+    template <class Basis> static void eval_basis(Vector<Real>& value, const Vector<Basis> X, Integer dof, Integer Nnds, const typename Basis::EvalOpType& EvalOp) {
+      Long Nelem = X.Dim() / dof;
+      SCTL_ASSERT(X.Dim() == Nelem * dof);
+
+      value.ReInit(Nelem*Nnds*dof);
+      Matrix<Real> X_(Nelem*dof, Nnds, value.begin(),false);
+      Basis::Eval(X_, X, EvalOp);
+      for (Long j = 0; j < Nelem; j++) { // Rearrange data
+        Matrix<Real> X(Nnds, dof, X_[j*dof], false);
+        X = Matrix<Real>(dof, Nnds, X_[j*dof], false).Transpose();
+      }
+    }
+
+    template <int CoordDim, int FnDim, class Real, class FnBasis, class ElemList> static void DiscretizeSurfaceFn(Vector<FnBasis>& U, const ElemList& elements, std::function<void(Real*,Real*,Real*)> fn) {
+      using CoordBasis = typename ElemList::CoordBasis;
+      const long Nelem = elements.NElem();
+      U.ReInit(Nelem * FnDim);
+
+      Matrix<Real> X, X_grad;
+      { // Set X, X_grad
+        Vector<CoordBasis> coord = elements.ElemVector();
+        Vector<CoordBasis> coord_grad;
+        CoordBasis::Grad(coord_grad, coord);
+
+        const auto Meval = CoordBasis::SetupEval(FnBasis::Nodes());
+        CoordBasis::Eval(X, coord, Meval);
+        CoordBasis::Eval(X_grad, coord_grad, Meval);
+      }
+
+      for (long i = 0; i < Nelem; i++) {
+        for (long j = 0; j < FnBasis::Size(); j++) {
+          Real X_[CoordDim], Xn[CoordDim], U_[FnDim];
+          for (long k = 0; k < CoordDim; k++) {
+            X_[k] = X[i*CoordDim+k][j];
+          }
+          { // Set Xn
+            Real Xu[CoordDim], Xv[CoordDim];
+            for (long k = 0; k < CoordDim; k++) {
+              Xu[k] = X_grad[(i*CoordDim+k)*2+0][j];
+              Xv[k] = X_grad[(i*CoordDim+k)*2+1][j];
+            }
+            Real dA = 0;
+            for (long k = 0; k < CoordDim; k++) {
+              Xn[k] = Xu[(k+1)%CoordDim] * Xv[(k+2)%CoordDim];
+              Xn[k] -= Xv[(k+1)%CoordDim] * Xu[(k+2)%CoordDim];
+              dA += Xn[k] * Xn[k];
+            }
+            dA = sqrt(dA);
+            for (long k = 0; k < CoordDim; k++) {
+              Xn[k] /= dA;
+            }
+          }
+          fn(U_, X_, Xn);
+          for (long k = 0; k < FnDim; k++) {
+            U[i*FnDim+k][j] = U_[k];
+          }
+        }
+      }
+    }
+
+    Vector<Real> Xt_;
+    Matrix<Real> M_singular;
+    Matrix<Real> M_near_singular;
+    Vector<Pair<Long,Long>> pair_lst;
+    Integer order_direct_;
+
+    Real period_length_;
+    Comm comm_;
+};
+
+}  // end namespace
+
+#endif  //_SCTL_BOUNDARY_QUADRATURE_HPP_

include/sctl/comm.txx

@@ -982,21 +982,22 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
   }
   srand(myrank);
 
-  Long totSize, nelem = arr_.Dim();
+  Long totSize;
   {                 // Local and global sizes. O(log p)
+    Long nelem = arr_.Dim();
     Allreduce<Long>(Ptr2ConstItr<Long>(&nelem, 1), Ptr2Itr<Long>(&totSize, 1), 1, CommOp::SUM);
   }
 
   if (npes == 1) {  // SortedElem <--- local_sort(arr_)
     SortedElem = arr_;
-    omp_par::merge_sort(SortedElem.begin(), SortedElem.begin() + nelem);
+    omp_par::merge_sort(SortedElem.begin(), SortedElem.end());
     return;
   }
 
   Vector<Type> arr;
   {  // arr <-- local_sort(arr_)
     arr = arr_;
-    omp_par::merge_sort(arr.begin(), arr.begin() + nelem);
+    omp_par::merge_sort(arr.begin(), arr.end());
   }
 
   Vector<Type> nbuff, nbuff_ext, rbuff, rbuff_ext;  // Allocate memory.
@@ -1012,6 +1013,7 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
       Vector<Type> glb_splitters;
       {  // Take random splitters. glb_splt_count = const = 100~1000
         Integer splt_count;
+        Long nelem = arr.Dim();
         {  // Set splt_coun. O( 1 ) -- Let p * splt_count = t
           splt_count = (100 * nelem) / totSize;
           if (npes > 100) splt_count = (drand48() * totSize) < (100 * nelem) ? 1 : 0;
@@ -1051,7 +1053,7 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
       {  // Compute local rank
 #pragma omp parallel for schedule(static)
         for (Integer i = 0; i < glb_splt_count; i++) {
-          lrank[i] = std::lower_bound(arr.begin(), arr.begin() + nelem, glb_splitters[i]) - arr.begin();
+          lrank[i] = std::lower_bound(arr.begin(), arr.end(), glb_splitters[i]) - arr.begin();
         }
       }
 
@@ -1061,20 +1063,20 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
       }
 
       {  // Determine split_key, totSize_new
-        ConstIterator<Long> split_disp = grank.begin();
+        Integer splitter_idx = 0;
         for (Integer i = 0; i < glb_splt_count; i++) {
-          if (labs(grank[i] - totSize / 2) < labs(*split_disp - totSize / 2)) {
-            split_disp = grank.begin() + i;
+          if (labs(grank[i] - totSize / 2) < labs(grank[splitter_idx] - totSize / 2)) {
+            splitter_idx = i;
           }
         }
-        split_key = glb_splitters[split_disp - grank.begin()];
+        split_key = glb_splitters[splitter_idx];
 
         if (myrank <= (npes - 1) / 2)
-          totSize_new = split_disp[0];
+          totSize_new = grank[splitter_idx];
         else
-          totSize_new = totSize - split_disp[0];
+          totSize_new = totSize - grank[splitter_idx];
 
-        // double err=(((double)*split_disp)/(totSize/2))-1.0;
+        // double err=(((double)grank[splitter_idx])/(totSize/2))-1.0;
         // if(fabs<double>(err)<0.01 || npes<=16) break;
         // else if(!myrank) std::cout<<err<<'\n';
       }
@@ -1082,24 +1084,21 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
 
     Integer split_id = (npes - 1) / 2;
     {  // Split problem into two. O( N/p )
-      Integer new_p0 = (myrank <= split_id ? 0 : split_id + 1);
-      Integer cmp_p0 = (myrank > split_id ? 0 : split_id + 1);
-
       Integer partner;
       {  // Set partner
-        partner = myrank + cmp_p0 - new_p0;
+        partner = myrank + (split_id+1) * (myrank<=split_id ? 1 : -1);
         if (partner >= npes) partner = npes - 1;
         assert(partner >= 0);
       }
-      bool extra_partner = (npes % 2 == 1 && npes - 1 == myrank);
+      bool extra_partner = (npes % 2 == 1 && myrank == npes - 1);
 
       Long ssize = 0, lsize = 0;
       ConstIterator<Type> sbuff, lbuff;
       {  // Set ssize, lsize, sbuff, lbuff
-        Long split_indx = std::lower_bound(arr.begin(), arr.begin() + nelem, split_key) - arr.begin();
-        ssize = (myrank > split_id ? split_indx : nelem - split_indx);
+        Long split_indx = std::lower_bound(arr.begin(), arr.end(), split_key) - arr.begin();
+        ssize = (myrank > split_id ? split_indx : arr.Dim() - split_indx);
         sbuff = (myrank > split_id ? arr.begin() : arr.begin() + split_indx);
-        lsize = (myrank <= split_id ? split_indx : nelem - split_indx);
+        lsize = (myrank <= split_id ? split_indx : arr.Dim() - split_indx);
         lbuff = (myrank <= split_id ? arr.begin() : arr.begin() + split_indx);
       }
 
@@ -1119,23 +1118,24 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
         if (extra_partner) MPI_Sendrecv(nullptr, 0, CommDatatype<Type>::value(), split_id, 0, (ext_rsize ? &rbuff_ext[0] : nullptr), ext_rsize, CommDatatype<Type>::value(), split_id, 0, comm, &status);
       }
 
-      Long nbuff_size = lsize + rsize + ext_rsize;
       {  // nbuff <-- merge(lbuff, rbuff, rbuff_ext)
         nbuff.ReInit(lsize + rsize);
         omp_par::merge<ConstIterator<Type>>(lbuff, (lbuff + lsize), rbuff.begin(), rbuff.begin() + rsize, nbuff.begin(), omp_p, std::less<Type>());
-        if (ext_rsize > 0 && nbuff.Dim() > 0) {
-          nbuff_ext.ReInit(nbuff_size);
-          omp_par::merge(nbuff.begin(), nbuff.begin() + (lsize + rsize), rbuff_ext.begin(), rbuff_ext.begin() + ext_rsize, nbuff_ext.begin(), omp_p, std::less<Type>());
-          nbuff.Swap(nbuff_ext);
-          nbuff_ext.ReInit(0);
+        if (ext_rsize > 0) {
+          if (nbuff.Dim() > 0) {
+            nbuff_ext.ReInit(lsize + rsize + ext_rsize);
+            omp_par::merge(nbuff.begin(), nbuff.begin() + (lsize + rsize), rbuff_ext.begin(), rbuff_ext.begin() + ext_rsize, nbuff_ext.begin(), omp_p, std::less<Type>());
+            nbuff.Swap(nbuff_ext);
+            nbuff_ext.ReInit(0);
+          } else {
+            nbuff.Swap(rbuff_ext);
+          }
         }
       }
 
       // Copy new data.
       totSize = totSize_new;
-      nelem = nbuff_size;
       arr.Swap(nbuff);
-      nbuff.ReInit(0);
     }
 
     {  // Split comm.  O( log(p) ) ??

include/sctl/kernel_functions.hpp

@@ -0,0 +1,228 @@
+#ifndef _SCTL_KERNEL_FUNCTIONS_HPP_
+#define _SCTL_KERNEL_FUNCTIONS_HPP_
+
+namespace SCTL_NAMESPACE {
+
+struct Laplace3D_FxU {
+  template <class Real> static constexpr Real ScaleFactor() {
+    return 1 / (4 * const_pi<Real>());
+  }
+  template <class Real> static void Eval(Real (&u)[1][1], const Real (&r)[3], const Real (&n)[3], void* ctx_ptr) {
+    Real r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
+    Real rinv = (r2>1e-16 ? 1/sqrt<Real>(r2) : 0);
+    u[0][0] = rinv;
+  }
+};
+struct Laplace3D_DxU {
+  template <class Real> static constexpr Real ScaleFactor() {
+    return 1 / (4 * const_pi<Real>());
+  }
+  template <class Real> static void Eval(Real (&u)[1][1], const Real (&r)[3], const Real (&n)[3], void* ctx_ptr) {
+    Real r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
+    Real rinv = (r2>1e-16 ? 1/sqrt<Real>(r2) : 0);
+    Real rdotn = r[0]*n[0] + r[1]*n[1] + r[2]*n[2];
+    Real rinv3 = rinv * rinv * rinv;
+    u[0][0] = rdotn * rinv3;
+  }
+};
+struct Laplace3D_FxdU{
+  template <class Real> static constexpr Real ScaleFactor() {
+    return 1 / (4 * const_pi<Real>());
+  }
+  template <class Real> static void Eval(Real (&u)[1][3], const Real (&r)[3], const Real (&n)[3], void* ctx_ptr) {
+    Real r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
+    Real rinv = (r2>1e-16 ? 1/sqrt<Real>(r2) : 0);
+    Real rinv3 = rinv * rinv * rinv;
+    u[0][0] = -r[0] * rinv3;
+    u[0][1] = -r[1] * rinv3;
+    u[0][2] = -r[2] * rinv3;
+  }
+};
+
+template <class uKernel> class GenericKernel {
+
+    template <class Real, Integer D, Integer K0, Integer K1> static constexpr Integer get_DIM  (void (*uKer)(Real (&u)[K0][K1], const Real (&r)[D], const Real (&n)[D], void* ctx_ptr)) { return D; }
+    template <class Real, Integer D, Integer K0, Integer K1> static constexpr Integer get_KDIM0(void (*uKer)(Real (&u)[K0][K1], const Real (&r)[D], const Real (&n)[D], void* ctx_ptr)) { return K0; }
+    template <class Real, Integer D, Integer K0, Integer K1> static constexpr Integer get_KDIM1(void (*uKer)(Real (&u)[K0][K1], const Real (&r)[D], const Real (&n)[D], void* ctx_ptr)) { return K1; }
+
+    static constexpr Integer DIM   = get_DIM  (uKernel::template Eval<double>);
+    static constexpr Integer KDIM0 = get_KDIM0(uKernel::template Eval<double>);
+    static constexpr Integer KDIM1 = get_KDIM1(uKernel::template Eval<double>);
+
+  public:
+
+    static constexpr Integer CoordDim() {
+      return DIM;
+    }
+    static constexpr Integer SrcDim() {
+      return KDIM0;
+    }
+    static constexpr Integer TrgDim() {
+      return KDIM1;
+    }
+    template <class Real> static constexpr Real ScaleFactor() {
+      return uKernel::template ScaleFactor<Real>();
+    }
+
+    template <class Real, Integer DOF> void uKernelEval(Iterator<Real> u, ConstIterator<Real> r, ConstIterator<Real> n, ConstIterator<Real> f) const {
+      Real M[KDIM0][KDIM1];
+      uKernel::Eval(M, *(Real (*)[DIM])r, *(Real (*)[DIM])n, ctx_ptr);
+      for (Integer i = 0; i < DOF; i++) {
+        for (Integer k0 = 0; k0 < KDIM0; k0++) {
+          for (Integer k1 = 0; k1 < KDIM1; k1++) {
+            u[i*KDIM1+k1] += M[k0][k1] * f[i*KDIM0+k0];
+          }
+        }
+      }
+    }
+
+    template <class Real, Integer DOF_ = 0> void Eval(Vector<Real>& U, const Vector<Real>& Xt, const Vector<Real>& Xs, const Vector<Real>& Xn, const Vector<Real>& F) const {
+      Long Ns = Xs.Dim() / DIM;
+      Long Nt = Xt.Dim() / DIM;
+      Long DOF = (DOF_ ? DOF_ : F.Dim() / Ns / KDIM0);
+      SCTL_ASSERT(Xs.Dim() == Ns * DIM);
+      SCTL_ASSERT(Xt.Dim() == Nt * DIM);
+      SCTL_ASSERT(F.Dim() == Ns * DOF * KDIM0);
+
+      if (!DOF_) {
+        if (DOF == 1) Eval<Real,1>(U,Xt,Xs,Xn,F);
+        if (DOF == 2) Eval<Real,2>(U,Xt,Xs,Xn,F);
+        if (DOF == 3) Eval<Real,3>(U,Xt,Xs,Xn,F);
+        if (DOF == 4) Eval<Real,4>(U,Xt,Xs,Xn,F);
+        if (DOF == 5) Eval<Real,5>(U,Xt,Xs,Xn,F);
+        if (DOF == 6) Eval<Real,6>(U,Xt,Xs,Xn,F);
+        if (DOF <= 6) return;
+      }
+
+      if (U.Dim() != Nt * DOF * KDIM1) {
+        U.ReInit(Nt * DOF * KDIM1);
+        U.SetZero();
+      }
+      #pragma omp parallel for schedule(static)
+      for (Long t = 0; t < Nt; t++) {
+        Iterator<Real> u;
+        Vector<Real> u_buff0;
+        StaticArray<Real,DOF_*KDIM1> u_buff1;
+        if (!DOF_) { // Set u
+          u_buff0.ReInit(DOF*KDIM1);
+          u = u_buff0.begin();
+        } else {
+          u = (Iterator<Real>)u_buff1;
+        }
+        for (Integer k = 0; k < DOF*KDIM1; k++) {
+          u[k] = 0;
+        }
+
+        StaticArray<Real,DIM> r;
+        for (Long s = 0; s < Ns; s++) {
+          for (Integer k = 0; k < DIM; k++) {
+            r[k] = Xs[s*DIM+k] - Xt[t*DIM+k];
+          }
+          ConstIterator<Real> n = Xn.begin() + s*DIM;
+          if (DOF_) {
+            ConstIterator<Real> f = F.begin() + s*DOF_*KDIM0;
+            this->template uKernelEval<Real,DOF_>(u,r,n,f);
+          } else {
+            ConstIterator<Real> f = F.begin() + s*DOF*KDIM0;
+            for (Integer k = 0; k < DOF; k++) {
+              this->template uKernelEval<Real,1>(u + k*KDIM1,r,n,f + k*KDIM0);
+            }
+          }
+        }
+        for (Integer k = 0; k < DOF*KDIM1; k++) {
+          U[t*DOF*KDIM1 + k] += u[k];
+        }
+      }
+    }
+
+    template <class Real> void KernelMatrix(Matrix<Real>& M, const Vector<Real>& Xt, const Vector<Real>& Xs, const Vector<Real>& Xn) const {
+      Long Ns = Xs.Dim() / DIM;
+      Long Nt = Xt.Dim() / DIM;
+      SCTL_ASSERT(Xs.Dim() == Ns * DIM);
+      SCTL_ASSERT(Xt.Dim() == Nt * DIM);
+
+      if (M.Dim(0) != Ns * KDIM0 || M.Dim(1) != Nt * KDIM1) {
+        M.ReInit(Ns * KDIM0, Nt * KDIM1);
+      }
+      #pragma omp parallel for schedule(static)
+      for (Long t = 0; t < Nt; t++) {
+        StaticArray<Real,DIM> r;
+        StaticArray<Real,KDIM0> f;
+        for (Integer k = 0; k < KDIM0; k++) {
+          f[k] = 0;
+        }
+        for (Long s = 0; s < Ns; s++) {
+          for (Integer k = 0; k < DIM; k++) {
+            r[k] = Xs[s*DIM+k] - Xt[t*DIM+k];
+          }
+          ConstIterator<Real> n = Xn.begin() + s*DIM;
+
+          for (Integer k = 0; k < KDIM0; k++) {
+            f[k] = 1;
+            Iterator<Real> u = M[s*KDIM0+k] + t*KDIM1;
+            for (Integer i = 0; i < KDIM1; i++) u[i] = 0;
+            this->template uKernelEval<Real,1>(u,r,n,f);
+            f[k] = 0;
+          }
+        }
+      }
+    }
+
+  private:
+    void* ctx_ptr;
+};
+
+template <class Real, Integer DIM> class ParticleFMM {
+  public:
+
+    template <class Kernel> static void Eval(Vector<Real>& U, const Vector<Real>& Xt, const Vector<Real>& Xs, const Vector<Real>& Xn, const Vector<Real>& F, const Kernel& kernel, const Comm& comm) {
+      #ifdef SCTL_HAVE_PVFMM
+      SCTL_ASSERT(false);
+      #else
+      constexpr Integer KDIM0 = Kernel::SrcDim();
+      constexpr Integer KDIM1 = Kernel::TrgDim();
+      const Integer rank = comm.Rank();
+      const Integer np = comm.Size();
+      const Long Ns = Xs.Dim() / DIM;
+      const Long Nt = Xt.Dim() / DIM;
+      SCTL_ASSERT(Xt.Dim() == Nt * DIM);
+
+      const Long dof = F.Dim() / (Ns * KDIM0);
+      SCTL_ASSERT(Ns && F.Dim() == Ns * dof * KDIM0);
+
+      Vector<Real> Xs_, Xn_, F_;
+      if (U.Dim() != Nt * dof * KDIM1) {
+        U.ReInit(Nt * dof * KDIM1);
+        U.SetZero();
+      }
+      for (Long i = 0; i < np; i++) {
+        auto send_recv_vec = [comm,rank,np](Vector<Real>& X_, const Vector<Real>& X, Integer offset){
+          Integer send_partner = (rank + offset) % np;
+          Integer recv_partner = (rank + np - offset) % np;
+
+          Long send_cnt = X.Dim(), recv_cnt;
+          void* recv_req = comm.Irecv(     Ptr2Itr<Long>(&recv_cnt,1), 1, recv_partner, offset);
+          void* send_req = comm.Isend(Ptr2ConstItr<Long>(&send_cnt,1), 1, send_partner, offset);
+          comm.Wait(recv_req);
+          comm.Wait(send_req);
+
+          X_.ReInit(recv_cnt);
+          recv_req = comm.Irecv(X_.begin(), recv_cnt, recv_partner, offset);
+          send_req = comm.Isend(X .begin(), send_cnt, send_partner, offset);
+          comm.Wait(recv_req);
+          comm.Wait(send_req);
+        };
+        send_recv_vec(Xs_, Xs, i);
+        send_recv_vec(Xn_, Xn, i);
+        send_recv_vec(F_ , F , i);
+        kernel.Eval(U, Xt, Xs_, Xn_, F_);
+      }
+      #endif
+    }
+
+  private:
+};
+
+}  // end namespace
+
+#endif  //_SCTL_KERNEL_FUNCTIONS_HPP_

include/sctl/morton.hpp

@@ -27,6 +27,10 @@ template <Integer DIM = 3> class Morton {
 
   static constexpr Integer MAX_DEPTH = SCTL_MAX_DEPTH;
 
+  static constexpr Integer MaxDepth() {
+    return MAX_DEPTH;
+  }
+
   Morton() {
     depth = 0;
     for (Integer i = 0; i < DIM; i++) x[i] = 0;
@@ -85,10 +89,14 @@ template <Integer DIM = 3> class Morton {
   }
 
   void NbrList(Vector<Morton> &nlst, uint8_t level, bool periodic) const {
-    nlst.ReInit(1);
+    static constexpr Integer MAX_NBRS = sctl::pow<DIM,Integer>(3);
+    StaticArray<Morton<DIM>,MAX_NBRS> nbrs;
+    Integer Nnbrs = 0;
+
     UINT_T mask = ~((((UINT_T)1) << (MAX_DEPTH - level)) - 1);
-    for (Integer i = 0; i < DIM; i++) nlst[0].x[i] = x[i] & mask;
-    nlst[0].depth = level;
+    for (Integer i = 0; i < DIM; i++) nbrs[0].x[i] = x[i] & mask;
+    nbrs[0].depth = level;
+    Nnbrs++;
 
     Morton m;
     Integer k = 1;
@@ -96,34 +104,45 @@ template <Integer DIM = 3> class Morton {
     if (periodic) {
       for (Integer i = 0; i < DIM; i++) {
         for (Integer j = 0; j < k; j++) {
-          m = nlst[j];
+          m = nbrs[j];
           m.x[i] = (m.x[i] + mask) & (maxCoord - 1);
-          nlst.PushBack(m);
+          nbrs[Nnbrs] = m;
+          Nnbrs++;
         }
         for (Integer j = 0; j < k; j++) {
-          m = nlst[j];
+          m = nbrs[j];
           m.x[i] = (m.x[i] - mask) & (maxCoord - 1);
-          nlst.PushBack(m);
+          nbrs[Nnbrs] = m;
+          Nnbrs++;
         }
-        k = nlst.Dim();
+        k = Nnbrs;
       }
     } else {
       for (Integer i = 0; i < DIM; i++) {
         for (Integer j = 0; j < k; j++) {
-          m = nlst[j];
+          m = nbrs[j];
           if (m.x[i] + mask < maxCoord) {
             m.x[i] += mask;
-            nlst.PushBack(m);
+            nbrs[Nnbrs] = m;
+            Nnbrs++;
           }
         }
         for (Integer j = 0; j < k; j++) {
-          m = nlst[j];
+          m = nbrs[j];
           if (m.x[i] >= mask) {
             m.x[i] -= mask;
-            nlst.PushBack(m);
+            nbrs[Nnbrs] = m;
+            Nnbrs++;
           }
         }
-        k = nlst.Dim();
+        k = Nnbrs;
+      }
+    }
+    if (nlst.Dim() != Nnbrs) {
+      nlst.ReInit(Nnbrs, nbrs);
+    } else {
+      for (Integer i = 0; i < Nnbrs; i++) {
+        nlst[i] = nbrs[i];
       }
     }
   }

include/sctl/sph_harm.hpp

@@ -493,7 +493,7 @@ template <class Real> class SphericalHarmonics{
     }
 };
 
-template class SphericalHarmonics<double>;
+//template class SphericalHarmonics<double>;
 
 }  // end namespace
 

include/sctl/tensor.hpp

@@ -0,0 +1,226 @@
+#ifndef _SCTL_TENSOR_HPP_
+#define _SCTL_TENSOR_HPP_
+
+#include SCTL_INCLUDE(mem_mgr.hpp)
+#include SCTL_INCLUDE(common.hpp)
+
+#include <iostream>
+
+namespace SCTL_NAMESPACE {
+
+template <class ValueType, bool own_data, Long... Args> class Tensor {
+
+    template <Long k> static constexpr Long SizeHelper() {
+      return 1;
+    }
+    template <Long k, Long d, Long... dd> static constexpr Long SizeHelper() {
+      return (k >= 0 ? d : 1) * SizeHelper<k+1, dd...>();
+    }
+
+    template <Long k> static constexpr Long DimHelper() {
+      return 1;
+    }
+    template <Long k, Long d0, Long... dd> static constexpr Long DimHelper() {
+      return k==0 ? d0 : DimHelper<k-1,dd...>();
+    }
+
+    template <typename T> static constexpr Long OrderHelper() {
+      return 0;
+    }
+    template <typename T, Long d, Long... dd> static constexpr Long OrderHelper() {
+      return 1 + OrderHelper<void, dd...>();
+    }
+
+    template <Long k, bool own_data_, Long... dd> struct RotateType {
+      using Value = Tensor<ValueType,own_data_,dd...>;
+    };
+    template <bool own_data_, Long d, Long... dd> struct RotateType<0,own_data_,d,dd...> {
+      using Value = Tensor<ValueType,own_data_,d,dd...>;
+    };
+    template <Long k, bool own_data_, Long d, Long... dd> struct RotateType<k,own_data_,d,dd...> {
+      using Value = typename RotateType<k-1,own_data_,dd...,d>::Value;
+    };
+
+  public:
+
+    static constexpr Long Order() {
+      return OrderHelper<void, Args...>();
+    }
+
+    static constexpr Long Size() {
+      return SizeHelper<0,Args...>();
+    }
+
+    template <Long k> static constexpr Long Dim() {
+      return DimHelper<k,Args...>();
+    }
+
+
+    Tensor(Iterator<ValueType> src_iter = NullIterator<ValueType>()) {
+      Init(src_iter);
+    }
+
+    Tensor(const Tensor &M) {
+      Init((Iterator<ValueType>)M.begin());
+    }
+
+    template <bool own_data_> Tensor(const Tensor<ValueType,own_data_,Args...> &M) {
+      Init((Iterator<ValueType>)M.begin());
+    }
+
+    Tensor &operator=(const Tensor &M) {
+      memcopy(begin(), M.begin(), Size());
+      return *this;
+    }
+
+
+    Iterator<ValueType> begin() {
+      return own_data ? (Iterator<ValueType>)buff : iter_[0];
+    }
+
+    ConstIterator<ValueType> begin() const {
+      return own_data ? (ConstIterator<ValueType>)buff : (ConstIterator<ValueType>)iter_[0];
+    }
+
+    Iterator<ValueType> end() {
+      return begin() + Size();
+    }
+
+    ConstIterator<ValueType> end() const {
+      return begin() + Size();
+    }
+
+
+    template <class ...PackedLong> ValueType& operator()(PackedLong... ii) {
+      return begin()[offset<0>(ii...)];
+    }
+
+    template <class ...PackedLong> ValueType operator()(PackedLong... ii) const {
+      return begin()[offset<0>(ii...)];
+    }
+
+
+    typename RotateType<1,true,Args...>::Value RotateLeft() const {
+      typename RotateType<1,true,Args...>::Value Tr;
+      const auto& T = *this;
+
+      constexpr Long N0 = Dim<0>();
+      constexpr Long N1 = Size() / N0;
+      for (Long i = 0; i < N0; i++) {
+        for (Long j = 0; j < N1; j++) {
+          Tr.begin()[j*N0+i] = T.begin()[i*N1+j];
+        }
+      }
+      return Tr;
+    }
+
+    typename RotateType<Order()-1,true,Args...>::Value RotateRight() const {
+      typename RotateType<Order()-1,true,Args...>::Value Tr;
+      const auto& T = *this;
+
+      constexpr Long N0 = Dim<Order()-1>();
+      constexpr Long N1 = Size() / N0;
+      for (Long i = 0; i < N0; i++) {
+        for (Long j = 0; j < N1; j++) {
+          Tr.begin()[i*N1+j] = T.begin()[j*N0+i];
+        }
+      }
+      return Tr;
+    }
+
+
+    Tensor<ValueType, true, Args...> operator*(const ValueType &s) const {
+      Tensor<ValueType, true, Args...> M0;
+      const auto &M1 = *this;
+
+      for (Long i = 0; i < Size(); i++) {
+          M0.begin()[i] = M1.begin()[i]*s;
+      }
+      return M0;
+    }
+
+    template <bool own_data_> Tensor<ValueType, true, Args...> operator+(const Tensor<ValueType, own_data_, Args...> &M2) const {
+      Tensor<ValueType, true, Args...> M0;
+      const auto &M1 = *this;
+
+      for (Long i = 0; i < Size(); i++) {
+          M0.begin()[i] = M1.begin()[i] + M2.begin()[i];
+      }
+      return M0;
+    }
+
+    template <bool own_data_> Tensor<ValueType, true, Args...> operator-(const Tensor<ValueType, own_data_, Args...> &M2) const {
+      Tensor<ValueType, true, Args...> M0;
+      const auto &M1 = *this;
+
+      for (Long i = 0; i < Size(); i++) {
+          M0.begin()[i] = M1.begin()[i] - M2.begin()[i];
+      }
+      return M0;
+    }
+
+    template <bool own_data_, Long N1, Long N2> Tensor<ValueType, true, Dim<0>(), N2> operator*(const Tensor<ValueType, own_data_, N1, N2> &M2) const {
+      static_assert(Order() == 2, "Multiplication is only defined for tensors of order two.");
+      static_assert(Dim<1>() == N1, "Tensor dimensions dont match for multiplication.");
+      Tensor<ValueType, true, Dim<0>(), N2> M0;
+      const auto &M1 = *this;
+
+      for (Long i = 0; i < Dim<0>(); i++) {
+        for (Long j = 0; j < N2; j++) {
+          ValueType Mij = 0;
+          for (Long k = 0; k < N1; k++) {
+            Mij += M1(i,k)*M2(k,j);
+          }
+          M0(i,j) = Mij;
+        }
+      }
+      return M0;
+    }
+
+  private:
+
+    template <Integer k> static Long offset() {
+      return 0;
+    }
+    template <Integer k, class ...PackedLong> static Long offset(Long i, PackedLong... ii) {
+      return i * SizeHelper<-(k+1),Args...>() + offset<k+1>(ii...);
+    }
+
+    void Init(Iterator<ValueType> src_iter) {
+      if (own_data) {
+        if (src_iter != NullIterator<ValueType>()) {
+          memcopy((Iterator<ValueType>)buff, src_iter, Size());
+        }
+      } else {
+        if (Size()) {
+          SCTL_UNUSED(src_iter[0]);
+          SCTL_UNUSED(src_iter[Size()-1]);
+          iter_[0] = Ptr2Itr<ValueType>(&src_iter[0], Size());
+        } else {
+          iter_[0] = NullIterator<ValueType>();
+        }
+      }
+    }
+
+    StaticArray<ValueType,own_data?Size():0> buff;
+    StaticArray<Iterator<ValueType>,own_data?0:1> iter_;
+};
+
+template <class ValueType, bool own_data, Long N1, Long N2> std::ostream& operator<<(std::ostream &output, const Tensor<ValueType, own_data, N1, N2> &M) {
+  std::ios::fmtflags f(std::cout.flags());
+  output << std::fixed << std::setprecision(4) << std::setiosflags(std::ios::left);
+  for (Long i = 0; i < N1; i++) {
+    for (Long j = 0; j < N2; j++) {
+      float f = ((float)M(i,j));
+      if (sctl::fabs<float>(f) < 1e-25) f = 0;
+      output << std::setw(10) << ((double)f) << ' ';
+    }
+    output << ";\n";
+  }
+  std::cout.flags(f);
+  return output;
+}
+
+}  // end namespace
+
+#endif  //_SCTL_TENSOR_HPP_

include/sctl/tree.hpp

@@ -833,6 +833,32 @@ template <Integer DIM> class Tree {
         }
       }
 
+      Broadcast<ValueType>(name);
+    }
+
+    template <class ValueType> void Broadcast(const std::string& name) {
+      Integer np = comm.Size();
+      Integer rank = comm.Rank();
+
+      Vector<Long> dsp;
+      Iterator<Vector<char>> data_;
+      Iterator<Vector<Long>> cnt_;
+      GetData_(data_, cnt_, name);
+      Vector<ValueType> data(data_->Dim()/sizeof(ValueType), (Iterator<ValueType>)data_->begin(), false);
+      Vector<Long>& cnt = *cnt_;
+      scan(dsp, cnt);
+
+      Long dof;
+      { // Set dof
+        StaticArray<Long,2> Nl, Ng;
+        Nl[0] = data.Dim();
+        Nl[1] = omp_par::reduce(cnt.begin(), cnt.Dim());
+        comm.Allreduce((ConstIterator<Long>)Nl, (Iterator<Long>)Ng, 2, Comm::CommOp::SUM);
+        dof = Ng[0] / std::max<Long>(Ng[1],1);
+        SCTL_ASSERT(Nl[0] == Nl[1] * dof);
+        SCTL_ASSERT(Ng[0] == Ng[1] * dof);
+      }
+
       { // Broadcast
         const Vector<Morton<DIM>>& send_mid = user_mid;
         const Vector<Long>& send_node_cnt = user_cnt;