1 files changed, 79 insertions, 0 deletions
diff --git a/extern/mantaflow/helper/util/integrator.h b/extern/mantaflow/helper/util/integrator.h
new file mode 100644
index 00000000000..5b1b02a5197
--- /dev/null
+++ b/extern/mantaflow/helper/util/integrator.h
@@ -0,0 +1,79 @@
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2011 Tobias Pfaff, Nils Thuerey
+ *
+ * This program is free software, distributed under the terms of the
+ * Apache License, Version 2.0
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Helper functions for simple integration
+ *
+ ******************************************************************************/
+
+#ifndef _INTEGRATE_H
+#define _INTEGRATE_H
+
+#include <vector>
+#include "vectorbase.h"
+#include "kernel.h"
+
+namespace Manta {
+
+enum IntegrationMode { IntEuler = 0, IntRK2, IntRK4 };
+
+//! Integrate a particle set with a given velocity kernel
+template<class VelKernel> void integratePointSet(VelKernel &k, int mode)
+{
+  typedef typename VelKernel::type0 PosType;
+  PosType &x = k.getArg0();
+  const std::vector<Vec3> &u = k.getRet();
+  const int N = x.size();
+
+  if (mode == IntEuler) {
+    for (int i = 0; i < N; i++)
+      x[i].pos += u[i];
+  }
+  else if (mode == IntRK2) {
+    PosType x0(x);
+
+    for (int i = 0; i < N; i++)
+      x[i].pos = x0[i].pos + 0.5 * u[i];
+
+    k.run();
+    for (int i = 0; i < N; i++)
+      x[i].pos = x0[i].pos + u[i];
+  }
+  else if (mode == IntRK4) {
+    PosType x0(x);
+    std::vector<Vec3> uTotal(u);
+
+    for (int i = 0; i < N; i++)
+      x[i].pos = x0[i].pos + 0.5 * u[i];
+
+    k.run();
+    for (int i = 0; i < N; i++) {
+      x[i].pos = x0[i].pos + 0.5 * u[i];
+      uTotal[i] += 2 * u[i];
+    }
+
+    k.run();
+    for (int i = 0; i < N; i++) {
+      x[i].pos = x0[i].pos + u[i];
+      uTotal[i] += 2 * u[i];
+    }
+
+    k.run();
+    for (int i = 0; i < N; i++)
+      x[i].pos = x0[i].pos + (Real)(1. / 6.) * (uTotal[i] + u[i]);
+  }
+  else
+    errMsg("unknown integration type");
+
+  // for(int i=0; i<N; i++) std::cout << x[i].pos.y-x[0].pos.y << std::endl;
+  // std::cout << "<><><>" << std::endl;
+}
+
+}  // namespace Manta
+
+#endif