Public Helmholtz FMM

Box.hpp

@@ -0,0 +1,53 @@
+#ifndef BOX_H
+#define BOX_H
+
+#include "General.hpp"
+#include "Vec3.hpp"
+#include "NFunction.hpp"
+
+class Box
+{
+  NFunction M;             // Multipole (outgoing) expansion
+  NFunction L;             // Local (incoming) expansion
+
+ public:
+  typedef vector<int> PointList;
+  typedef PointList::iterator pointIter;
+
+  int n;                   // The box number
+  Vec3 center;             // The box center
+
+  Box* parent;             // A pointer to this box's parent
+  PointList pointIndex;    // If leaf, contains indices of points it contains
+
+  Box() {}
+  Box( int n_, const Vec3& center_ ) : n(n_), center(center_) {}
+  ~Box() {}
+
+  inline void addPointIndex( int index ) {
+    pointIndex.push_back( index );
+  }
+
+  inline void makeMultipole( Quadrature* q ) {
+    M = NFunction(q);
+  }
+
+  inline void makeLocal( Quadrature* q ) {
+    L = NFunction(q);
+  }
+
+  inline NFunction& getMultipole() {
+    return M;
+  }
+
+  inline NFunction& getLocal() {
+    return L;
+  }
+
+ private:
+  // Disable Copy and Assignment
+  Box(const Box& S) { (void) S; }
+  void operator=(const Box& S) { (void) S; }
+};
+
+#endif

Direct.hpp

@@ -0,0 +1,21 @@
+#ifndef DIRECT_H
+#define DIRECT_H
+
+template <typename FUNC>
+void Direct(FUNC K, const vector<Vec3>& p,
+	    const vector<complex>& psi, vector<complex>& omega)
+{
+  int N = p.size();
+  for( int m = 0; m < N; ++m ) {
+    const Vec3& pm = p[m];
+    omega[m] = 0;
+    for( int n = 0; n != m; ++n ) {
+      double r = pm.dist(p[n]);
+      complex E = K(r);
+      omega[m] += psi[n] * E;
+      omega[n] += psi[m] * E;
+    }
+  }
+}
+
+#endif