Mantaflow [Part 1]: Added preprocessed Mantaflow source files

Includes preprocessed Mantaflow source files for both OpenMP and TBB (if OpenMP is not present, TBB files will be used instead).

These files come directly from the Mantaflow repository. Future updates to the core fluid solver will take place by updating the files.

Reviewed By: sergey, mont29

Maniphest Tasks: T59995

Differential Revision: https://developer.blender.org/D3850

15 files changed, 4977 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
diff --git a/extern/mantaflow/helper/util/interpol.h b/extern/mantaflow/helper/util/interpol.h
new file mode 100644
index 00000000000..24d5d2ada06
--- /dev/null
+++ b/extern/mantaflow/helper/util/interpol.h
@@ -0,0 +1,324 @@
+/******************************************************************************
+ *
+ * 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 interpolation
+ *
+ ******************************************************************************/
+
+#ifndef _INTERPOL_H
+#define _INTERPOL_H
+
+#include "vectorbase.h"
+
+// Grid values are stored at i+0.5, j+0.5, k+0.5
+// MAC grid values are stored at i,j+0.5,k+0.5 (for x) ...
+
+namespace Manta {
+
+inline Vec3 fdTangent(const Vec3 &p0, const Vec3 &p1, const Vec3 &p2)
+{
+  return 0.5 * (getNormalized(p2 - p1) + getNormalized(p1 - p0));
+}
+
+inline Vec3 crTangent(const Vec3 &p0, const Vec3 &p1, const Vec3 &p2)
+{
+  return 0.5 * (p2 - p0);
+}
+
+inline Vec3 hermiteSpline(const Vec3 &p0, const Vec3 &p1, const Vec3 &m0, const Vec3 &m1, Real t)
+{
+  const Real t2 = t * t, t3 = t2 * t;
+  return (2.0 * t3 - 3.0 * t2 + 1.0) * p0 + (t3 - 2.0 * t2 + t) * m0 +
+         (-2.0 * t3 + 3.0 * t2) * p1 + (t3 - t2) * m1;
+}
+
+static inline void checkIndexInterpol(const Vec3i &size, IndexInt idx)
+{
+  if (idx < 0 || idx > (IndexInt)size.x * size.y * size.z) {
+    std::ostringstream s;
+    s << "Grid interpol dim " << size << " : index " << idx << " out of bound ";
+    errMsg(s.str());
+  }
+}
+
+// ----------------------------------------------------------------------
+// Grid interpolators
+// ----------------------------------------------------------------------
+
+#define BUILD_INDEX \
+  Real px = pos.x - 0.5f, py = pos.y - 0.5f, pz = pos.z - 0.5f; \
+  int xi = (int)px; \
+  int yi = (int)py; \
+  int zi = (int)pz; \
+  Real s1 = px - (Real)xi, s0 = 1. - s1; \
+  Real t1 = py - (Real)yi, t0 = 1. - t1; \
+  Real f1 = pz - (Real)zi, f0 = 1. - f1; \
+  /* clamp to border */ \
+  if (px < 0.) { \
+    xi = 0; \
+    s0 = 1.0; \
+    s1 = 0.0; \
+  } \
+  if (py < 0.) { \
+    yi = 0; \
+    t0 = 1.0; \
+    t1 = 0.0; \
+  } \
+  if (pz < 0.) { \
+    zi = 0; \
+    f0 = 1.0; \
+    f1 = 0.0; \
+  } \
+  if (xi >= size.x - 1) { \
+    xi = size.x - 2; \
+    s0 = 0.0; \
+    s1 = 1.0; \
+  } \
+  if (yi >= size.y - 1) { \
+    yi = size.y - 2; \
+    t0 = 0.0; \
+    t1 = 1.0; \
+  } \
+  if (size.z > 1) { \
+    if (zi >= size.z - 1) { \
+      zi = size.z - 2; \
+      f0 = 0.0; \
+      f1 = 1.0; \
+    } \
+  } \
+  const int X = 1; \
+  const int Y = size.x;
+
+template<class T> inline T interpol(const T *data, const Vec3i &size, const int Z, const Vec3 &pos)
+{
+  BUILD_INDEX
+  IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi;
+  DEBUG_ONLY(checkIndexInterpol(size, idx));
+  DEBUG_ONLY(checkIndexInterpol(size, idx + X + Y + Z));
+
+  return ((data[idx] * t0 + data[idx + Y] * t1) * s0 +
+          (data[idx + X] * t0 + data[idx + X + Y] * t1) * s1) *
+             f0 +
+         ((data[idx + Z] * t0 + data[idx + Y + Z] * t1) * s0 +
+          (data[idx + X + Z] * t0 + data[idx + X + Y + Z] * t1) * s1) *
+             f1;
+}
+
+template<int c>
+inline Real interpolComponent(const Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos)
+{
+  BUILD_INDEX
+  IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi;
+  DEBUG_ONLY(checkIndexInterpol(size, idx));
+  DEBUG_ONLY(checkIndexInterpol(size, idx + X + Y + Z));
+
+  return ((data[idx][c] * t0 + data[idx + Y][c] * t1) * s0 +
+          (data[idx + X][c] * t0 + data[idx + X + Y][c] * t1) * s1) *
+             f0 +
+         ((data[idx + Z][c] * t0 + data[idx + Y + Z][c] * t1) * s0 +
+          (data[idx + X + Z][c] * t0 + data[idx + X + Y + Z][c] * t1) * s1) *
+             f1;
+}
+
+template<class T>
+inline void setInterpol(
+    T *data, const Vec3i &size, const int Z, const Vec3 &pos, const T &v, Real *sumBuffer)
+{
+  BUILD_INDEX
+  IndexInt idx = (IndexInt)xi + (IndexInt)Y * yi + (IndexInt)Z * zi;
+  DEBUG_ONLY(checkIndexInterpol(size, idx));
+  DEBUG_ONLY(checkIndexInterpol(size, idx + X + Y + Z));
+
+  T *ref = &data[idx];
+  Real *sum = &sumBuffer[idx];
+  Real s0f0 = s0 * f0, s1f0 = s1 * f0, s0f1 = s0 * f1, s1f1 = s1 * f1;
+  Real w0 = t0 * s0f0, wx = t0 * s1f0, wy = t1 * s0f0, wxy = t1 * s1f0;
+  Real wz = t0 * s0f1, wxz = t0 * s1f1, wyz = t1 * s0f1, wxyz = t1 * s1f1;
+
+  sum[Z] += wz;
+  sum[X + Z] += wxz;
+  sum[Y + Z] += wyz;
+  sum[X + Y + Z] += wxyz;
+  ref[Z] += wz * v;
+  ref[X + Z] += wxz * v;
+  ref[Y + Z] += wyz * v;
+  ref[X + Y + Z] += wxyz * v;
+  sum[0] += w0;
+  sum[X] += wx;
+  sum[Y] += wy;
+  sum[X + Y] += wxy;
+  ref[0] += w0 * v;
+  ref[X] += wx * v;
+  ref[Y] += wy * v;
+  ref[X + Y] += wxy * v;
+}
+
+#define BUILD_INDEX_SHIFT \
+  BUILD_INDEX \
+  /* shifted coords */ \
+  int s_xi = (int)pos.x, s_yi = (int)pos.y, s_zi = (int)pos.z; \
+  Real s_s1 = pos.x - (Real)s_xi, s_s0 = 1. - s_s1; \
+  Real s_t1 = pos.y - (Real)s_yi, s_t0 = 1. - s_t1; \
+  Real s_f1 = pos.z - (Real)s_zi, s_f0 = 1. - s_f1; \
+  /* clamp to border */ \
+  if (pos.x < 0) { \
+    s_xi = 0; \
+    s_s0 = 1.0; \
+    s_s1 = 0.0; \
+  } \
+  if (pos.y < 0) { \
+    s_yi = 0; \
+    s_t0 = 1.0; \
+    s_t1 = 0.0; \
+  } \
+  if (pos.z < 0) { \
+    s_zi = 0; \
+    s_f0 = 1.0; \
+    s_f1 = 0.0; \
+  } \
+  if (s_xi >= size.x - 1) { \
+    s_xi = size.x - 2; \
+    s_s0 = 0.0; \
+    s_s1 = 1.0; \
+  } \
+  if (s_yi >= size.y - 1) { \
+    s_yi = size.y - 2; \
+    s_t0 = 0.0; \
+    s_t1 = 1.0; \
+  } \
+  if (size.z > 1) { \
+    if (s_zi >= size.z - 1) { \
+      s_zi = size.z - 2; \
+      s_f0 = 0.0; \
+      s_f1 = 1.0; \
+    } \
+  }
+
+inline Vec3 interpolMAC(const Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos)
+{
+  BUILD_INDEX_SHIFT;
+  DEBUG_ONLY(checkIndexInterpol(size, (zi * (IndexInt)size.y + yi) * (IndexInt)size.x + xi));
+  DEBUG_ONLY(checkIndexInterpol(
+      size, (s_zi * (IndexInt)size.y + s_yi) * (IndexInt)size.x + s_xi + X + Y + Z));
+
+  // process individual components
+  Vec3 ret(0.);
+  {  // X
+    const Vec3 *ref = &data[((zi * size.y + yi) * size.x + s_xi)];
+    ret.x = f0 * ((ref[0].x * t0 + ref[Y].x * t1) * s_s0 +
+                  (ref[X].x * t0 + ref[X + Y].x * t1) * s_s1) +
+            f1 * ((ref[Z].x * t0 + ref[Z + Y].x * t1) * s_s0 +
+                  (ref[X + Z].x * t0 + ref[X + Y + Z].x * t1) * s_s1);
+  }
+  {  // Y
+    const Vec3 *ref = &data[((zi * size.y + s_yi) * size.x + xi)];
+    ret.y = f0 * ((ref[0].y * s_t0 + ref[Y].y * s_t1) * s0 +
+                  (ref[X].y * s_t0 + ref[X + Y].y * s_t1) * s1) +
+            f1 * ((ref[Z].y * s_t0 + ref[Z + Y].y * s_t1) * s0 +
+                  (ref[X + Z].y * s_t0 + ref[X + Y + Z].y * s_t1) * s1);
+  }
+  {  // Z
+    const Vec3 *ref = &data[((s_zi * size.y + yi) * size.x + xi)];
+    ret.z = s_f0 *
+                ((ref[0].z * t0 + ref[Y].z * t1) * s0 + (ref[X].z * t0 + ref[X + Y].z * t1) * s1) +
+            s_f1 * ((ref[Z].z * t0 + ref[Z + Y].z * t1) * s0 +
+                    (ref[X + Z].z * t0 + ref[X + Y + Z].z * t1) * s1);
+  }
+  return ret;
+}
+
+inline void setInterpolMAC(
+    Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos, const Vec3 &val, Vec3 *sumBuffer)
+{
+  BUILD_INDEX_SHIFT;
+  DEBUG_ONLY(checkIndexInterpol(size, (zi * (IndexInt)size.y + yi) * (IndexInt)size.x + xi));
+  DEBUG_ONLY(checkIndexInterpol(
+      size, (s_zi * (IndexInt)size.y + s_yi) * (IndexInt)size.x + s_xi + X + Y + Z));
+
+  // process individual components
+  {  // X
+    const IndexInt idx = (IndexInt)(zi * size.y + yi) * size.x + s_xi;
+    Vec3 *ref = &data[idx], *sum = &sumBuffer[idx];
+    Real s0f0 = s_s0 * f0, s1f0 = s_s1 * f0, s0f1 = s_s0 * f1, s1f1 = s_s1 * f1;
+    Real w0 = t0 * s0f0, wx = t0 * s1f0, wy = t1 * s0f0, wxy = t1 * s1f0;
+    Real wz = t0 * s0f1, wxz = t0 * s1f1, wyz = t1 * s0f1, wxyz = t1 * s1f1;
+
+    sum[Z].x += wz;
+    sum[X + Z].x += wxz;
+    sum[Y + Z].x += wyz;
+    sum[X + Y + Z].x += wxyz;
+    ref[Z].x += wz * val.x;
+    ref[X + Z].x += wxz * val.x;
+    ref[Y + Z].x += wyz * val.x;
+    ref[X + Y + Z].x += wxyz * val.x;
+    sum[0].x += w0;
+    sum[X].x += wx;
+    sum[Y].x += wy;
+    sum[X + Y].x += wxy;
+    ref[0].x += w0 * val.x;
+    ref[X].x += wx * val.x;
+    ref[Y].x += wy * val.x;
+    ref[X + Y].x += wxy * val.x;
+  }
+  {  // Y
+    const IndexInt idx = (IndexInt)(zi * size.y + s_yi) * size.x + xi;
+    Vec3 *ref = &data[idx], *sum = &sumBuffer[idx];
+    Real s0f0 = s0 * f0, s1f0 = s1 * f0, s0f1 = s0 * f1, s1f1 = s1 * f1;
+    Real w0 = s_t0 * s0f0, wx = s_t0 * s1f0, wy = s_t1 * s0f0, wxy = s_t1 * s1f0;
+    Real wz = s_t0 * s0f1, wxz = s_t0 * s1f1, wyz = s_t1 * s0f1, wxyz = s_t1 * s1f1;
+
+    sum[Z].y += wz;
+    sum[X + Z].y += wxz;
+    sum[Y + Z].y += wyz;
+    sum[X + Y + Z].y += wxyz;
+    ref[Z].y += wz * val.y;
+    ref[X + Z].y += wxz * val.y;
+    ref[Y + Z].y += wyz * val.y;
+    ref[X + Y + Z].y += wxyz * val.y;
+    sum[0].y += w0;
+    sum[X].y += wx;
+    sum[Y].y += wy;
+    sum[X + Y].y += wxy;
+    ref[0].y += w0 * val.y;
+    ref[X].y += wx * val.y;
+    ref[Y].y += wy * val.y;
+    ref[X + Y].y += wxy * val.y;
+  }
+  {  // Z
+    const IndexInt idx = (IndexInt)(s_zi * size.y + yi) * size.x + xi;
+    Vec3 *ref = &data[idx], *sum = &sumBuffer[idx];
+    Real s0f0 = s0 * s_f0, s1f0 = s1 * s_f0, s0f1 = s0 * s_f1, s1f1 = s1 * s_f1;
+    Real w0 = t0 * s0f0, wx = t0 * s1f0, wy = t1 * s0f0, wxy = t1 * s1f0;
+    Real wz = t0 * s0f1, wxz = t0 * s1f1, wyz = t1 * s0f1, wxyz = t1 * s1f1;
+
+    sum[0].z += w0;
+    sum[X].z += wx;
+    sum[Y].z += wy;
+    sum[X + Y].z += wxy;
+    sum[Z].z += wz;
+    sum[X + Z].z += wxz;
+    sum[Y + Z].z += wyz;
+    sum[X + Y + Z].z += wxyz;
+    ref[0].z += w0 * val.z;
+    ref[X].z += wx * val.z;
+    ref[Y].z += wy * val.z;
+    ref[X + Y].z += wxy * val.z;
+    ref[Z].z += wz * val.z;
+    ref[X + Z].z += wxz * val.z;
+    ref[Y + Z].z += wyz * val.z;
+    ref[X + Y + Z].z += wxyz * val.z;
+  }
+}
+
+#undef BUILD_INDEX
+#undef BUILD_INDEX_SHIFT
+
+}  // namespace Manta
+
+#endif
diff --git a/extern/mantaflow/helper/util/interpolHigh.h b/extern/mantaflow/helper/util/interpolHigh.h
new file mode 100644
index 00000000000..c2a77442b9c
--- /dev/null
+++ b/extern/mantaflow/helper/util/interpolHigh.h
@@ -0,0 +1,204 @@
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2014 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 higher order interpolation
+ *
+ ******************************************************************************/
+
+#ifndef _INTERPOLHIGH_H
+#define _INTERPOLHIGH_H
+
+#include "vectorbase.h"
+
+namespace Manta {
+
+template<class T> inline T cubicInterp(const Real interp, const T *points)
+{
+  T d0 = (points[2] - points[0]) * 0.5;
+  T d1 = (points[3] - points[1]) * 0.5;
+  T deltak = (points[2] - points[1]);
+
+  // disabled: if (deltak * d0 < 0.0) d0 = 0;
+  // disabled: if (deltak * d1 < 0.0) d1 = 0;
+
+  T a0 = points[1];
+  T a1 = d0;
+  T a2 = 3.0 * deltak - 2.0 * d0 - d1;
+  T a3 = -2.0 * deltak + d0 + d1;
+
+  Real squared = interp * interp;
+  Real cubed = squared * interp;
+  return a3 * cubed + a2 * squared + a1 * interp + a0;
+}
+
+template<class T> inline T interpolCubic2D(const T *data, const Vec3i &size, const Vec3 &pos)
+{
+  const Real px = pos.x - 0.5f, py = pos.y - 0.5f;
+
+  const int x1 = (int)px;
+  const int x2 = x1 + 1;
+  const int x3 = x1 + 2;
+  const int x0 = x1 - 1;
+
+  const int y1 = (int)py;
+  const int y2 = y1 + 1;
+  const int y3 = y1 + 2;
+  const int y0 = y1 - 1;
+
+  if (x0 < 0 || y0 < 0 || x3 >= size[0] || y3 >= size[1]) {
+    return interpol(data, size, 0, pos);
+  }
+
+  const Real xInterp = px - x1;
+  const Real yInterp = py - y1;
+
+  const int y0x = y0 * size[0];
+  const int y1x = y1 * size[0];
+  const int y2x = y2 * size[0];
+  const int y3x = y3 * size[0];
+
+  const T p0[] = {data[x0 + y0x], data[x1 + y0x], data[x2 + y0x], data[x3 + y0x]};
+  const T p1[] = {data[x0 + y1x], data[x1 + y1x], data[x2 + y1x], data[x3 + y1x]};
+  const T p2[] = {data[x0 + y2x], data[x1 + y2x], data[x2 + y2x], data[x3 + y2x]};
+  const T p3[] = {data[x0 + y3x], data[x1 + y3x], data[x2 + y3x], data[x3 + y3x]};
+
+  const T finalPoints[] = {cubicInterp(xInterp, p0),
+                           cubicInterp(xInterp, p1),
+                           cubicInterp(xInterp, p2),
+                           cubicInterp(xInterp, p3)};
+
+  return cubicInterp(yInterp, finalPoints);
+}
+
+template<class T>
+inline T interpolCubic(const T *data, const Vec3i &size, const int Z, const Vec3 &pos)
+{
+  if (Z == 0)
+    return interpolCubic2D(data, size, pos);
+
+  const Real px = pos.x - 0.5f, py = pos.y - 0.5f, pz = pos.z - 0.5f;
+
+  const int x1 = (int)px;
+  const int x2 = x1 + 1;
+  const int x3 = x1 + 2;
+  const int x0 = x1 - 1;
+
+  const int y1 = (int)py;
+  const int y2 = y1 + 1;
+  const int y3 = y1 + 2;
+  const int y0 = y1 - 1;
+
+  const int z1 = (int)pz;
+  const int z2 = z1 + 1;
+  const int z3 = z1 + 2;
+  const int z0 = z1 - 1;
+
+  if (x0 < 0 || y0 < 0 || z0 < 0 || x3 >= size[0] || y3 >= size[1] || z3 >= size[2]) {
+    return interpol(data, size, Z, pos);
+  }
+
+  const Real xInterp = px - x1;
+  const Real yInterp = py - y1;
+  const Real zInterp = pz - z1;
+
+  const int slabsize = size[0] * size[1];
+  const int z0Slab = z0 * slabsize;
+  const int z1Slab = z1 * slabsize;
+  const int z2Slab = z2 * slabsize;
+  const int z3Slab = z3 * slabsize;
+
+  const int y0x = y0 * size[0];
+  const int y1x = y1 * size[0];
+  const int y2x = y2 * size[0];
+  const int y3x = y3 * size[0];
+
+  const int y0z0 = y0x + z0Slab;
+  const int y1z0 = y1x + z0Slab;
+  const int y2z0 = y2x + z0Slab;
+  const int y3z0 = y3x + z0Slab;
+
+  const int y0z1 = y0x + z1Slab;
+  const int y1z1 = y1x + z1Slab;
+  const int y2z1 = y2x + z1Slab;
+  const int y3z1 = y3x + z1Slab;
+
+  const int y0z2 = y0x + z2Slab;
+  const int y1z2 = y1x + z2Slab;
+  const int y2z2 = y2x + z2Slab;
+  const int y3z2 = y3x + z2Slab;
+
+  const int y0z3 = y0x + z3Slab;
+  const int y1z3 = y1x + z3Slab;
+  const int y2z3 = y2x + z3Slab;
+  const int y3z3 = y3x + z3Slab;
+
+  // get the z0 slice
+  const T p0[] = {data[x0 + y0z0], data[x1 + y0z0], data[x2 + y0z0], data[x3 + y0z0]};
+  const T p1[] = {data[x0 + y1z0], data[x1 + y1z0], data[x2 + y1z0], data[x3 + y1z0]};
+  const T p2[] = {data[x0 + y2z0], data[x1 + y2z0], data[x2 + y2z0], data[x3 + y2z0]};
+  const T p3[] = {data[x0 + y3z0], data[x1 + y3z0], data[x2 + y3z0], data[x3 + y3z0]};
+
+  // get the z1 slice
+  const T p4[] = {data[x0 + y0z1], data[x1 + y0z1], data[x2 + y0z1], data[x3 + y0z1]};
+  const T p5[] = {data[x0 + y1z1], data[x1 + y1z1], data[x2 + y1z1], data[x3 + y1z1]};
+  const T p6[] = {data[x0 + y2z1], data[x1 + y2z1], data[x2 + y2z1], data[x3 + y2z1]};
+  const T p7[] = {data[x0 + y3z1], data[x1 + y3z1], data[x2 + y3z1], data[x3 + y3z1]};
+
+  // get the z2 slice
+  const T p8[] = {data[x0 + y0z2], data[x1 + y0z2], data[x2 + y0z2], data[x3 + y0z2]};
+  const T p9[] = {data[x0 + y1z2], data[x1 + y1z2], data[x2 + y1z2], data[x3 + y1z2]};
+  const T p10[] = {data[x0 + y2z2], data[x1 + y2z2], data[x2 + y2z2], data[x3 + y2z2]};
+  const T p11[] = {data[x0 + y3z2], data[x1 + y3z2], data[x2 + y3z2], data[x3 + y3z2]};
+
+  // get the z3 slice
+  const T p12[] = {data[x0 + y0z3], data[x1 + y0z3], data[x2 + y0z3], data[x3 + y0z3]};
+  const T p13[] = {data[x0 + y1z3], data[x1 + y1z3], data[x2 + y1z3], data[x3 + y1z3]};
+  const T p14[] = {data[x0 + y2z3], data[x1 + y2z3], data[x2 + y2z3], data[x3 + y2z3]};
+  const T p15[] = {data[x0 + y3z3], data[x1 + y3z3], data[x2 + y3z3], data[x3 + y3z3]};
+
+  // interpolate
+  const T z0Points[] = {cubicInterp(xInterp, p0),
+                        cubicInterp(xInterp, p1),
+                        cubicInterp(xInterp, p2),
+                        cubicInterp(xInterp, p3)};
+  const T z1Points[] = {cubicInterp(xInterp, p4),
+                        cubicInterp(xInterp, p5),
+                        cubicInterp(xInterp, p6),
+                        cubicInterp(xInterp, p7)};
+  const T z2Points[] = {cubicInterp(xInterp, p8),
+                        cubicInterp(xInterp, p9),
+                        cubicInterp(xInterp, p10),
+                        cubicInterp(xInterp, p11)};
+  const T z3Points[] = {cubicInterp(xInterp, p12),
+                        cubicInterp(xInterp, p13),
+                        cubicInterp(xInterp, p14),
+                        cubicInterp(xInterp, p15)};
+
+  const T finalPoints[] = {cubicInterp(yInterp, z0Points),
+                           cubicInterp(yInterp, z1Points),
+                           cubicInterp(yInterp, z2Points),
+                           cubicInterp(yInterp, z3Points)};
+
+  return cubicInterp(zInterp, finalPoints);
+}
+
+inline Vec3 interpolCubicMAC(const Vec3 *data, const Vec3i &size, const int Z, const Vec3 &pos)
+{
+  // warning - not yet optimized...
+  Real vx = interpolCubic<Vec3>(data, size, Z, pos + Vec3(0.5, 0, 0))[0];
+  Real vy = interpolCubic<Vec3>(data, size, Z, pos + Vec3(0, 0.5, 0))[1];
+  Real vz = 0.f;
+  if (Z != 0)
+    vz = interpolCubic<Vec3>(data, size, Z, pos + Vec3(0, 0, 0.5))[2];
+  return Vec3(vx, vy, vz);
+}
+
+}  // namespace Manta
+
+#endif
diff --git a/extern/mantaflow/helper/util/matrixbase.h b/extern/mantaflow/helper/util/matrixbase.h
new file mode 100644
index 00000000000..9875998f0be
--- /dev/null
+++ b/extern/mantaflow/helper/util/matrixbase.h
@@ -0,0 +1,394 @@
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2015 Kiwon Um, Nils Thuerey
+ *
+ * This program is free software, distributed under the terms of the
+ * GNU General Public License (GPL)
+ * http://www.gnu.org/licenses
+ *
+ * Matrix (3x3) class
+ *
+ ******************************************************************************/
+
+#ifndef MATRIXBASE_H
+#define MATRIXBASE_H
+
+#include "vectorbase.h"
+
+namespace Manta {
+
+template<typename T> class Matrix3x3 {
+ public:
+  // NOTE: default is the identity matrix!
+  explicit Matrix3x3(const T &p00 = 1,
+                     const T &p01 = 0,
+                     const T &p02 = 0,
+                     const T &p10 = 0,
+                     const T &p11 = 1,
+                     const T &p12 = 0,
+                     const T &p20 = 0,
+                     const T &p21 = 0,
+                     const T &p22 = 1)
+  {
+    v[0][0] = p00;
+    v[0][1] = p01;
+    v[0][2] = p02;
+    v[1][0] = p10;
+    v[1][1] = p11;
+    v[1][2] = p12;
+    v[2][0] = p20;
+    v[2][1] = p21;
+    v[2][2] = p22;
+  }
+
+  explicit Matrix3x3(const Vector3D<T> &diag)
+  {
+    v[0][0] = diag.x;
+    v[0][1] = 0;
+    v[0][2] = 0;
+    v[1][0] = 0;
+    v[1][1] = diag.y;
+    v[1][2] = 0;
+    v[2][0] = 0;
+    v[2][1] = 0;
+    v[2][2] = diag.z;
+  }
+
+  Matrix3x3(const Vector3D<T> &c0, const Vector3D<T> &c1, const Vector3D<T> &c2)
+  {
+    v[0][0] = c0.x;
+    v[0][1] = c1.x;
+    v[0][2] = c2.x;
+    v[1][0] = c0.y;
+    v[1][1] = c1.y;
+    v[1][2] = c2.y;
+    v[2][0] = c0.z;
+    v[2][1] = c1.z;
+    v[2][2] = c2.z;
+  }
+
+  // assignment operators
+  Matrix3x3 &operator+=(const Matrix3x3 &m)
+  {
+    v00 += m.v00;
+    v01 += m.v01;
+    v02 += m.v02;
+    v10 += m.v10;
+    v11 += m.v11;
+    v12 += m.v12;
+    v20 += m.v20;
+    v21 += m.v21;
+    v22 += m.v22;
+    return *this;
+  }
+  Matrix3x3 &operator-=(const Matrix3x3 &m)
+  {
+    v00 -= m.v00;
+    v01 -= m.v01;
+    v02 -= m.v02;
+    v10 -= m.v10;
+    v11 -= m.v11;
+    v12 -= m.v12;
+    v20 -= m.v20;
+    v21 -= m.v21;
+    v22 -= m.v22;
+    return *this;
+  }
+  Matrix3x3 &operator*=(const T s)
+  {
+    v00 *= s;
+    v01 *= s;
+    v02 *= s;
+    v10 *= s;
+    v11 *= s;
+    v12 *= s;
+    v20 *= s;
+    v21 *= s;
+    v22 *= s;
+    return *this;
+  }
+  Matrix3x3 &operator/=(const T s)
+  {
+    v00 /= s;
+    v01 /= s;
+    v02 /= s;
+    v10 /= s;
+    v11 /= s;
+    v12 /= s;
+    v20 /= s;
+    v21 /= s;
+    v22 /= s;
+    return *this;
+  }
+
+  // binary operators
+  Matrix3x3 operator+(const Matrix3x3 &m) const
+  {
+    return Matrix3x3(*this) += m;
+  }
+  Matrix3x3 operator-(const Matrix3x3 &m) const
+  {
+    return Matrix3x3(*this) -= m;
+  }
+  Matrix3x3 operator*(const Matrix3x3 &m) const
+  {
+    return Matrix3x3(v00 * m.v00 + v01 * m.v10 + v02 * m.v20,
+                     v00 * m.v01 + v01 * m.v11 + v02 * m.v21,
+                     v00 * m.v02 + v01 * m.v12 + v02 * m.v22,
+
+                     v10 * m.v00 + v11 * m.v10 + v12 * m.v20,
+                     v10 * m.v01 + v11 * m.v11 + v12 * m.v21,
+                     v10 * m.v02 + v11 * m.v12 + v12 * m.v22,
+
+                     v20 * m.v00 + v21 * m.v10 + v22 * m.v20,
+                     v20 * m.v01 + v21 * m.v11 + v22 * m.v21,
+                     v20 * m.v02 + v21 * m.v12 + v22 * m.v22);
+  }
+  Matrix3x3 operator*(const T s) const
+  {
+    return Matrix3x3(*this) *= s;
+  }
+  Vector3D<T> operator*(const Vector3D<T> &v) const
+  {
+    return Vector3D<T>(v00 * v.x + v01 * v.y + v02 * v.z,
+                       v10 * v.x + v11 * v.y + v12 * v.z,
+                       v20 * v.x + v21 * v.y + v22 * v.z);
+  }
+  Vector3D<T> transposedMul(const Vector3D<T> &v) const
+  {
+    // M^T*v
+    return Vector3D<T>(v00 * v.x + v10 * v.y + v20 * v.z,
+                       v01 * v.x + v11 * v.y + v21 * v.z,
+                       v02 * v.x + v12 * v.y + v22 * v.z);
+  }
+  Matrix3x3 transposedMul(const Matrix3x3 &m) const
+  {
+    // M^T*M
+    return Matrix3x3(v00 * m.v00 + v10 * m.v10 + v20 * m.v20,
+                     v00 * m.v01 + v10 * m.v11 + v20 * m.v21,
+                     v00 * m.v02 + v10 * m.v12 + v20 * m.v22,
+
+                     v01 * m.v00 + v11 * m.v10 + v21 * m.v20,
+                     v01 * m.v01 + v11 * m.v11 + v21 * m.v21,
+                     v01 * m.v02 + v11 * m.v12 + v21 * m.v22,
+
+                     v02 * m.v00 + v12 * m.v10 + v22 * m.v20,
+                     v02 * m.v01 + v12 * m.v11 + v22 * m.v21,
+                     v02 * m.v02 + v12 * m.v12 + v22 * m.v22);
+  }
+  Matrix3x3 mulTranspose(const Matrix3x3 &m) const
+  {
+    // M*m^T
+    return Matrix3x3(v00 * m.v00 + v01 * m.v01 + v02 * m.v02,
+                     v00 * m.v10 + v01 * m.v11 + v02 * m.v12,
+                     v00 * m.v20 + v01 * m.v21 + v02 * m.v22,
+
+                     v10 * m.v00 + v11 * m.v01 + v12 * m.v02,
+                     v10 * m.v10 + v11 * m.v11 + v12 * m.v12,
+                     v10 * m.v20 + v11 * m.v21 + v12 * m.v22,
+
+                     v20 * m.v00 + v21 * m.v01 + v22 * m.v02,
+                     v20 * m.v10 + v21 * m.v11 + v22 * m.v12,
+                     v20 * m.v20 + v21 * m.v21 + v22 * m.v22);
+  }
+
+  bool operator==(const Matrix3x3 &m) const
+  {
+    return (v00 == m.v00 && v01 == m.v01 && v02 == m.v02 && v10 == m.v10 && v11 == m.v11 &&
+            v12 == m.v12 && v20 == m.v20 && v21 == m.v21 && v22 == m.v22);
+  }
+
+  const T &operator()(const int r, const int c) const
+  {
+    return v[r][c];
+  }
+  T &operator()(const int r, const int c)
+  {
+    return const_cast<T &>(const_cast<const Matrix3x3 &>(*this)(r, c));
+  }
+
+  T trace() const
+  {
+    return v00 + v11 + v22;
+  }
+  T sumSqr() const
+  {
+    return (v00 * v00 + v01 * v01 + v02 * v02 + v10 * v10 + v11 * v11 + v12 * v12 + v20 * v20 +
+            v21 * v21 + v22 * v22);
+  }
+
+  Real determinant() const
+  {
+    return (v00 * v11 * v22 - v00 * v12 * v21 + v01 * v12 * v20 - v01 * v10 * v22 +
+            v02 * v10 * v21 - v02 * v11 * v20);
+  }
+  Matrix3x3 &transpose()
+  {
+    return *this = transposed();
+  }
+  Matrix3x3 transposed() const
+  {
+    return Matrix3x3(v00, v10, v20, v01, v11, v21, v02, v12, v22);
+  }
+  Matrix3x3 &invert()
+  {
+    return *this = inverse();
+  }
+  Matrix3x3 inverse() const
+  {
+    const Real det = determinant();  // FIXME: assert(det);
+    const Real idet = 1e0 / det;
+    return Matrix3x3(idet * (v11 * v22 - v12 * v21),
+                     idet * (v02 * v21 - v01 * v22),
+                     idet * (v01 * v12 - v02 * v11),
+                     idet * (v12 * v20 - v10 * v22),
+                     idet * (v00 * v22 - v02 * v20),
+                     idet * (v02 * v10 - v00 * v12),
+                     idet * (v10 * v21 - v11 * v20),
+                     idet * (v01 * v20 - v00 * v21),
+                     idet * (v00 * v11 - v01 * v10));
+  }
+  bool getInverse(Matrix3x3 &inv) const
+  {
+    const Real det = determinant();
+    if (det == 0e0)
+      return false;  // FIXME: is it likely to happen the floating error?
+
+    const Real idet = 1e0 / det;
+    inv.v00 = idet * (v11 * v22 - v12 * v21);
+    inv.v01 = idet * (v02 * v21 - v01 * v22);
+    inv.v02 = idet * (v01 * v12 - v02 * v11);
+
+    inv.v10 = idet * (v12 * v20 - v10 * v22);
+    inv.v11 = idet * (v00 * v22 - v02 * v20);
+    inv.v12 = idet * (v02 * v10 - v00 * v12);
+
+    inv.v20 = idet * (v10 * v21 - v11 * v20);
+    inv.v21 = idet * (v01 * v20 - v00 * v21);
+    inv.v22 = idet * (v00 * v11 - v01 * v10);
+
+    return true;
+  }
+
+  Real normOne() const
+  {
+    // the maximum absolute column sum of the matrix
+    return max(std::fabs(v00) + std::fabs(v10) + std::fabs(v20),
+               std::fabs(v01) + std::fabs(v11) + std::fabs(v21),
+               std::fabs(v02) + std::fabs(v12) + std::fabs(v22));
+  }
+  Real normInf() const
+  {
+    // the maximum absolute row sum of the matrix
+    return max(std::fabs(v00) + std::fabs(v01) + std::fabs(v02),
+               std::fabs(v10) + std::fabs(v11) + std::fabs(v12),
+               std::fabs(v20) + std::fabs(v21) + std::fabs(v22));
+  }
+
+  Vector3D<T> eigenvalues() const
+  {
+    Vector3D<T> eigen;
+
+    const Real b = -v00 - v11 - v22;
+    const Real c = v00 * (v11 + v22) + v11 * v22 - v12 * v21 - v01 * v10 - v02 * v20;
+    Real d = -v00 * (v11 * v22 - v12 * v21) - v20 * (v01 * v12 - v11 * v02) -
+             v10 * (v02 * v21 - v22 * v01);
+    const Real f = (3.0 * c - b * b) / 3.0;
+    const Real g = (2.0 * b * b * b - 9.0 * b * c + 27.0 * d) / 27.0;
+    const Real h = g * g / 4.0 + f * f * f / 27.0;
+
+    Real sign;
+    if (h > 0) {
+      Real r = -g / 2.0 + std::sqrt(h);
+      if (r < 0) {
+        r = -r;
+        sign = -1.0;
+      }
+      else
+        sign = 1.0;
+      Real s = sign * std::pow(r, 1.0 / 3.0);
+      Real t = -g / 2.0 - std::sqrt(h);
+      if (t < 0) {
+        t = -t;
+        sign = -1.0;
+      }
+      else
+        sign = 1.0;
+      Real u = sign * std::pow(t, 1.0 / 3.0);
+      eigen[0] = (s + u) - b / 3.0;
+      eigen[1] = eigen[2] = 0;
+    }
+    else if (h == 0) {
+      if (d < 0) {
+        d = -d;
+        sign = -1.0;
+      }
+      sign = 1.0;
+      eigen[0] = -1.0 * sign * std::pow(d, 1.0 / 3.0);
+      eigen[1] = eigen[2] = 0;
+    }
+    else {
+      const Real i = std::sqrt(g * g / 4.0 - h);
+      const Real j = std::pow(i, 1.0 / 3.0);
+      const Real k = std::acos(-g / (2.0 * i));
+      const Real l = -j;
+      const Real m = std::cos(k / 3.0);
+      const Real n = std::sqrt(3.0) * std::sin(k / 3.0);
+      const Real p = -b / 3.0;
+      eigen[0] = 2e0 * j * m + p;
+      eigen[1] = l * (m + n) + p;
+      eigen[2] = l * (m - n) + p;
+    }
+
+    return eigen;
+  }
+
+  static Matrix3x3 I()
+  {
+    return Matrix3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
+  }
+
+#ifdef _WIN32
+#  pragma warning(disable : 4201)
+#endif
+  union {
+    struct {
+      T v00, v01, v02, v10, v11, v12, v20, v21, v22;
+    };
+    T v[3][3];
+    T v1[9];
+  };
+#ifdef _WIN32
+#  pragma warning(default : 4201)
+#endif
+};
+
+template<typename T1, typename T> inline Matrix3x3<T> operator*(const T1 s, const Matrix3x3<T> &m)
+{
+  return m * static_cast<T>(s);
+}
+
+template<typename T> inline Matrix3x3<T> crossProductMatrix(const Vector3D<T> &v)
+{
+  return Matrix3x3<T>(0, -v.z, v.y, v.z, 0, -v.x, -v.y, v.x, 0);
+}
+
+template<typename T> inline Matrix3x3<T> outerProduct(const Vector3D<T> &a, const Vector3D<T> &b)
+{
+  return Matrix3x3<T>(a.x * b.x,
+                      a.x * b.y,
+                      a.x * b.z,
+                      a.y * b.x,
+                      a.y * b.y,
+                      a.y * b.z,
+                      a.z * b.x,
+                      a.z * b.y,
+                      a.z * b.z);
+}
+
+typedef Matrix3x3<Real> Matrix3x3f;
+
+}  // namespace Manta
+
+#endif /* MATRIXBASE_H */
diff --git a/extern/mantaflow/helper/util/mcubes.h b/extern/mantaflow/helper/util/mcubes.h
new file mode 100644
index 00000000000..bd1c780e932
--- /dev/null
+++ b/extern/mantaflow/helper/util/mcubes.h
@@ -0,0 +1,308 @@
+/******************************************************************************
+ *
+ * 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
+ *
+ * Marching cubes lookup indices
+ *
+ ******************************************************************************/
+
+#ifndef _MCUBES_H_
+#define _MCUBES_H_
+
+static const int mcEdges[24] = {0, 1, 1, 2, 2, 3, 3, 0, 4, 5, 5, 6,
+                                6, 7, 7, 4, 0, 4, 1, 5, 2, 6, 3, 7};
+
+static const int cubieOffsetX[8] = {0, 1, 1, 0, 0, 1, 1, 0};
+static const int cubieOffsetY[8] = {0, 0, 1, 1, 0, 0, 1, 1};
+static const int cubieOffsetZ[8] = {0, 0, 0, 0, 1, 1, 1, 1};
+
+/* which edges are needed ? */
+/* cf. http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/ */
+static const short mcEdgeTable[256] = {
+    0x0,   0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a,
+    0xd03, 0xe09, 0xf00, 0x190, 0x99,  0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895,
+    0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, 0x230, 0x339, 0x33,  0x13a, 0x636, 0x73f, 0x435,
+    0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9, 0x1a3, 0xaa,
+    0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, 0x460,
+    0x569, 0x663, 0x76a, 0x66,  0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963,
+    0xa69, 0xb60, 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff,  0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff,
+    0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55,  0x15c,
+    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6,
+    0x2cf, 0x1c5, 0xcc,  0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9,
+    0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc,  0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9,
+    0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55,  0x35f, 0x256,
+    0x55a, 0x453, 0x759, 0x650, 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc,
+    0x3f5, 0xff,  0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f,
+    0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66,  0x76a, 0x663, 0x569, 0x460, 0xca0, 0xda9, 0xea3,
+    0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa,  0x1a3, 0x2a9, 0x3a0,
+    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a,
+    0x33,  0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795,
+    0x49f, 0x596, 0x29a, 0x393, 0x99,  0x190, 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905,
+    0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0};
+
+/* triangles for the 256 intersection possibilities */
+/* cf. http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/ */
+static const short mcTriTable[256][16] = {
+    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
+    {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
+    {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+    {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
+    {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
+    {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
+    {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
+    {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
+    {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
+    {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
+    {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
+    {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
+    {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
+    {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+    {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
+    {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
+    {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
+    {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
+    {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
+    {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
+    {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
+    {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
+    {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
+    {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
+    {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
+    {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
+    {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
+    {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
+    {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
+    {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
+    {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
+    {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
+    {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
+    {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
+    {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
+    {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
+    {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
+    {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
+    {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
+    {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
+    {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
+    {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
+    {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
+    {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
+    {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
+    {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
+    {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
+    {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
+    {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
+    {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
+    {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
+    {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
+    {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
+    {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
+    {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
+    {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
+    {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
+    {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
+    {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
+    {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
+    {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
+    {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
+    {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
+    {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
+    {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
+    {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
+    {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
+    {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
+    {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
+    {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
+    {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
+    {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
+    {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
+    {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
+    {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
+    {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
+    {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
+    {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
+    {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
+    {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
+    {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
+    {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
+    {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
+    {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
+    {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
+    {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
+    {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
+    {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
+    {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
+    {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
+    {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
+    {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
+    {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
+    {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
+    {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
+    {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
+    {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
+    {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
+    {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
+    {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
+    {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
+    {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
+    {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
+    {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
+    {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
+    {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
+    {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
+    {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
+    {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
+    {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
+    {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
+    {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
+    {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
+    {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
+    {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
+    {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
+    {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
+    {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
+    {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
+    {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
+    {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
+    {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
+    {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
+    {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
+    {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
+    {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
+    {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
+    {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
+    {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
+    {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
+    {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
+    {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
+    {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
+    {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
+    {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
+    {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
+    {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
+    {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
+    {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
+    {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
+    {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
+    {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
+    {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
+    {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
+    {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
+    {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
+    {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
+    {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
+    {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
+    {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
+    {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
+    {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
+    {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
+    {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
+    {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
+    {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
+    {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
+    {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
+    {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
+    {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+    {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
+    {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
+    {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
+    {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
+    {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
+    {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
+    {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
+    {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
+    {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
+    {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
+    {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
+    {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
+    {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
+    {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+    {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
+    {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
+    {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
+    {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
+    {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
+
+#endif
\ No newline at end of file
diff --git a/extern/mantaflow/helper/util/quaternion.h b/extern/mantaflow/helper/util/quaternion.h
new file mode 100644
index 00000000000..c4e161baee2
--- /dev/null
+++ b/extern/mantaflow/helper/util/quaternion.h
@@ -0,0 +1,103 @@
+/******************************************************************************
+ *
+ * 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
+ *
+ * Basic quaternion class
+ *
+ ******************************************************************************/
+
+#ifndef _QUATERNION_H
+#define _QUATERNION_H
+
+#include "vectorbase.h"
+
+namespace Manta {
+
+//! Very basic quaternion class
+class Quaternion {
+ public:
+  //! default constructor
+  Quaternion() : x(0), y(0), z(0), w(0)
+  {
+  }
+
+  //! copy constructor
+  Quaternion(const Quaternion &q) : x(q.x), y(q.y), z(q.z), w(q.w)
+  {
+  }
+
+  //! construct a quaternion from members
+  Quaternion(Real _x, Real _y, Real _z, Real _w) : x(_x), y(_y), z(_z), w(_w)
+  {
+  }
+
+  //! construct a quaternion from imag/real parts
+  Quaternion(Vec3 i, Real r) : x(i.x), y(i.y), z(i.z), w(r)
+  {
+  }
+
+  //! Assign operator
+  inline Quaternion &operator=(const Quaternion &q)
+  {
+    x = q.x;
+    y = q.y;
+    z = q.z;
+    w = q.w;
+    return *this;
+  }
+
+  //! Assign multiplication operator
+  inline Quaternion &operator*=(const Real a)
+  {
+    x *= a;
+    y *= a;
+    z *= a;
+    w *= a;
+    return *this;
+  }
+
+  //! return inverse quaternion
+  inline Quaternion inverse() const
+  {
+    Real mag = 1.0 / (x * x + y * y + z * z + w * w);
+    return Quaternion(-x * mag, -y * mag, -z * mag, w * mag);
+  }
+
+  //! imaginary part accessor
+  inline Vec3 imag()
+  {
+    return Vec3(x, y, z);
+  }
+
+  // imaginary part
+  Real x;
+  Real y;
+  Real z;
+
+  // real part
+  Real w;
+};
+
+//! Multiplication operator
+inline Quaternion operator*(const Quaternion &q1, const Quaternion &q2)
+{
+  return Quaternion(q2.w * q1.x + q2.x * q1.w + q2.y * q1.z - q2.z * q1.y,
+                    q2.w * q1.y + q2.y * q1.w + q2.z * q1.x - q2.x * q1.z,
+                    q2.w * q1.z + q2.z * q1.w + q2.x * q1.y - q2.y * q1.x,
+                    q2.w * q1.w - q2.x * q1.x - q2.y * q1.y - q2.z * q1.z);
+}
+
+//! Multiplication operator
+inline Quaternion operator*(const Quaternion &q, const Real a)
+{
+  return Quaternion(q.x * a, q.y * a, q.z * a, q.w * a);
+}
+
+}  // namespace Manta
+
+#endif
diff --git a/extern/mantaflow/helper/util/randomstream.h b/extern/mantaflow/helper/util/randomstream.h
new file mode 100644
index 00000000000..35b9c7d8858
--- /dev/null
+++ b/extern/mantaflow/helper/util/randomstream.h
@@ -0,0 +1,429 @@
+/******************************************************************************
+ *
+ * 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
+ *
+ * Random numbers
+ *
+ * Based on an example by Makoto Matsumoto, Takuji Nishimura, Shawn Cokus, and Richard J. Wagner
+ *
+ ******************************************************************************/
+
+#ifndef _RANDOMSTREAM_H
+#define _RANDOMSTREAM_H
+
+namespace Manta {
+
+#include <iostream>
+#include <stdio.h>
+#include <time.h>
+#include "vectorbase.h"
+
+class MTRand {
+  // Data
+ public:
+  typedef unsigned long uint32;  // unsigned integer type, at least 32 bits
+
+  enum { N = 624 };       // length of state vector
+  enum { SAVE = N + 1 };  // length of array for save()
+
+ protected:
+  enum { M = 397 };  // period parameter
+
+  uint32 state[N];  // internal state
+  uint32 *pNext;    // next value to get from state
+  int left;         // number of values left before reload needed
+
+  // Methods
+ public:
+  MTRand(const uint32 &oneSeed);                               // initialize with a simple uint32
+  MTRand(uint32 *const bigSeed, uint32 const seedLength = N);  // or an array
+  MTRand();  // auto-initialize with /dev/urandom or time() and clock()
+
+  // Do NOT use for CRYPTOGRAPHY without securely hashing several returned
+  // values together, otherwise the generator state can be learned after
+  // reading 624 consecutive values.
+
+  // Access to 32-bit random numbers
+  double rand();                       // real number in [0,1]
+  double rand(const double &n);        // real number in [0,n]
+  double randExc();                    // real number in [0,1)
+  double randExc(const double &n);     // real number in [0,n)
+  double randDblExc();                 // real number in (0,1)
+  double randDblExc(const double &n);  // real number in (0,n)
+  uint32 randInt();                    // integer in [0,2^32-1]
+  uint32 randInt(const uint32 &n);     // integer in [0,n] for n < 2^32
+  double operator()()
+  {
+    return rand();
+  }  // same as rand()
+
+  // Access to 53-bit random numbers (capacity of IEEE double precision)
+  double rand53();  // real number in [0,1)
+
+  // Access to nonuniform random number distributions
+  double randNorm(const double &mean = 0.0, const double &variance = 1.0);
+
+  // Re-seeding functions with same behavior as initializers
+  void seed(const uint32 oneSeed);
+  void seed(uint32 *const bigSeed, const uint32 seedLength = N);
+  void seed();
+
+  // Saving and loading generator state
+  void save(uint32 *saveArray) const;  // to array of size SAVE
+  void load(uint32 *const loadArray);  // from such array
+  friend std::ostream &operator<<(std::ostream &os, const MTRand &mtrand);
+  friend std::istream &operator>>(std::istream &is, MTRand &mtrand);
+
+ protected:
+  void initialize(const uint32 oneSeed);
+  void reload();
+  uint32 hiBit(const uint32 &u) const
+  {
+    return u & 0x80000000UL;
+  }
+  uint32 loBit(const uint32 &u) const
+  {
+    return u & 0x00000001UL;
+  }
+  uint32 loBits(const uint32 &u) const
+  {
+    return u & 0x7fffffffUL;
+  }
+  uint32 mixBits(const uint32 &u, const uint32 &v) const
+  {
+    return hiBit(u) | loBits(v);
+  }
+  uint32 twist(const uint32 &m, const uint32 &s0, const uint32 &s1) const
+  {
+    return m ^ (mixBits(s0, s1) >> 1) ^ (-loBit(s1) & 0x9908b0dfUL);
+  }
+  static uint32 hash(time_t t, clock_t c);
+};
+
+inline MTRand::MTRand(const uint32 &oneSeed)
+{
+  seed(oneSeed);
+}
+
+inline MTRand::MTRand(uint32 *const bigSeed, const uint32 seedLength)
+{
+  seed(bigSeed, seedLength);
+}
+
+inline MTRand::MTRand()
+{
+  seed();
+}
+
+inline double MTRand::rand()
+{
+  return double(randInt()) * (1.0 / 4294967295.0);
+}
+
+inline double MTRand::rand(const double &n)
+{
+  return rand() * n;
+}
+
+inline double MTRand::randExc()
+{
+  return double(randInt()) * (1.0 / 4294967296.0);
+}
+
+inline double MTRand::randExc(const double &n)
+{
+  return randExc() * n;
+}
+
+inline double MTRand::randDblExc()
+{
+  return (double(randInt()) + 0.5) * (1.0 / 4294967296.0);
+}
+
+inline double MTRand::randDblExc(const double &n)
+{
+  return randDblExc() * n;
+}
+
+inline double MTRand::rand53()
+{
+  uint32 a = randInt() >> 5, b = randInt() >> 6;
+  return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);  // by Isaku Wada
+}
+
+inline double MTRand::randNorm(const double &mean, const double &variance)
+{
+  // Return a real number from a normal (Gaussian) distribution with given
+  // mean and variance by Box-Muller method
+  double r = sqrt(-2.0 * log(1.0 - randDblExc())) * variance;
+  double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
+  return mean + r * cos(phi);
+}
+
+inline MTRand::uint32 MTRand::randInt()
+{
+  // Pull a 32-bit integer from the generator state
+  // Every other access function simply transforms the numbers extracted here
+
+  if (left == 0)
+    reload();
+  --left;
+
+  uint32 s1;
+  s1 = *pNext++;
+  s1 ^= (s1 >> 11);
+  s1 ^= (s1 << 7) & 0x9d2c5680UL;
+  s1 ^= (s1 << 15) & 0xefc60000UL;
+  return (s1 ^ (s1 >> 18));
+}
+
+inline MTRand::uint32 MTRand::randInt(const uint32 &n)
+{
+  // Find which bits are used in n
+  // Optimized by Magnus Jonsson (magnus@smartelectronix.com)
+  uint32 used = n;
+  used |= used >> 1;
+  used |= used >> 2;
+  used |= used >> 4;
+  used |= used >> 8;
+  used |= used >> 16;
+
+  // Draw numbers until one is found in [0,n]
+  uint32 i;
+  do
+    i = randInt() & used;  // toss unused bits to shorten search
+  while (i > n);
+  return i;
+}
+
+inline void MTRand::seed(const uint32 oneSeed)
+{
+  // Seed the generator with a simple uint32
+  initialize(oneSeed);
+  reload();
+}
+
+inline void MTRand::seed(uint32 *const bigSeed, const uint32 seedLength)
+{
+  // Seed the generator with an array of uint32's
+  // There are 2^19937-1 possible initial states.  This function allows
+  // all of those to be accessed by providing at least 19937 bits (with a
+  // default seed length of N = 624 uint32's).  Any bits above the lower 32
+  // in each element are discarded.
+  // Just call seed() if you want to get array from /dev/urandom
+  initialize(19650218UL);
+  const unsigned int Nenum = N;
+  int i = 1;
+  uint32 j = 0;
+  int k = (Nenum > seedLength ? Nenum : seedLength);
+  for (; k; --k) {
+    state[i] = state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL);
+    state[i] += (bigSeed[j] & 0xffffffffUL) + j;
+    state[i] &= 0xffffffffUL;
+    ++i;
+    ++j;
+    if (i >= N) {
+      state[0] = state[N - 1];
+      i = 1;
+    }
+    if (j >= seedLength)
+      j = 0;
+  }
+  for (k = N - 1; k; --k) {
+    state[i] = state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL);
+    state[i] -= i;
+    state[i] &= 0xffffffffUL;
+    ++i;
+    if (i >= N) {
+      state[0] = state[N - 1];
+      i = 1;
+    }
+  }
+  state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
+  reload();
+}
+
+inline void MTRand::seed()
+{
+  // Seed the generator with an array from /dev/urandom if available
+  // Otherwise use a hash of time() and clock() values
+
+  // First try getting an array from /dev/urandom
+  FILE *urandom = fopen("/dev/urandom", "rb");
+  if (urandom) {
+    uint32 bigSeed[N];
+    uint32 *s = bigSeed;
+    int i = N;
+    bool success = true;
+    while (success && i--)
+      success = fread(s++, sizeof(uint32), 1, urandom);
+    fclose(urandom);
+    if (success) {
+      seed(bigSeed, N);
+      return;
+    }
+  }
+
+  // Was not successful, so use time() and clock() instead
+  seed(hash(time(NULL), clock()));
+}
+
+inline void MTRand::initialize(const uint32 intseed)
+{
+  // Initialize generator state with seed
+  // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
+  // In previous versions, most significant bits (MSBs) of the seed affect
+  // only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
+  uint32 *s = state;
+  uint32 *r = state;
+  int i = 1;
+  *s++ = intseed & 0xffffffffUL;
+  for (; i < N; ++i) {
+    *s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL;
+    r++;
+  }
+}
+
+inline void MTRand::reload()
+{
+  // Generate N new values in state
+  // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
+  uint32 *p = state;
+  int i;
+  for (i = N - M; i--; ++p)
+    *p = twist(p[M], p[0], p[1]);
+  for (i = M; --i; ++p)
+    *p = twist(p[M - N], p[0], p[1]);
+  *p = twist(p[M - N], p[0], state[0]);
+
+  left = N, pNext = state;
+}
+
+inline MTRand::uint32 MTRand::hash(time_t t, clock_t c)
+{
+  // Get a uint32 from t and c
+  // Better than uint32(x) in case x is floating point in [0,1]
+  // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
+
+  static uint32 differ = 0;  // guarantee time-based seeds will change
+
+  uint32 h1 = 0;
+  unsigned char *p = (unsigned char *)&t;
+  for (size_t i = 0; i < sizeof(t); ++i) {
+    h1 *= std::numeric_limits<unsigned char>::max() + 2U;
+    h1 += p[i];
+  }
+  uint32 h2 = 0;
+  p = (unsigned char *)&c;
+  for (size_t j = 0; j < sizeof(c); ++j) {
+    h2 *= std::numeric_limits<unsigned char>::max() + 2U;
+    h2 += p[j];
+  }
+  return (h1 + differ++) ^ h2;
+}
+
+inline void MTRand::save(uint32 *saveArray) const
+{
+  uint32 *sa = saveArray;
+  const uint32 *s = state;
+  int i = N;
+  for (; i--; *sa++ = *s++) {
+  }
+  *sa = left;
+}
+
+inline void MTRand::load(uint32 *const loadArray)
+{
+  uint32 *s = state;
+  uint32 *la = loadArray;
+  int i = N;
+  for (; i--; *s++ = *la++) {
+  }
+  left = *la;
+  pNext = &state[N - left];
+}
+
+inline std::ostream &operator<<(std::ostream &os, const MTRand &mtrand)
+{
+  const MTRand::uint32 *s = mtrand.state;
+  int i = mtrand.N;
+  for (; i--; os << *s++ << "\t") {
+  }
+  return os << mtrand.left;
+}
+
+inline std::istream &operator>>(std::istream &is, MTRand &mtrand)
+{
+  MTRand::uint32 *s = mtrand.state;
+  int i = mtrand.N;
+  for (; i--; is >> *s++) {
+  }
+  is >> mtrand.left;
+  mtrand.pNext = &mtrand.state[mtrand.N - mtrand.left];
+  return is;
+}
+
+// simple interface to mersenne twister
+class RandomStream {
+ public:
+  inline RandomStream(long seed) : mtr(seed){};
+  ~RandomStream()
+  {
+  }
+
+  /*! get a random number from the stream */
+  inline double getDouble(void)
+  {
+    return mtr.rand();
+  };
+  inline float getFloat(void)
+  {
+    return (float)mtr.rand();
+  };
+
+  inline float getFloat(float min, float max)
+  {
+    return mtr.rand(max - min) + min;
+  };
+  inline float getRandNorm(float mean, float var)
+  {
+    return mtr.randNorm(mean, var);
+  };
+
+#if FLOATINGPOINT_PRECISION == 1
+  inline Real getReal()
+  {
+    return getFloat();
+  }
+
+#else
+  inline Real getReal()
+  {
+    return getDouble();
+  }
+#endif
+
+  inline Vec3 getVec3()
+  {
+    Real a = getReal(), b = getReal(), c = getReal();
+    return Vec3(a, b, c);
+  }
+  inline Vec3 getVec3Norm()
+  {
+    Vec3 a = getVec3();
+    normalize(a);
+    return a;
+  }
+
+ private:
+  MTRand mtr;
+};
+
+}  // namespace Manta
+
+#endif
diff --git a/extern/mantaflow/helper/util/rcmatrix.h b/extern/mantaflow/helper/util/rcmatrix.h
new file mode 100644
index 00000000000..39951cece2e
--- /dev/null
+++ b/extern/mantaflow/helper/util/rcmatrix.h
@@ -0,0 +1,1112 @@
+//
+//  Helper matrix class, RCMatrix.h
+//  Required for PD optimizations (guiding)
+//  Thanks to Ryoichi Ando, and Robert Bridson
+//
+
+#ifndef RCMATRIX3_H
+#define RCMATRIX3_H
+
+#include <iterator>
+#include <cassert>
+#include <vector>
+#include <fstream>
+
+// index type
+#define int_index long long
+
+// link to omp & tbb for now
+#if OPENMP == 1 || TBB == 1
+#  define MANTA_ENABLE_PARALLEL 0
+// allow the preconditioner to be computed in parallel? (can lead to slightly non-deterministic
+// results)
+#  define MANTA_ENABLE_PARALLEL_PC 0
+// use c++11 code?
+#  define MANTA_USE_CPP11 1
+#else
+#  define MANTA_ENABLE_PARALLEL 0
+#  define MANTA_ENABLE_PARALLEL_PC 0
+#  define MANTA_USE_CPP11 0
+#endif
+
+#if MANTA_ENABLE_PARALLEL == 1
+#  include <thread>
+#  include <algorithm>
+
+static const int manta_num_threads = std::thread::hardware_concurrency();
+
+// For clang
+#  define parallel_for(total_size) \
+    { \
+      int_index parallel_array_size = (total_size); \
+      std::vector<std::thread> threads(manta_num_threads); \
+      for (int thread_number = 0; thread_number < manta_num_threads; thread_number++) { \
+  threads[thread_number] = std::thread([&](int_index parallel_array_size, int_index thread_number ) { \
+		for( int_index parallel_index=thread_number; parallel_index < parallel_array_size; parallel_index += manta_num_threads ) {
+
+#  define parallel_end \
+    } \
+    },parallel_array_size,thread_number); \
+    } \
+    for (auto &thread : threads) \
+      thread.join(); \
+    }
+
+#  define parallel_block \
+    { \
+      std::vector<std::thread> threads; \
+      {
+
+#  define do_parallel threads.push_back( std::thread([&]() {
+#  define do_end \
+    } ) );
+
+#  define block_end \
+    } \
+    for (auto &thread : threads) { \
+      thread.join(); \
+    } \
+    }
+
+#else
+
+#  define parallel_for(size) \
+    { \
+      int thread_number = 0; \
+      int_index parallel_index = 0; \
+      for (int_index parallel_index = 0; parallel_index < (int_index)size; parallel_index++) {
+#  define parallel_end \
+    } \
+    thread_number = parallel_index = 0; \
+    }
+
+#  define parallel_block
+#  define do_parallel
+#  define do_end
+#  define block_end
+
+#endif
+
+#include "vectorbase.h"
+
+namespace Manta {
+
+static const unsigned default_expected_none_zeros = 7;
+
+template<class N, class T> struct RCMatrix {
+  struct RowEntry {
+    std::vector<N> index;
+    std::vector<T> value;
+  };
+  RCMatrix() : n(0), expected_none_zeros(default_expected_none_zeros)
+  {
+  }
+  RCMatrix(N size, N expected_none_zeros = default_expected_none_zeros)
+      : n(0), expected_none_zeros(expected_none_zeros)
+  {
+    resize(size);
+  }
+  RCMatrix(const RCMatrix &m) : n(0), expected_none_zeros(default_expected_none_zeros)
+  {
+    init(m);
+  }
+  RCMatrix &operator=(const RCMatrix &m)
+  {
+    expected_none_zeros = m.expected_none_zeros;
+    init(m);
+    return *this;
+  }
+  RCMatrix &operator=(RCMatrix &&m)
+  {
+    matrix = m.matrix;
+    offsets = m.offsets;
+    expected_none_zeros = m.expected_none_zeros;
+    n = m.n;
+    m.n = 0;
+    m.matrix.clear();
+    m.offsets.clear();
+    return *this;
+  }
+  RCMatrix(RCMatrix &&m)
+      : n(m.n), expected_none_zeros(m.expected_none_zeros), matrix(m.matrix), offsets(m.offsets)
+  {
+    m.n = 0;
+    m.matrix.clear();
+    m.offsets.clear();
+  }
+  void init(const RCMatrix &m)
+  {
+    expected_none_zeros = m.expected_none_zeros;
+    resize(m.n);
+    parallel_for(n)
+    {
+      N i = parallel_index;
+      if (m.matrix[i]) {
+        alloc_row(i);
+        matrix[i]->index = m.matrix[i]->index;
+        matrix[i]->value = m.matrix[i]->value;
+      }
+      else {
+        dealloc_row(i);
+      }
+    }
+    parallel_end
+  }
+  ~RCMatrix()
+  {
+    clear();
+  }
+  void clear()
+  {
+    for (N i = 0; i < n; i++) {
+      dealloc_row(i);
+      matrix[i] = NULL;
+      if (offsets.size())
+        offsets[i] = 0;
+    }
+  };
+  bool empty(N i) const
+  {
+    return matrix[i] == NULL;
+  }
+  N row_nonzero_size(N i) const
+  {
+    return matrix[i] == NULL ? 0 : matrix[i]->index.size();
+  }
+  void resize(N size, N expected_none_zeros = 0)
+  {
+    if (!expected_none_zeros) {
+      expected_none_zeros = this->expected_none_zeros;
+    }
+    if (n > size) {
+      // Shrinking
+      for (N i = size ? size - 1 : 0; i < n; i++)
+        dealloc_row(i);
+      matrix.resize(size);
+    }
+    else if (n < size) {
+      // Expanding
+      matrix.resize(size);
+      for (N i = n; i < size; i++) {
+        matrix[i] = NULL;
+        if (offsets.size())
+          offsets[i] = 0;
+      }
+    }
+    n = size;
+  }
+  void alloc_row(N i)
+  {
+    assert(i < n);
+    if (!matrix[i]) {
+      matrix[i] = new RowEntry;
+      matrix[i]->index.reserve(expected_none_zeros);
+      matrix[i]->value.reserve(expected_none_zeros);
+      if (offsets.size())
+        offsets[i] = 0;
+    }
+  }
+  void dealloc_row(N i)
+  {
+    assert(i < n);
+    if (matrix[i]) {
+      if (offsets.empty() || !offsets[i])
+        delete matrix[i];
+      matrix[i] = NULL;
+      if (offsets.size())
+        offsets[i] = 0;
+    }
+  }
+  T operator()(N i, N j) const
+  {
+    assert(i < n);
+    for (Iterator it = row_begin(i); it; ++it) {
+      if (it.index() == j)
+        return it.value();
+    }
+    return T(0.0);
+  }
+  void add_to_element_checked(N i, N j, T val)
+  {
+    if ((i < 0) || (j < 0) || (i >= n) || (j >= n))
+      return;
+    add_to_element(i, j, val);
+  }
+  void add_to_element(N i, N j, T increment_value)
+  {
+    if (std::abs(increment_value) > VECTOR_EPSILON) {
+      assert(i < n);
+      assert(offsets.empty() || offsets[i] == 0);
+      alloc_row(i);
+      std::vector<N> &index = matrix[i]->index;
+      std::vector<T> &value = matrix[i]->value;
+      for (N k = 0; k < (N)index.size(); ++k) {
+        if (index[k] == j) {
+          value[k] += increment_value;
+          return;
+        }
+        else if (index[k] > j) {
+          index.insert(index.begin() + k, j);
+          value.insert(value.begin() + k, increment_value);
+          return;
+        }
+      }
+      index.push_back(j);
+      value.push_back(increment_value);
+    }
+  }
+
+  void set_element(N i, N j, T v)
+  {
+    if (std::abs(v) > VECTOR_EPSILON) {
+      assert(i < n);
+      assert(offsets.empty() || offsets[i] == 0);
+      alloc_row(i);
+      std::vector<N> &index = matrix[i]->index;
+      std::vector<T> &value = matrix[i]->value;
+      for (N k = 0; k < (N)index.size(); ++k) {
+        if (index[k] == j) {
+          value[k] = v;
+          return;
+        }
+        else if (index[k] > j) {
+          index.insert(index.begin() + k, j);
+          value.insert(value.begin() + k, v);
+          return;
+        }
+      }
+      index.push_back(j);
+      value.push_back(v);
+    }
+  }
+
+  // Make sure that j is the biggest column in the row, no duplication allowed
+  void fix_element(N i, N j, T v)
+  {
+    if (std::abs(v) > VECTOR_EPSILON) {
+      assert(i < n);
+      assert(offsets.empty() || offsets[i] == 0);
+      alloc_row(i);
+      std::vector<N> &index = matrix[i]->index;
+      std::vector<T> &value = matrix[i]->value;
+      index.push_back(j);
+      value.push_back(v);
+    }
+  }
+  int_index trim_zero_entries(double e = VECTOR_EPSILON)
+  {
+    std::vector<int_index> deleted_entries(n, 0);
+    parallel_for(n)
+    {
+      N i = parallel_index;
+      if (matrix[i]) {
+        std::vector<N> &index = matrix[i]->index;
+        std::vector<T> &value = matrix[i]->value;
+        N head = 0;
+        N k = 0;
+        for (k = 0; k < index.size(); ++k) {
+          if (std::abs(value[k]) > e) {
+            index[head] = index[k];
+            value[head] = value[k];
+            ++head;
+          }
+        }
+        if (head != k) {
+          index.erase(index.begin() + head, index.end());
+          value.erase(value.begin() + head, value.end());
+          deleted_entries[i] += k - head;
+        }
+        if (!offsets.size() && !head) {
+          remove_row(i);
+        }
+      }
+    }
+    parallel_end
+        //
+        int_index sum_deleted(0);
+    for (int_index i = 0; i < n; i++)
+      sum_deleted += deleted_entries[i];
+    return sum_deleted;
+  }
+  void remove_reference(N i)
+  {
+    if (offsets.size() && offsets[i] && matrix[i]) {
+      RowEntry *save = matrix[i];
+      matrix[i] = new RowEntry;
+      *matrix[i] = *save;
+      for (N &index : matrix[i]->index)
+        index += offsets[i];
+      offsets[i] = 0;
+    }
+  }
+  void remove_row(N i)
+  {
+    dealloc_row(i);
+  }
+  bool is_symmetric(double e = VECTOR_EPSILON) const
+  {
+    std::vector<bool> flags(n, true);
+    parallel_for(n)
+    {
+      N i = parallel_index;
+      bool flag = true;
+      for (Iterator it = row_begin(i); it; ++it) {
+        N index = it.index();
+        T value = it.value();
+        if (std::abs(value) > e) {
+          bool found_entry = false;
+          for (Iterator it_i = row_begin(index); it_i; ++it_i) {
+            if (it_i.index() == i) {
+              found_entry = true;
+              if (std::abs(value - it_i.value()) > e) {
+                flag = false;
+                break;
+              }
+            }
+          }
+          if (!found_entry)
+            flag = false;
+          if (!flag)
+            break;
+        }
+      }
+      flags[i] = flag;
+    }
+    parallel_end for (N i = 0; i < matrix.size(); ++i)
+    {
+      if (!flags[i])
+        return false;
+    }
+    return true;
+  }
+
+  void expand()
+  {
+    if (offsets.empty())
+      return;
+    for (N i = 1; i < n; i++) {
+      if (offsets[i]) {
+        RowEntry *ref = matrix[i];
+        matrix[i] = new RowEntry;
+        *matrix[i] = *ref;
+        for (N j = 0; j < (N)matrix[i]->index.size(); j++) {
+          matrix[i]->index[j] += offsets[i];
+        }
+      }
+    }
+    offsets.resize(0);
+  }
+
+  N column(N i) const
+  {
+    return empty(i) ? 0 : row_begin(i, row_nonzero_size(i) - 1).index();
+  }
+  N getColumnSize() const
+  {
+    N max_column(0);
+    auto column = [&](N i) {
+      N max_column(0);
+      for (Iterator it = row_begin(i); it; ++it)
+        max_column = std::max(max_column, it.index());
+      return max_column + 1;
+    };
+    for (N i = 0; i < n; i++)
+      max_column = std::max(max_column, column(i));
+    return max_column;
+  }
+  N getNonzeroSize() const
+  {
+    N nonzeros(0);
+    for (N i = 0; i < n; ++i) {
+      nonzeros += row_nonzero_size(i);
+    }
+    return nonzeros;
+  }
+  class Iterator : std::iterator<std::input_iterator_tag, T> {
+   public:
+    Iterator(const RowEntry *rowEntry, N k, N offset) : rowEntry(rowEntry), k(k), offset(offset)
+    {
+    }
+    operator bool() const
+    {
+      return rowEntry != NULL && k < (N)rowEntry->index.size();
+    }
+    Iterator &operator++()
+    {
+      ++k;
+      return *this;
+    }
+    T value() const
+    {
+      return rowEntry->value[k];
+    }
+    N index() const
+    {
+      return rowEntry->index[k] + offset;
+    }
+    N index_raw() const
+    {
+      return rowEntry->index[k];
+    }
+    N size() const
+    {
+      return rowEntry == NULL ? 0 : rowEntry->index.size();
+    }
+
+   protected:
+    const RowEntry *rowEntry;
+    N k, offset;
+  };
+  Iterator row_begin(N n, N k = 0) const
+  {
+    return Iterator(matrix[n], k, offsets.size() ? offsets[n] : 0);
+  }
+  class DynamicIterator : public Iterator {
+   public:
+    DynamicIterator(RowEntry *rowEntry, N k, N offset)
+        : rowEntry(rowEntry), Iterator(rowEntry, k, offset)
+    {
+    }
+    void setValue(T value)
+    {
+      rowEntry->value[Iterator::k] = value;
+    }
+    void setIndex(N index)
+    {
+      rowEntry->index[Iterator::k] = index;
+    }
+
+   protected:
+    RowEntry *rowEntry;
+  };
+  DynamicIterator dynamic_row_begin(N n, N k = 0)
+  {
+    N offset = offsets.size() ? offsets[n] : 0;
+    if (offset) {
+      printf("---- Warning ----\n");
+      printf("Dynamic iterator is not allowed for referenced rows.\n");
+      printf("You should be very careful otherwise this causes some bugs.\n");
+      printf(
+          "We encourage you that you convert this row into a raw format, then loop over it...\n");
+      printf("-----------------\n");
+      exit(0);
+    }
+    return DynamicIterator(matrix[n], k, offset);
+  }
+  RCMatrix transpose(N rowsize = 0,
+                     unsigned expected_none_zeros = default_expected_none_zeros) const
+  {
+    if (!rowsize)
+      rowsize = getColumnSize();
+    RCMatrix result(rowsize, expected_none_zeros);
+    for (N i = 0; i < n; i++) {
+      for (Iterator it = row_begin(i); it; ++it)
+        result.fix_element(it.index(), i, it.value());
+    }
+#if MANTA_USE_CPP11 == 1
+    return std::move(result);
+#else
+    return result;
+#endif
+  }
+
+  RCMatrix getKtK() const
+  {
+    RCMatrix m = transpose();
+    RCMatrix result(n, expected_none_zeros);
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      for (Iterator it_A = m.row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        assert(j < n);
+        T a = it_A.value();
+        if (std::abs(a) > VECTOR_EPSILON) {
+          for (Iterator it_B = row_begin(j); it_B; ++it_B) {
+            // result.add_to_element(i,it_B.index(),it_B.value()*a);
+            double value = it_B.value() * a;
+            if (std::abs(value) > VECTOR_EPSILON)
+              result.add_to_element(i, it_B.index(), value);
+          }
+        }
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix operator*(const RCMatrix &m) const
+  {
+    RCMatrix result(n, expected_none_zeros);
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      for (Iterator it_A = row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        assert(j < m.n);
+        T a = it_A.value();
+        if (std::abs(a) > VECTOR_EPSILON) {
+          for (Iterator it_B = m.row_begin(j); it_B; ++it_B) {
+            // result.add_to_element(i,it_B.index(),it_B.value()*a);
+            double value = it_B.value() * a;
+            if (std::abs(value) > VECTOR_EPSILON)
+              result.add_to_element(i, it_B.index(), value);
+          }
+        }
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix sqrt() const
+  {
+    RCMatrix result(n, expected_none_zeros);
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      for (Iterator it_A = row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        result.set_element(i, j, std::sqrt(it_A.value()));
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix operator*(const double k) const
+  {
+    RCMatrix result(n, expected_none_zeros);
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      for (Iterator it_A = row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        result.add_to_element(i, j, it_A.value() * k);
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix applyKernel(const RCMatrix &kernel, const int nx, const int ny) const
+  {
+    RCMatrix result(n, expected_none_zeros);
+    // find center position of kernel (half of kernel size)
+    int kCols = kernel.n, kRows = kernel.n, rows = nx, cols = ny;
+    int kCenterX = kCols / 2;
+    int kCenterY = kRows / 2;
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      if (i >= rows)
+        break;
+      for (Iterator it_A = row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        if (j >= cols)
+          break;
+        for (int m = 0; m < kRows; ++m) {    // kernel rows
+          int mm = kRows - 1 - m;            // row index of flipped kernel
+          for (int n = 0; n < kCols; ++n) {  // kernel columns
+            int nn = kCols - 1 - n;          // column index of flipped kernel
+            // index of input signal, used for checking boundary
+            int ii = i + (m - kCenterY);
+            int jj = j + (n - kCenterX);
+            // ignore input samples which are out of bound
+            if (ii >= 0 && ii < rows && jj >= 0 && jj < cols)
+              result.add_to_element(i, j, (*this)(ii, jj) * kernel(mm, nn));
+          }
+        }
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix applyHorizontalKernel(const RCMatrix &kernel, const int nx, const int ny) const
+  {
+    RCMatrix result(n, expected_none_zeros);
+    // find center position of kernel (half of kernel size)
+    int kCols = kernel.n, kRows = 1, rows = nx, cols = ny;
+    int kCenterX = kCols / 2;
+    int kCenterY = kRows / 2;
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      if (i >= rows)
+        break;
+      for (Iterator it_A = row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        if (j >= cols)
+          break;
+        for (int m = 0; m < kRows; ++m) {    // kernel rows
+          int mm = kRows - 1 - m;            // row index of flipped kernel
+          for (int n = 0; n < kCols; ++n) {  // kernel columns
+            int nn = kCols - 1 - n;          // column index of flipped kernel
+            // index of input signal, used for checking boundary
+            int ii = i + (m - kCenterY);
+            int jj = j + (n - kCenterX);
+            // ignore input samples which are out of bound
+            if (ii >= 0 && ii < rows && jj >= 0 && jj < cols)
+              result.add_to_element(i, j, (*this)(ii, jj) * kernel(mm, nn));
+          }
+        }
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix applyVerticalKernel(const RCMatrix &kernel, const int nx, const int ny) const
+  {
+    RCMatrix result(n, expected_none_zeros);
+    // find center position of kernel (half of kernel size)
+    int kCols = 1, kRows = kernel.n, rows = nx, cols = ny;
+    int kCenterX = kCols / 2;
+    int kCenterY = kRows / 2;
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      if (i >= rows)
+        break;
+      for (Iterator it_A = row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        if (j >= cols)
+          break;
+        for (int m = 0; m < kRows; ++m) {    // kernel rows
+          int mm = kRows - 1 - m;            // row index of flipped kernel
+          for (int n = 0; n < kCols; ++n) {  // kernel columns
+            int nn = kCols - 1 - n;          // column index of flipped kernel
+            // index of input signal, used for checking boundary
+            int ii = i + (m - kCenterY);
+            int jj = j + (n - kCenterX);
+            // ignore input samples which are out of bound
+            if (ii >= 0 && ii < rows && jj >= 0 && jj < cols)
+              result.add_to_element(i, j, (*this)(ii, jj) * kernel(mm, nn));
+          }
+        }
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix applySeparableKernel(const RCMatrix &kernelH,
+                                const RCMatrix &kernelV,
+                                const int nx,
+                                const int ny) const
+  {
+    return applyHorizontalKernel(kernelH, nx, ny).applyVerticalKernel(kernelV, nx, ny);
+  }
+
+  RCMatrix applySeparableKernelTwice(const RCMatrix &kernelH,
+                                     const RCMatrix &kernelV,
+                                     const int nx,
+                                     const int ny) const
+  {
+    return applySeparableKernel(kernelH, kernelV, nx, ny)
+        .applySeparableKernel(kernelH, kernelV, nx, ny);
+  }
+
+  std::vector<T> operator*(const std::vector<T> &rhs) const
+  {
+    std::vector<T> result(n, 0.0);
+    multiply(rhs, result);
+#if MANTA_USE_CPP11 == 1
+    return std::move(result);
+#else
+    return result;
+#endif
+  }
+  void multiply(const std::vector<T> &rhs, std::vector<T> &result) const
+  {
+    result.resize(n);
+    for (N i = 0; i < n; i++) {
+      T new_value = 0.0;
+      for (Iterator it = row_begin(i); it; ++it) {
+        N j_index = it.index();
+        assert(j_index < rhs.size());
+        new_value += rhs[j_index] * it.value();
+      }
+      result[i] = new_value;
+    }
+  }
+  RCMatrix operator+(const RCMatrix &m) const
+  {
+    RCMatrix A(*this);
+    return std::move(A.add(m));
+  }
+  RCMatrix &add(const RCMatrix &m)
+  {
+    if (m.n > n)
+      resize(m.n);
+    parallel_for(m.n)
+    {
+      N i = parallel_index;
+      for (Iterator it = m.row_begin(i); it; ++it) {
+        add_to_element(i, it.index(), it.value());
+      }
+    }
+    parallel_end return *this;
+  }
+  RCMatrix operator-(const RCMatrix &m) const
+  {
+    RCMatrix A(*this);
+    return std::move(A.sub(m));
+  }
+  RCMatrix &sub(const RCMatrix &m)
+  {
+    if (m.n > n)
+      resize(m.n);
+    parallel_for(m.n)
+    {
+      N i = parallel_index;
+      for (Iterator it = m.row_begin(i); it; ++it) {
+        add_to_element(i, it.index(), -it.value());
+      }
+    }
+    parallel_end return *this;
+  }
+  RCMatrix &replace(const RCMatrix &m, int rowInd, int colInd)
+  {
+    if (m.n > n)
+      resize(m.n);
+    parallel_for(m.n)
+    {
+      N i = parallel_index;
+      for (Iterator it = m.row_begin(i); it; ++it) {
+        set_element(i + rowInd, it.index() + colInd, it.value());
+      }
+    }
+    parallel_end return *this;
+  }
+  Real max_residual(const std::vector<T> &lhs, const std::vector<T> &rhs) const
+  {
+    std::vector<T> r = operator*(lhs);
+    Real max_residual = 0.0;
+    for (N i = 0; i < rhs.size(); i++) {
+      if (!empty(i))
+        max_residual = std::max(max_residual, std::abs(r[i] - rhs[i]));
+    }
+    return max_residual;
+  }
+  std::vector<T> residual_vector(const std::vector<T> &lhs, const std::vector<T> &rhs) const
+  {
+    std::vector<T> result = operator*(lhs);
+    assert(result.size() == rhs.size());
+    for (N i = 0; i < result.size(); i++) {
+      result[i] = std::abs(result[i] - rhs[i]);
+    }
+#if MANTA_USE_CPP11 == 1
+    return std::move(result);
+#else
+    return result;
+#endif
+  }
+  T norm() const
+  {
+    T result(0.0);
+    for (N i = 0; i < n; ++i) {
+      for (Iterator it = row_begin(i); it; ++it) {
+        result = std::max(result, std::abs(it.value()));
+      }
+    }
+    return result;
+  }
+
+  T norm_L2_sqr() const
+  {
+    T result(0.0);
+    for (N i = 0; i < n; ++i) {
+      for (Iterator it = row_begin(i); it; ++it) {
+        result += (it.value()) * (it.value());
+      }
+    }
+    return result;
+  }
+
+  void write_matlab(std::ostream &output,
+                    unsigned int rows,
+                    unsigned int columns,
+                    const char *variable_name)
+  {
+    output << variable_name << "=sparse([";
+    for (N i = 0; i < n; ++i) {
+      if (matrix[i]) {
+        const std::vector<N> &index = matrix[i]->index;
+        for (N j = 0; j < (N)index.size(); ++j) {
+          output << i + 1 << " ";
+        }
+      }
+    }
+    output << "],...\n  [";
+    for (N i = 0; i < n; ++i) {
+      if (matrix[i]) {
+        const std::vector<N> &index = matrix[i]->index;
+        for (N j = 0; j < (N)index.size(); ++j) {
+          output << index[j] + (offsets.empty() ? 0 : offsets[i]) + 1 << " ";
+        }
+      }
+    }
+    output << "],...\n  [";
+    for (N i = 0; i < n; ++i) {
+      if (matrix[i]) {
+        const std::vector<T> &value = matrix[i]->value;
+        for (N j = 0; j < value.size(); ++j) {
+          output << value[j] << " ";
+        }
+      }
+    }
+    output << "], " << rows << ", " << columns << ");" << std::endl;
+  };
+  void export_matlab(std::string filename, std::string name)
+  {
+    // Export this matrix
+    std::ofstream file;
+    file.open(filename.c_str());
+    write_matlab(file, n, getColumnSize(), name.c_str());
+    file.close();
+  }
+  void print_readable(std::string name, bool printNonZero = true)
+  {
+    std::cout << name << " \n";
+    for (int i = 0; i < n; ++i) {
+      for (int j = 0; j < n; ++j) {
+        if (printNonZero) {
+          if ((*this)(i, j) == 0) {
+            std::cout << "  .";
+            continue;
+          }
+        }
+        else {
+          if ((*this)(i, j) == 0) {
+            continue;
+          }
+        }
+
+        if ((*this)(i, j) >= 0)
+          std::cout << " ";
+        std::cout << " " << (*this)(i, j);
+      }
+      std::cout << " \n";
+    }
+  }
+  ///
+  N n;
+  N expected_none_zeros;
+  std::vector<RowEntry *> matrix;
+  std::vector<int> offsets;
+};
+
+template<class N, class T>
+static inline RCMatrix<N, T> operator*(const std::vector<T> &diagonal, const RCMatrix<N, T> &A)
+{
+  RCMatrix<N, T> result(A);
+  parallel_for(result.n)
+  {
+    N row(parallel_index);
+    for (auto it = result.dynamic_row_begin(row); it; ++it) {
+      it.setValue(it.value() * diagonal[row]);
+    }
+  }
+  parallel_end return std::move(result);
+}
+
+template<class N, class T> struct RCFixedMatrix {
+  std::vector<N> rowstart;
+  std::vector<N> index;
+  std::vector<T> value;
+  N n;
+  N max_rowlength;
+  //
+  RCFixedMatrix() : n(0), max_rowlength(0)
+  {
+  }
+  RCFixedMatrix(const RCMatrix<N, T> &matrix)
+  {
+    n = matrix.n;
+    rowstart.resize(n + 1);
+    rowstart[0] = 0;
+    max_rowlength = 0;
+    for (N i = 0; i < n; i++) {
+      if (!matrix.empty(i)) {
+        rowstart[i + 1] = rowstart[i] + matrix.row_nonzero_size(i);
+        max_rowlength = std::max(max_rowlength, rowstart[i + 1] - rowstart[i]);
+      }
+      else {
+        rowstart[i + 1] = rowstart[i];
+      }
+    }
+    value.resize(rowstart[n]);
+    index.resize(rowstart[n]);
+    N j = 0;
+    for (N i = 0; i < n; i++) {
+      for (typename RCMatrix<N, T>::Iterator it = matrix.row_begin(i); it; ++it) {
+        value[j] = it.value();
+        index[j] = it.index();
+        ++j;
+      }
+    }
+  }
+  class Iterator : std::iterator<std::input_iterator_tag, T> {
+   public:
+    Iterator(N start, N end, const std::vector<N> &index, const std::vector<T> &value)
+        : index_array(index), value_array(value), k(start), start(start), end(end)
+    {
+    }
+    operator bool() const
+    {
+      return k < end;
+    }
+    Iterator &operator++()
+    {
+      ++k;
+      return *this;
+    }
+    T value() const
+    {
+      return value_array[k];
+    }
+    N index() const
+    {
+      return index_array[k];
+    }
+    N size() const
+    {
+      return end - start;
+    }
+
+   private:
+    const std::vector<N> &index_array;
+    const std::vector<T> &value_array;
+    N k, start, end;
+  };
+  Iterator row_begin(N n) const
+  {
+    return Iterator(rowstart[n], rowstart[n + 1], index, value);
+  }
+  std::vector<T> operator*(const std::vector<T> &rhs) const
+  {
+    std::vector<T> result(n, 0.0);
+    multiply(rhs, result);
+#if MANTA_USE_CPP11 == 1
+    return std::move(result);
+#else
+    return result;
+#endif
+  }
+  void multiply(const std::vector<T> &rhs, std::vector<T> &result) const
+  {
+    result.resize(n);
+    parallel_for(n)
+    {
+      N i = parallel_index;
+      T new_value = 0.0;
+      for (Iterator it = row_begin(i); it; ++it) {
+        N j_index = it.index();
+        assert(j_index < rhs.size());
+        new_value += rhs[j_index] * it.value();
+      }
+      result[i] = new_value;
+    }
+    parallel_end
+  }
+  RCMatrix<N, T> operator*(const RCFixedMatrix &m) const
+  {
+    RCMatrix<N, T> result(n, max_rowlength);
+    // Run in parallel
+    parallel_for(result.n)
+    {
+      N i = parallel_index;
+      for (Iterator it_A = row_begin(i); it_A; ++it_A) {
+        N j = it_A.index();
+        assert(j < m.n);
+        T a = it_A.value();
+        if (std::abs(a) > VECTOR_EPSILON) {
+          for (Iterator it_B = m.row_begin(j); it_B; ++it_B) {
+            result.add_to_element(i, it_B.index(), it_B.value() * a);
+          }
+        }
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+
+  RCMatrix<N, T> toRCMatrix() const
+  {
+    RCMatrix<N, T> result(n, 0);
+    parallel_for(n)
+    {
+      N i = parallel_index;
+      N size = rowstart[i + 1] - rowstart[i];
+      result.matrix[i] = new typename RCMatrix<N, T>::RowEntry;
+      result.matrix[i]->index.resize(size);
+      result.matrix[i]->value.resize(size);
+      for (N j = 0; j < size; j++) {
+        result.matrix[i]->index[j] = index[rowstart[i] + j];
+        result.matrix[i]->value[j] = value[rowstart[i] + j];
+      }
+    }
+    parallel_end
+#if MANTA_USE_CPP11 == 1
+        return std::move(result);
+#else
+        return result;
+#endif
+  }
+};
+
+typedef RCMatrix<int, Real> Matrix;
+typedef RCFixedMatrix<int, Real> FixedMatrix;
+
+}  // namespace Manta
+
+#undef parallel_for
+#undef parallel_end
+
+#undef parallel_block
+#undef do_parallel
+#undef do_end
+#undef block_end
+
+#endif
diff --git a/extern/mantaflow/helper/util/simpleimage.cpp b/extern/mantaflow/helper/util/simpleimage.cpp
new file mode 100644
index 00000000000..9846fa5bd96
--- /dev/null
+++ b/extern/mantaflow/helper/util/simpleimage.cpp
@@ -0,0 +1,312 @@
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2014 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
+ *
+ * Simple image IO
+ *
+ ******************************************************************************/
+
+#include "vectorbase.h"
+#include "simpleimage.h"
+
+namespace Manta {
+
+// write rectangle to ppm
+bool SimpleImage::writePpm(
+    std::string filename, int minx, int miny, int maxx, int maxy, bool invertXY)
+{
+  int w = maxx - minx;
+  int h = maxy - miny;
+
+  if (w <= 0 || h <= 0 || w > mSize[0] || h > mSize[1]) {
+    errMsg("SimpleImage::WritePPM Invalid rect: w="
+           << w << ", h=" << h << ", size=" << mSize[0] << "," << mSize[1] << " min/max: " << minx
+           << "," << miny << " to " << maxx << "," << maxy << ", resetting... ");
+    minx = miny = 0;
+    maxx = mSize[0] - 1;
+    maxy = mSize[1] - 1;
+    w = mSize[0] - 1;
+    h = mSize[1] - 1;
+  }
+
+  FILE *fp = fopen(filename.c_str(), "wb");
+  if (fp == NULL) {
+    errMsg("SimpleImage::WritePPM Unable to open '" << filename << "' for writing");
+    return false;
+  }
+  fprintf(fp, "P6\n%d %d\n255\n", w, h);
+
+  int pixCnt = 0;
+  for (int j = maxy - 1; j >= miny; j--)
+    for (int i = minx; i < maxx; i++) {
+      unsigned char col[3];
+      for (int l = 0; l < 3; l++) {
+        float val;
+        if (invertXY)
+          val = (float)get(j, i)[l];
+        else
+          val = (float)get(i, j)[l];
+
+        val = clamp(val, (float)0., (float)1.);
+        col[l] = (unsigned char)(255. * val);
+      }
+      // col[1] = col[2] = col[0];
+      // if (fwrite(col,1,3, fp) != 3) errMsg("SimpleImage::writePpm fwrite failed");
+      fwrite(col, 1, 3, fp);
+      pixCnt++;
+      // fprintf(stderr,"%d %d %d \n",col[0],i,j);
+    }
+
+  fclose(fp);
+  // debMsg("WritePPM Wrote '"<<filename<<"', region="<<minx<<","<<miny<<" to
+  // "<<maxx<<","<<maxy<<"; "<<pixCnt, 1);
+
+  return true;
+}
+
+bool SimpleImage::writePpm(std::string filename)
+{
+  return writePpm(filename, 0, 0, getSize()[0], getSize()[1]);
+}
+
+// read in a ppm file, and init the image accordingly
+bool SimpleImage::initFromPpm(std::string filename)
+{
+  // maximum length of a line of text
+  const int MAXLINE = 1024;
+
+  int filetype = 0;
+  enum { PGM, PPM };  // possible file types
+
+  FILE *fp;
+  char line[MAXLINE];
+  int size, rowsize;
+
+  // Read in file type
+  fp = fopen(filename.c_str(), "rb");
+  if (!fp) {
+    if (mAbortOnError)
+      debMsg("SimpleImage Error - unable to open file '" << filename << "' for reading", 1);
+    return 0;
+  }
+
+  // 1st line: PPM or PGM
+  if (fgets(line, MAXLINE, fp) == NULL) {
+    if (mAbortOnError)
+      debMsg("SimpleImage::initFromPpm fgets failed", 1);
+    return 0;
+  }
+
+  if (line[1] == '5')
+    filetype = PGM;
+  else if (line[1] == '6')
+    filetype = PPM;
+  else {
+    if (mAbortOnError)
+      debMsg("SimpleImage Error: need PPM or PGM file as input!", 1);
+    return 0;
+  }
+
+  // Read in width and height, & allocate space
+  // 2nd line: width height
+  if (fgets(line, MAXLINE, fp) == NULL) {
+    if (mAbortOnError)
+      errMsg("SimpleImage::initFromPpm fgets failed");
+    return 0;
+  }
+  int windW = 0, windH = 0;  // size of the window on the screen
+  int intsFound = sscanf(line, "%d %d", &windW, &windH);
+  if (intsFound == 1) {
+    // only X found, search on next line as well for Y...
+    if (sscanf(line, "%d", &windH) != 1) {
+      if (mAbortOnError)
+        errMsg("initFromPpm Ppm dimensions not found!" << windW << "," << windH);
+      return 0;
+    }
+    else {
+      // ok, found 2 lines
+      // debMsg("initFromPpm Ppm dimensions found!"<<windW<<","<<windH, 1);
+    }
+  }
+  else if (intsFound == 2) {
+    // ok!
+  }
+  else {
+    if (mAbortOnError)
+      errMsg("initFromPpm Ppm dimensions not found at all!" << windW << "," << windH);
+    return 0;
+  }
+
+  if (filetype == PGM) {
+    size = windH * windW;  // greymap: 1 byte per pixel
+    rowsize = windW;
+  }
+  else {
+    // filetype == PPM
+    size = windH * windW * 3;  // pixmap: 3 bytes per pixel
+    rowsize = windW * 3;
+  }
+
+  unsigned char *pic = new unsigned char[size];  // (GLubyte *)malloc (size);
+
+  // Read in maximum value (ignore) , could be scanned with sscanf as well, but this should be
+  // 255... 3rd line
+  if (fgets(line, MAXLINE, fp) == NULL) {
+    if (mAbortOnError)
+      errMsg("SimpleImage::initFromPpm fgets failed");
+    return 0;
+  }
+
+  // Read in the pixel array row-by-row: 1st row = top scanline */
+  unsigned char *ptr = NULL;
+  ptr = &pic[(windH - 1) * rowsize];
+  for (int i = windH; i > 0; i--) {
+    assertMsg(fread((void *)ptr, 1, rowsize, fp) == rowsize,
+              "SimpleImage::initFromPpm couldn't read data");
+    ptr -= rowsize;
+  }
+
+  // init image
+  this->init(windW, windH);
+  if (filetype == PGM) {
+    // grayscale
+    for (int i = 0; i < windW; i++) {
+      for (int j = 0; j < windH; j++) {
+        double r = (double)pic[(j * windW + i) * 1 + 0] / 255.;
+        (*this)(i, j) = Vec3(r, r, r);
+      }
+    }
+  }
+  else {
+    // convert grid to RGB vec's
+    for (int i = 0; i < windW; i++) {
+      for (int j = 0; j < windH; j++) {
+        // return mpData[y*mSize[0]+x];
+        double r = (double)pic[(j * windW + i) * 3 + 0] / 255.;
+        double g = (double)pic[(j * windW + i) * 3 + 1] / 255.;
+        double b = (double)pic[(j * windW + i) * 3 + 2] / 255.;
+
+        //(*this)(i,j) = Vec3(r,g,b);
+
+        // RGB values have to be rotated to get the right colors!?
+        // this might also be an artifact of photoshop export...?
+        (*this)(i, j) = Vec3(g, b, r);
+      }
+    }
+  }
+
+  delete[] pic;
+  fclose(fp);
+  return 1;
+}
+
+// check index is valid
+bool SimpleImage::indexIsValid(int i, int j)
+{
+  if (i < 0)
+    return false;
+  if (j < 0)
+    return false;
+  if (i >= mSize[0])
+    return false;
+  if (j >= mSize[1])
+    return false;
+  return true;
+}
+
+};  // namespace Manta
+
+//*****************************************************************************
+
+#include "grid.h"
+namespace Manta {
+
+// simple shaded output , note requires grid functionality!
+static void gridPrecompLight(const Grid<Real> &density, Grid<Real> &L, Vec3 light = Vec3(1, 1, 1))
+{
+  FOR_IJK(density)
+  {
+    Vec3 n = getGradient(density, i, j, k) * -1.;
+    normalize(n);
+
+    Real d = dot(light, n);
+    L(i, j, k) = d;
+  }
+}
+
+// simple shading with pre-computed gradient
+static inline void shadeCell(
+    Vec3 &dst, int shadeMode, Real src, Real light, int depthPos, Real depthInv)
+{
+  switch (shadeMode) {
+
+    case 1: {
+      // surfaces
+      Vec3 ambient = Vec3(0.1, 0.1, 0.1);
+      Vec3 diffuse = Vec3(0.9, 0.9, 0.9);
+      Real alpha = src;
+
+      // different color for depth?
+      diffuse[0] *= ((Real)depthPos * depthInv) * 0.7 + 0.3;
+      diffuse[1] *= ((Real)depthPos * depthInv) * 0.7 + 0.3;
+
+      Vec3 col = ambient + diffuse * light;
+
+      // img( 0+i, j ) = (1.-alpha) * img( 0+i, j ) + alpha * col;
+      dst = (1. - alpha) * dst + alpha * col;
+    } break;
+
+    default: {
+      // volumetrics / smoke
+      dst += depthInv * Vec3(src, src, src);
+    } break;
+  }
+}
+
+//! helper to project a grid intro an image (used for ppm export and GUI displauy)
+void projectImg(SimpleImage &img, const Grid<Real> &val, int shadeMode = 0, Real scale = 1.)
+{
+  Vec3i s = val.getSize();
+  Vec3 si = Vec3(1. / (Real)s[0], 1. / (Real)s[1], 1. / (Real)s[2]);
+
+  // init image size
+  int imgSx = s[0];
+  if (val.is3D())
+    imgSx += s[2] + s[0];  // mult views in 3D
+  img.init(imgSx, std::max(s[0], std::max(s[1], s[2])));
+
+  // precompute lighting
+  Grid<Real> L(val);
+  gridPrecompLight(val, L, Vec3(1, 1, 1));
+
+  FOR_IJK(val)
+  {
+    Vec3i idx(i, j, k);
+    shadeCell(img(0 + i, j), shadeMode, val(idx), L(idx), k, si[2]);
+  }
+
+  if (val.is3D()) {
+
+    FOR_IJK(val)
+    {
+      Vec3i idx(i, j, k);
+      shadeCell(img(s[0] + k, j), shadeMode, val(idx), L(idx), i, si[0]);
+    }
+
+    FOR_IJK(val)
+    {
+      Vec3i idx(i, j, k);
+      shadeCell(img(s[0] + s[2] + i, k), shadeMode, val(idx), L(idx), j, si[1]);
+    }
+
+  }  // 3d
+
+  img.mapRange(1. / scale);
+}
+
+};  // namespace Manta
diff --git a/extern/mantaflow/helper/util/simpleimage.h b/extern/mantaflow/helper/util/simpleimage.h
new file mode 100644
index 00000000000..d7e88b83f74
--- /dev/null
+++ b/extern/mantaflow/helper/util/simpleimage.h
@@ -0,0 +1,205 @@
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2014 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
+ *
+ * Simple image IO
+ *
+ ******************************************************************************/
+
+#ifndef MANTA_SIMPLEIMAGE_H
+#define MANTA_SIMPLEIMAGE_H
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include "manta.h"
+#include "vectorbase.h"
+
+namespace Manta {
+
+//*****************************************************************************
+// simple 2d image class
+// template<class Scalar>
+class SimpleImage {
+ public:
+  // cons/des
+  SimpleImage() : mSize(-1), mpData(NULL), mAbortOnError(true){};
+  virtual ~SimpleImage()
+  {
+    if (mpData)
+      delete[] mpData;
+  };
+
+  //! set to constant
+  void reset(Real val = 0.)
+  {
+    const Vec3 v = Vec3(val);
+    for (int i = 0; i < mSize[0] * mSize[1]; i++)
+      mpData[i] = v;
+  }
+  //! init memory & reset to zero
+  void init(int x, int y)
+  {
+    mSize = Vec3i(x, y, 0);
+    mpData = new Vec3[x * y];
+    reset();
+  };
+
+  inline bool checkIndex(int x, int y)
+  {
+    if ((x < 0) || (y < 0) || (x > mSize[0] - 1) || (y > mSize[1] - 1)) {
+      errMsg("SimpleImage::operator() Invalid access to " << x << "," << y << ", size=" << mSize);
+      return false;
+    }
+    return true;
+  }
+
+  // access element
+  inline Vec3 &operator()(int x, int y)
+  {
+    DEBUG_ONLY(checkIndex(x, y));
+    return mpData[y * mSize[0] + x];
+  };
+  inline Vec3 &get(int x, int y)
+  {
+    return (*this)(x, y);
+  }
+  inline Vec3 &getMap(int x, int y, int z, int axis)
+  {
+    int i = x, j = y;
+    if (axis == 1)
+      j = z;
+    if (axis == 0) {
+      i = y;
+      j = z;
+    }
+    return get(i, j);
+  }
+
+  // output as string, debug
+  std::string toString()
+  {
+    std::ostringstream out;
+
+    for (int j = 0; j < mSize[1]; j++) {
+      for (int i = 0; i < mSize[0]; i++) {
+        // normal zyx order */
+        out << (*this)(i, j);
+        out << " ";
+      }
+      // if (format)
+      out << std::endl;
+    }
+
+    return out.str();
+  }
+
+  // multiply all values by f
+  void add(Vec3 f)
+  {
+    for (int j = 0; j < mSize[1]; j++)
+      for (int i = 0; i < mSize[0]; i++) {
+        get(i, j) += f;
+      }
+  }
+  // multiply all values by f
+  void multiply(Real f)
+  {
+    for (int j = 0; j < mSize[1]; j++)
+      for (int i = 0; i < mSize[0]; i++) {
+        get(i, j) *= f;
+      }
+  }
+  // map 0-f to 0-1 range, clamp
+  void mapRange(Real f)
+  {
+    for (int j = 0; j < mSize[1]; j++)
+      for (int i = 0; i < mSize[0]; i++) {
+        get(i, j) /= f;
+        for (int c = 0; c < 3; ++c)
+          get(i, j)[c] = clamp(get(i, j)[c], (Real)0., (Real)1.);
+      }
+  }
+
+  // normalize max values
+  void normalizeMax()
+  {
+    Real max = normSquare(get(0, 0));
+    for (int j = 0; j < mSize[1]; j++)
+      for (int i = 0; i < mSize[0]; i++) {
+        if (normSquare(get(i, j)) > max)
+          max = normSquare(get(i, j));
+      }
+    max = sqrt(max);
+    Real invMax = 1. / max;
+    for (int j = 0; j < mSize[1]; j++)
+      for (int i = 0; i < mSize[0]; i++) {
+        get(i, j) *= invMax;
+      }
+  };
+
+  // normalize min and max values
+  void normalizeMinMax()
+  {
+    Real max = normSquare(get(0, 0));
+    Real min = max;
+    for (int j = 0; j < mSize[1]; j++)
+      for (int i = 0; i < mSize[0]; i++) {
+        if (normSquare(get(i, j)) > max)
+          max = normSquare(get(i, j));
+        if (normSquare(get(i, j)) < min)
+          min = normSquare(get(i, j));
+      }
+    max = sqrt(max);
+    min = sqrt(min);
+    Real factor = 1. / (max - min);
+    for (int j = 0; j < mSize[1]; j++)
+      for (int i = 0; i < mSize[0]; i++) {
+        get(i, j) -= min;
+        get(i, j) *= factor;
+      }
+  };
+
+  void setAbortOnError(bool set)
+  {
+    mAbortOnError = set;
+  }
+
+  // ppm in/output
+
+  // write whole image
+  bool writePpm(std::string filename);
+  // write rectangle to ppm
+  bool writePpm(
+      std::string filename, int minx, int miny, int maxx, int maxy, bool invertXY = false);
+  // read in a ppm file, and init the image accordingly
+  bool initFromPpm(std::string filename);
+
+  // check index is valid
+  bool indexIsValid(int i, int j);
+
+  //! access
+  inline Vec3i getSize() const
+  {
+    return mSize;
+  }
+
+ protected:
+  //! size
+  Vec3i mSize;
+  //! data
+  Vec3 *mpData;
+  // make errors fatal, or continue?
+  bool mAbortOnError;
+
+};  // SimpleImage
+
+};  // namespace Manta
+
+#endif
diff --git a/extern/mantaflow/helper/util/solvana.h b/extern/mantaflow/helper/util/solvana.h
new file mode 100644
index 00000000000..9dc1ec83654
--- /dev/null
+++ b/extern/mantaflow/helper/util/solvana.h
@@ -0,0 +1,214 @@
+/******************************************************************************
+ *
+ * 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
+ *
+ * Analytical solutions to some problems
+ * generated using MATLAB symbolic math ccode
+ *
+ ******************************************************************************/
+
+#ifndef _SOLVANA_H
+#define _SOLVANA_H
+
+//! solves the equation [e1 e2 e3; 1 1 1]*x = g using least squares
+inline void SolveOverconstraint34(float e1x,
+                                  float e1y,
+                                  float e1z,
+                                  float e2x,
+                                  float e2y,
+                                  float e2z,
+                                  float e3x,
+                                  float e3y,
+                                  float e3z,
+                                  float g1,
+                                  float g2,
+                                  float g3,
+                                  float &x1,
+                                  float &x2,
+                                  float &x3)
+{
+  float e1x2 = e1x * e1x, e1y2 = e1y * e1y, e1z2 = e1z * e1z;
+  float e2x2 = e2x * e2x, e2y2 = e2y * e2y, e2z2 = e2z * e2z;
+  float e3x2 = e3x * e3x, e3y2 = e3y * e3y, e3z2 = e3z * e3z;
+  float e1xy = e1x * e1y, e1xz = e1x * e1z, e1yz = e1y * e1z;
+  float e2xy = e2x * e2y, e2xz = e2x * e2z, e2yz = e2y * e2z;
+  float e3xy = e3x * e3y, e3xz = e3x * e3z, e3yz = e3y * e3z;
+  float e12x = e1x * e2x, e12y = e1y * e2y, e12z = e1z * e2z;
+  float e13x = e1x * e3x, e13y = e1y * e3y, e13z = e1z * e3z;
+  float e23x = e2x * e3x, e23y = e2y * e3y, e23z = e2z * e3z;
+  float t1543 = e3y2 * e2x2;
+  float t1544 = e3x2 * e2y2;
+  float t1545 = e3z2 * e2x2;
+  float t1546 = e3x2 * e2z2;
+  float t1547 = e3z2 * e2y2;
+  float t1548 = e3y2 * e2z2;
+  float t1549 = e2y2 * e1x2;
+  float t1550 = e2x2 * e1y2;
+  float t1551 = e2z2 * e1x2;
+  float t1552 = e2x2 * e1z2;
+  float t1553 = e2z2 * e1y2;
+  float t1554 = e2y2 * e1z2;
+  float t1555 = e3y2 * e1x2;
+  float t1556 = e3x2 * e1y2;
+  float t1557 = e3z2 * e1x2;
+  float t1558 = e3x2 * e1z2;
+  float t1559 = e3z2 * e1y2;
+  float t1560 = e3y2 * e1z2;
+  float t1561 = e3z2 * e2y2 * e1x2;
+  float t1562 = e3y2 * e2z2 * e1x2;
+  float t1563 = e3z2 * e2x2 * e1y2;
+  float t1564 = e3x2 * e2z2 * e1y2;
+  float t1565 = e3y2 * e2x2 * e1z2;
+  float t1566 = e3x2 * e2y2 * e1z2;
+  float t1567 = e1xy * e2x * e3y * 2.0;
+  float t1568 = e1xy * e2y * e3x * 2.0;
+  float t1569 = e1xz * e2x * e3z * 2.0;
+  float t1570 = e1xz * e2z * e3x * 2.0;
+  float t1571 = e1yz * e2y * e3z * 2.0;
+  float t1572 = e1yz * e2z * e3y * 2.0;
+  float t1573 = e1x * e2xy * e3y * 2.0;
+  float t1574 = e1y * e2xy * e3x * 2.0;
+  float t1575 = e1x * e2xz * e3z * 2.0;
+  float t1576 = e1z * e2xz * e3x * 2.0;
+  float t1577 = e1y * e2yz * e3z * 2.0;
+  float t1578 = e1z * e2yz * e3y * 2.0;
+  float t1579 = e1x * e2y * e3xy * 2.0;
+  float t1580 = e1y * e2x * e3xy * 2.0;
+  float t1581 = e1x * e2z * e3xz * 2.0;
+  float t1582 = e1z * e2x * e3xz * 2.0;
+  float t1583 = e1y * e2z * e3yz * 2.0;
+  float t1584 = e1z * e2y * e3yz * 2.0;
+  float t1585 = e1xy * e2xz * e3yz * 2.0;
+  float t1586 = e1xy * e2yz * e3xz * 2.0;
+  float t1587 = e1xz * e2xy * e3yz * 2.0;
+  float t1588 = e1xz * e2yz * e3xy * 2.0;
+  float t1589 = e1yz * e2xy * e3xz * 2.0;
+  float t1590 = e1yz * e2xz * e3xy * 2.0;
+  float t1596 = e12x * e3y2 * 2.0;
+  float t1597 = e13x * e2y2 * 2.0;
+  float t1598 = e23x * e1y2 * 2.0;
+  float t1599 = e12x * e3z2 * 2.0;
+  float t1600 = e13x * e2z2 * 2.0;
+  float t1601 = e12y * e3x2 * 2.0;
+  float t1602 = e13y * e2x2 * 2.0;
+  float t1603 = e23y * e1x2 * 2.0;
+  float t1604 = e23x * e1z2 * 2.0;
+  float t1605 = e12y * e3z2 * 2.0;
+  float t1606 = e13y * e2z2 * 2.0;
+  float t1607 = e12z * e3x2 * 2.0;
+  float t1608 = e13z * e2x2 * 2.0;
+  float t1609 = e23z * e1x2 * 2.0;
+  float t1610 = e23y * e1z2 * 2.0;
+  float t1611 = e12z * e3y2 * 2.0;
+  float t1612 = e13z * e2y2 * 2.0;
+  float t1613 = e23z * e1y2 * 2.0;
+  float t1614 = e1xy * e2xy * 2.0;
+  float t1615 = e1xz * e2xz * 2.0;
+  float t1616 = e1yz * e2yz * 2.0;
+  float t1617 = e1xy * e3xy * 2.0;
+  float t1618 = e1xz * e3xz * 2.0;
+  float t1619 = e1yz * e3yz * 2.0;
+  float t1620 = e2xy * e3xy * 2.0;
+  float t1621 = e2xz * e3xz * 2.0;
+  float t1622 = e2yz * e3yz * 2.0;
+  float t1623 = e1xy * e2xy * e3z2 * 2.0;
+  float t1624 = e1xz * e2xz * e3y2 * 2.0;
+  float t1625 = e1yz * e2yz * e3x2 * 2.0;
+  float t1626 = e1xy * e3xy * e2z2 * 2.0;
+  float t1627 = e1xz * e3xz * e2y2 * 2.0;
+  float t1628 = e1yz * e3yz * e2x2 * 2.0;
+  float t1629 = e2xy * e3xy * e1z2 * 2.0;
+  float t1630 = e2xz * e3xz * e1y2 * 2.0;
+  float t1631 = e2yz * e3yz * e1x2 * 2.0;
+  float t1591 = t1550 + t1551 + t1560 + t1543 + t1552 + t1561 + t1570 + t1544 + t1553 + t1562 +
+                t1571 + t1580 + t1545 + t1554 + t1563 + t1572 + t1581 + t1590 + t1546 + t1555 +
+                t1564 + t1573 + t1582 + t1547 + t1556 + t1565 + t1574 + t1583 + t1548 + t1557 +
+                t1566 + t1575 + t1584 + t1549 + t1558 + t1567 + t1576 + t1585 + t1559 + t1568 +
+                t1577 + t1586 + t1569 + t1578 + t1587 - t1596 + t1579 + t1588 - t1597 + t1589 -
+                t1598 - t1599 - t1600 - t1601 - t1610 - t1602 - t1611 - t1620 - t1603 - t1612 -
+                t1621 - t1630 - t1604 - t1613 - t1622 - t1631 - t1605 - t1614 - t1623 - t1606 -
+                t1615 - t1624 - t1607 - t1616 - t1625 - t1608 - t1617 - t1626 - t1609 - t1618 -
+                t1627 - t1619 - t1628 - t1629;
+  float t1592 = 1.0 / t1591;
+  float t1635 = e13x * e2y2;
+  float t1636 = e13x * e2z2;
+  float t1637 = e13y * e2x2;
+  float t1638 = e13y * e2z2;
+  float t1639 = e13z * e2x2;
+  float t1640 = e13z * e2y2;
+  float t1653 = e23x * 2.0;
+  float t1654 = e23y * 2.0;
+  float t1655 = e23z * 2.0;
+  float t1641 = e3x2 + e3z2 + e3y2 + e2y2 + t1543 + e2z2 + t1544 + e2x2 + t1545 + t1546 + t1547 +
+                t1548 - t1620 - t1621 - t1622 - t1653 - t1654 - t1655;
+  float t1642 = e12x * e3y2;
+  float t1643 = e12x * e3z2;
+  float t1644 = e12y * e3x2;
+  float t1645 = e12y * e3z2;
+  float t1646 = e12z * e3x2;
+  float t1647 = e12z * e3y2;
+  float t1656 = e1x * e2y * e3xy;
+  float t1657 = e1y * e2x * e3xy;
+  float t1658 = e1x * e2z * e3xz;
+  float t1659 = e1z * e2x * e3xz;
+  float t1660 = e1y * e2z * e3yz;
+  float t1661 = e1z * e2y * e3yz;
+  float t1648 = e3x2 + e3z2 + e3y2 - e13x - e13y - e13z + e12x - e23y + e12y + t1642 - e23z -
+                t1660 + e12z + t1643 - t1661 + t1644 + t1645 + t1646 + t1647 - t1656 - t1657 -
+                e23x - t1658 - t1659;
+  float t1679 = e1x * e2xy * e3y;
+  float t1680 = e1y * e2xy * e3x;
+  float t1681 = e1x * e2xz * e3z;
+  float t1682 = e1z * e2xz * e3x;
+  float t1683 = e1y * e2yz * e3z;
+  float t1684 = e1z * e2yz * e3y;
+  float t1652 = e2y2 + e2z2 + e2x2 + e13x + e13y + e13z + t1640 - e12x - e23y - e12y - e23z -
+                e12z + t1635 - t1680 + t1636 - t1681 + t1637 - t1682 + t1638 - t1683 + t1639 -
+                t1684 - e23x - t1679;
+  float t1662 = e23x * e1y2;
+  float t1663 = e23y * e1x2;
+  float t1664 = e23x * e1z2;
+  float t1665 = e23z * e1x2;
+  float t1666 = e23y * e1z2;
+  float t1667 = e23z * e1y2;
+  float t1670 = e1xy * e2x * e3y;
+  float t1671 = e1xy * e2y * e3x;
+  float t1672 = e1xz * e2x * e3z;
+  float t1673 = e1xz * e2z * e3x;
+  float t1674 = e1yz * e2y * e3z;
+  float t1675 = e1yz * e2z * e3y;
+  float t1668 = e1x2 + e1y2 + e1z2 - e13x - e13y - e13z - e12x + e23y - e12y + e23z - e12z -
+                t1670 + t1662 - t1671 + t1663 - t1672 + t1664 - t1673 + t1665 - t1674 + t1666 -
+                t1675 + e23x + t1667;
+  float t1676 = e13x * 2.0;
+  float t1677 = e13y * 2.0;
+  float t1678 = e13z * 2.0;
+  float t1669 = e3x2 + e3z2 + e3y2 + t1560 + e1x2 + t1555 + e1y2 + t1556 + e1z2 + t1557 + t1558 +
+                t1559 - t1617 - t1618 - t1619 - t1676 - t1677 - t1678;
+  float t1686 = e12x * 2.0;
+  float t1687 = e12y * 2.0;
+  float t1688 = e12z * 2.0;
+  float t1685 = t1550 + t1551 + e2y2 + t1552 + e2z2 + t1553 + e2x2 + t1554 + e1x2 + e1y2 + e1z2 +
+                t1549 - t1614 - t1615 - t1616 - t1686 - t1687 - t1688;
+  x1 = -g2 * (-e1y * t1592 * t1641 + e2y * t1592 * t1648 + e3y * t1592 * t1652) -
+       g3 * (-e1z * t1592 * t1641 + e2z * t1592 * t1648 + e3z * t1592 * t1652) -
+       g1 * (-e1x * t1592 * t1641 + e2x * t1592 * t1648 +
+             e3x * t1592 *
+                 (e2y2 + e2z2 + e2x2 + e13x + e13y + e13z + t1640 + t1635 + t1636 + t1637 + t1638 +
+                  t1639 - e12x - e12y - e12z - e23x - e23y - e23z - e1x * e2xy * e3y -
+                  e1y * e2xy * e3x - e1x * e2xz * e3z - e1z * e2xz * e3x - e1y * e2yz * e3z -
+                  e1z * e2yz * e3y));
+  x2 = -g1 * (e1x * t1592 * t1648 - e2x * t1592 * t1669 + e3x * t1592 * t1668) -
+       g2 * (e1y * t1592 * t1648 - e2y * t1592 * t1669 + e3y * t1592 * t1668) -
+       g3 * (e1z * t1592 * t1648 - e2z * t1592 * t1669 + e3z * t1592 * t1668);
+  x3 = -g1 * (e1x * t1592 * t1652 + e2x * t1592 * t1668 - e3x * t1592 * t1685) -
+       g2 * (e1y * t1592 * t1652 + e2y * t1592 * t1668 - e3y * t1592 * t1685) -
+       g3 * (e1z * t1592 * t1652 + e2z * t1592 * t1668 - e3z * t1592 * t1685);
+}
+
+#endif
\ No newline at end of file
diff --git a/extern/mantaflow/helper/util/vector4d.cpp b/extern/mantaflow/helper/util/vector4d.cpp
new file mode 100644
index 00000000000..d342df607f5
--- /dev/null
+++ b/extern/mantaflow/helper/util/vector4d.cpp
@@ -0,0 +1,50 @@
+/******************************************************************************
+ *
+ * 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
+ *
+ * Basic vector class
+ *
+ ******************************************************************************/
+
+#include "vector4d.h"
+
+using namespace std;
+
+namespace Manta {
+
+template<> const Vector4D<int> Vector4D<int>::Zero(0, 0, 0, 0);
+template<> const Vector4D<float> Vector4D<float>::Zero(0.f, 0.f, 0.f, 0.f);
+template<> const Vector4D<double> Vector4D<double>::Zero(0., 0., 0., 0.);
+template<>
+const Vector4D<float> Vector4D<float>::Invalid(numeric_limits<float>::quiet_NaN(),
+                                               numeric_limits<float>::quiet_NaN(),
+                                               numeric_limits<float>::quiet_NaN(),
+                                               numeric_limits<float>::quiet_NaN());
+template<>
+const Vector4D<double> Vector4D<double>::Invalid(numeric_limits<double>::quiet_NaN(),
+                                                 numeric_limits<double>::quiet_NaN(),
+                                                 numeric_limits<double>::quiet_NaN(),
+                                                 numeric_limits<double>::quiet_NaN());
+template<> bool Vector4D<float>::isValid() const
+{
+  return !c_isnan(x) && !c_isnan(y) && !c_isnan(z) && !c_isnan(t);
+}
+template<> bool Vector4D<double>::isValid() const
+{
+  return !c_isnan(x) && !c_isnan(y) && !c_isnan(z) && !c_isnan(t);
+}
+
+//! Specialization for readable ints
+template<> std::string Vector4D<int>::toString() const
+{
+  char buf[256];
+  snprintf(buf, 256, "[%d,%d,%d,%d]", (*this)[0], (*this)[1], (*this)[2], (*this)[3]);
+  return std::string(buf);
+}
+
+}  // namespace Manta
diff --git a/extern/mantaflow/helper/util/vector4d.h b/extern/mantaflow/helper/util/vector4d.h
new file mode 100644
index 00000000000..c3d72ac8aff
--- /dev/null
+++ b/extern/mantaflow/helper/util/vector4d.h
@@ -0,0 +1,515 @@
+/******************************************************************************
+ *
+ * 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
+ *
+ * 4D vector class
+ *
+ ******************************************************************************/
+
+#ifndef _VECTOR4D_H
+#define _VECTOR4D_H
+
+#include "vectorbase.h"
+
+namespace Manta {
+
+//! Basic inlined vector class
+template<class S> class Vector4D {
+ public:
+  //! Constructor
+  inline Vector4D() : x(0), y(0), z(0), t(0)
+  {
+  }
+
+  //! Copy-Constructor
+  inline Vector4D(const Vector4D<S> &v) : x(v.x), y(v.y), z(v.z), t(v.t)
+  {
+  }
+
+  //! Copy-Constructor
+  inline Vector4D(const float *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3])
+  {
+  }
+
+  //! Copy-Constructor
+  inline Vector4D(const double *v) : x((S)v[0]), y((S)v[1]), z((S)v[2]), t((S)v[3])
+  {
+  }
+
+  //! Construct a vector from one S
+  inline Vector4D(S v) : x(v), y(v), z(v), t(v)
+  {
+  }
+
+  //! Construct a vector from three Ss
+  inline Vector4D(S vx, S vy, S vz, S vw) : x(vx), y(vy), z(vz), t(vw)
+  {
+  }
+
+  // Operators
+
+  //! Assignment operator
+  inline const Vector4D<S> &operator=(const Vector4D<S> &v)
+  {
+    x = v.x;
+    y = v.y;
+    z = v.z;
+    t = v.t;
+    return *this;
+  }
+  //! Assignment operator
+  inline const Vector4D<S> &operator=(S s)
+  {
+    x = y = z = t = s;
+    return *this;
+  }
+  //! Assign and add operator
+  inline const Vector4D<S> &operator+=(const Vector4D<S> &v)
+  {
+    x += v.x;
+    y += v.y;
+    z += v.z;
+    t += v.t;
+    return *this;
+  }
+  //! Assign and add operator
+  inline const Vector4D<S> &operator+=(S s)
+  {
+    x += s;
+    y += s;
+    z += s;
+    t += s;
+    return *this;
+  }
+  //! Assign and sub operator
+  inline const Vector4D<S> &operator-=(const Vector4D<S> &v)
+  {
+    x -= v.x;
+    y -= v.y;
+    z -= v.z;
+    t -= v.t;
+    return *this;
+  }
+  //! Assign and sub operator
+  inline const Vector4D<S> &operator-=(S s)
+  {
+    x -= s;
+    y -= s;
+    z -= s;
+    t -= s;
+    return *this;
+  }
+  //! Assign and mult operator
+  inline const Vector4D<S> &operator*=(const Vector4D<S> &v)
+  {
+    x *= v.x;
+    y *= v.y;
+    z *= v.z;
+    t *= v.t;
+    return *this;
+  }
+  //! Assign and mult operator
+  inline const Vector4D<S> &operator*=(S s)
+  {
+    x *= s;
+    y *= s;
+    z *= s;
+    t *= s;
+    return *this;
+  }
+  //! Assign and div operator
+  inline const Vector4D<S> &operator/=(const Vector4D<S> &v)
+  {
+    x /= v.x;
+    y /= v.y;
+    z /= v.z;
+    t /= v.t;
+    return *this;
+  }
+  //! Assign and div operator
+  inline const Vector4D<S> &operator/=(S s)
+  {
+    x /= s;
+    y /= s;
+    z /= s;
+    t /= s;
+    return *this;
+  }
+  //! Negation operator
+  inline Vector4D<S> operator-() const
+  {
+    return Vector4D<S>(-x, -y, -z, -t);
+  }
+
+  //! Get smallest component
+  // inline S min() const { return ( x<y ) ? ( ( x<z ) ? x:z ) : ( ( y<z ) ? y:z ); }
+  //! Get biggest component
+  // inline S max() const { return ( x>y ) ? ( ( x>z ) ? x:z ) : ( ( y>z ) ? y:z ); }
+
+  //! Test if all components are zero
+  inline bool empty()
+  {
+    return x == 0 && y == 0 && z == 0 && t == 0;
+  }
+
+  //! access operator
+  inline S &operator[](unsigned int i)
+  {
+    return value[i];
+  }
+  //! constant access operator
+  inline const S &operator[](unsigned int i) const
+  {
+    return value[i];
+  }
+
+  //! debug output vector to a string
+  std::string toString() const;
+
+  //! test if nans are present
+  bool isValid() const;
+
+  //! actual values
+  union {
+    S value[4];
+    struct {
+      S x;
+      S y;
+      S z;
+      S t;
+    };
+    struct {
+      S X;
+      S Y;
+      S Z;
+      S T;
+    };
+  };
+
+  // zero element
+  static const Vector4D<S> Zero, Invalid;
+
+ protected:
+};
+
+//************************************************************************
+// Additional operators
+//************************************************************************
+
+//! Addition operator
+template<class S> inline Vector4D<S> operator+(const Vector4D<S> &v1, const Vector4D<S> &v2)
+{
+  return Vector4D<S>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z, v1.t + v2.t);
+}
+//! Addition operator
+template<class S, class S2> inline Vector4D<S> operator+(const Vector4D<S> &v, S2 s)
+{
+  return Vector4D<S>(v.x + s, v.y + s, v.z + s, v.t + s);
+}
+//! Addition operator
+template<class S, class S2> inline Vector4D<S> operator+(S2 s, const Vector4D<S> &v)
+{
+  return Vector4D<S>(v.x + s, v.y + s, v.z + s, v.t + s);
+}
+
+//! Subtraction operator
+template<class S> inline Vector4D<S> operator-(const Vector4D<S> &v1, const Vector4D<S> &v2)
+{
+  return Vector4D<S>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.t - v2.t);
+}
+//! Subtraction operator
+template<class S, class S2> inline Vector4D<S> operator-(const Vector4D<S> &v, S2 s)
+{
+  return Vector4D<S>(v.x - s, v.y - s, v.z - s, v.t - s);
+}
+//! Subtraction operator
+template<class S, class S2> inline Vector4D<S> operator-(S2 s, const Vector4D<S> &v)
+{
+  return Vector4D<S>(s - v.x, s - v.y, s - v.z, s - v.t);
+}
+
+//! Multiplication operator
+template<class S> inline Vector4D<S> operator*(const Vector4D<S> &v1, const Vector4D<S> &v2)
+{
+  return Vector4D<S>(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z, v1.t * v2.t);
+}
+//! Multiplication operator
+template<class S, class S2> inline Vector4D<S> operator*(const Vector4D<S> &v, S2 s)
+{
+  return Vector4D<S>(v.x * s, v.y * s, v.z * s, v.t * s);
+}
+//! Multiplication operator
+template<class S, class S2> inline Vector4D<S> operator*(S2 s, const Vector4D<S> &v)
+{
+  return Vector4D<S>(s * v.x, s * v.y, s * v.z, s * v.t);
+}
+
+//! Division operator
+template<class S> inline Vector4D<S> operator/(const Vector4D<S> &v1, const Vector4D<S> &v2)
+{
+  return Vector4D<S>(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z, v1.t / v2.t);
+}
+//! Division operator
+template<class S, class S2> inline Vector4D<S> operator/(const Vector4D<S> &v, S2 s)
+{
+  return Vector4D<S>(v.x / s, v.y / s, v.z / s, v.t / s);
+}
+//! Division operator
+template<class S, class S2> inline Vector4D<S> operator/(S2 s, const Vector4D<S> &v)
+{
+  return Vector4D<S>(s / v.x, s / v.y, s / v.z, s / v.t);
+}
+
+//! Comparison operator
+template<class S> inline bool operator==(const Vector4D<S> &s1, const Vector4D<S> &s2)
+{
+  return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z && s1.t == s2.t;
+}
+
+//! Comparison operator
+template<class S> inline bool operator!=(const Vector4D<S> &s1, const Vector4D<S> &s2)
+{
+  return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z || s1.t != s2.t;
+}
+
+//************************************************************************
+// External functions
+//************************************************************************
+
+//! Dot product
+template<class S> inline S dot(const Vector4D<S> &t, const Vector4D<S> &v)
+{
+  return t.x * v.x + t.y * v.y + t.z * v.z + t.t * v.t;
+}
+
+//! Cross product
+/*template<class S>
+inline Vector4D<S> cross ( const Vector4D<S> &t, const Vector4D<S> &v ) {
+  NYI Vector4D<S> cp (
+    ( ( t.y*v.z ) - ( t.z*v.y ) ),
+    ( ( t.z*v.x ) - ( t.x*v.z ) ),
+    ( ( t.x*v.y ) - ( t.y*v.x ) ) );
+  return cp;
+}*/
+
+//! Compute the magnitude (length) of the vector
+template<class S> inline S norm(const Vector4D<S> &v)
+{
+  S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t;
+  return (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ? 1. : sqrt(l);
+}
+
+//! Compute squared magnitude
+template<class S> inline S normSquare(const Vector4D<S> &v)
+{
+  return v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t;
+}
+
+//! Returns a normalized vector
+template<class S> inline Vector4D<S> getNormalized(const Vector4D<S> &v)
+{
+  S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t;
+  if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON)
+    return v; /* normalized "enough"... */
+  else if (l > VECTOR_EPSILON * VECTOR_EPSILON) {
+    S fac = 1. / sqrt(l);
+    return Vector4D<S>(v.x * fac, v.y * fac, v.z * fac, v.t * fac);
+  }
+  else
+    return Vector4D<S>((S)0);
+}
+
+//! Compute the norm of the vector and normalize it.
+/*! \return The value of the norm */
+template<class S> inline S normalize(Vector4D<S> &v)
+{
+  S norm;
+  S l = v.x * v.x + v.y * v.y + v.z * v.z + v.t * v.t;
+  if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) {
+    norm = 1.;
+  }
+  else if (l > VECTOR_EPSILON * VECTOR_EPSILON) {
+    norm = sqrt(l);
+    v *= 1. / norm;
+  }
+  else {
+    v = Vector4D<S>::Zero;
+    norm = 0.;
+  }
+  return (S)norm;
+}
+
+//! Outputs the object in human readable form as string
+template<class S> std::string Vector4D<S>::toString() const
+{
+  char buf[256];
+  snprintf(buf,
+           256,
+           "[%+4.6f,%+4.6f,%+4.6f,%+4.6f]",
+           (double)(*this)[0],
+           (double)(*this)[1],
+           (double)(*this)[2],
+           (double)(*this)[3]);
+  // for debugging, optionally increase precision:
+  // snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) (
+  // *this ) [1], ( double ) ( *this ) [2], ( double ) ( *this ) [3] );
+  return std::string(buf);
+}
+
+template<> std::string Vector4D<int>::toString() const;
+
+//! Outputs the object in human readable form to stream
+template<class S> std::ostream &operator<<(std::ostream &os, const Vector4D<S> &i)
+{
+  os << i.toString();
+  return os;
+}
+
+//! Reads the contents of the object from a stream
+template<class S> std::istream &operator>>(std::istream &is, Vector4D<S> &i)
+{
+  char c;
+  char dummy[4];
+  is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> dummy >> i[3] >> c;
+  return is;
+}
+
+/**************************************************************************/
+// Define default vector alias
+/**************************************************************************/
+
+//! 3D vector class of type Real (typically float)
+typedef Vector4D<Real> Vec4;
+
+//! 3D vector class of type int
+typedef Vector4D<int> Vec4i;
+
+//! convert to Real Vector
+template<class T> inline Vec4 toVec4(T v)
+{
+  return Vec4(v[0], v[1], v[2], v[3]);
+}
+template<class T> inline Vec4i toVec4i(T v)
+{
+  return Vec4i(v[0], v[1], v[2], v[3]);
+}
+
+/**************************************************************************/
+// Specializations for common math functions
+/**************************************************************************/
+
+template<> inline Vec4 clamp<Vec4>(const Vec4 &a, const Vec4 &b, const Vec4 &c)
+{
+  return Vec4(
+      clamp(a.x, b.x, c.x), clamp(a.y, b.y, c.y), clamp(a.z, b.z, c.z), clamp(a.t, b.t, c.t));
+}
+template<> inline Vec4 safeDivide<Vec4>(const Vec4 &a, const Vec4 &b)
+{
+  return Vec4(
+      safeDivide(a.x, b.x), safeDivide(a.y, b.y), safeDivide(a.z, b.z), safeDivide(a.t, b.t));
+}
+template<> inline Vec4 nmod<Vec4>(const Vec4 &a, const Vec4 &b)
+{
+  return Vec4(nmod(a.x, b.x), nmod(a.y, b.y), nmod(a.z, b.z), nmod(a.t, b.t));
+}
+
+/**************************************************************************/
+// 4d interpolation (note only 4d here, 2d/3d interpolations are in interpol.h)
+/**************************************************************************/
+
+#define BUILD_INDEX_4D \
+  Real px = pos.x - 0.5f, py = pos.y - 0.5f, pz = pos.z - 0.5f, pt = pos.t - 0.5f; \
+  int xi = (int)px; \
+  int yi = (int)py; \
+  int zi = (int)pz; \
+  int ti = (int)pt; \
+  Real s1 = px - (Real)xi, s0 = 1. - s1; \
+  Real t1 = py - (Real)yi, t0 = 1. - t1; \
+  Real f1 = pz - (Real)zi, f0 = 1. - f1; \
+  Real g1 = pt - (Real)ti, g0 = 1. - g1; \
+  /* clamp to border */ \
+  if (px < 0.) { \
+    xi = 0; \
+    s0 = 1.0; \
+    s1 = 0.0; \
+  } \
+  if (py < 0.) { \
+    yi = 0; \
+    t0 = 1.0; \
+    t1 = 0.0; \
+  } \
+  if (pz < 0.) { \
+    zi = 0; \
+    f0 = 1.0; \
+    f1 = 0.0; \
+  } \
+  if (pt < 0.) { \
+    ti = 0; \
+    g0 = 1.0; \
+    g1 = 0.0; \
+  } \
+  if (xi >= size.x - 1) { \
+    xi = size.x - 2; \
+    s0 = 0.0; \
+    s1 = 1.0; \
+  } \
+  if (yi >= size.y - 1) { \
+    yi = size.y - 2; \
+    t0 = 0.0; \
+    t1 = 1.0; \
+  } \
+  if (zi >= size.z - 1) { \
+    zi = size.z - 2; \
+    f0 = 0.0; \
+    f1 = 1.0; \
+  } \
+  if (ti >= size.t - 1) { \
+    ti = size.t - 2; \
+    g0 = 0.0; \
+    g1 = 1.0; \
+  } \
+  const int sX = 1; \
+  const int sY = size.x;
+
+static inline void checkIndexInterpol4d(const Vec4i &size, int idx)
+{
+  if (idx < 0 || idx > size.x * size.y * size.z * size.t) {
+    std::ostringstream s;
+    s << "Grid interpol4d dim " << size << " : index " << idx << " out of bound ";
+    errMsg(s.str());
+  }
+}
+
+template<class T>
+inline T interpol4d(
+    const T *data, const Vec4i &size, const IndexInt sZ, const IndexInt sT, const Vec4 &pos)
+{
+  BUILD_INDEX_4D
+  IndexInt idx = (IndexInt)xi + sY * (IndexInt)yi + sZ * (IndexInt)zi + sT * (IndexInt)ti;
+  DEBUG_ONLY(checkIndexInterpol4d(size, idx));
+  DEBUG_ONLY(checkIndexInterpol4d(size, idx + sX + sY + sZ + sT));
+
+  return (((data[idx] * t0 + data[idx + sY] * t1) * s0 +
+           (data[idx + sX] * t0 + data[idx + sX + sY] * t1) * s1) *
+              f0 +
+          ((data[idx + sZ] * t0 + data[idx + sY + sZ] * t1) * s0 +
+           (data[idx + sX + sZ] * t0 + data[idx + sX + sY + sZ] * t1) * s1) *
+              f1) *
+             g0 +
+         (((data[idx + sT] * t0 + data[idx + sT + sY] * t1) * s0 +
+           (data[idx + sT + sX] * t0 + data[idx + sT + sX + sY] * t1) * s1) *
+              f0 +
+          ((data[idx + sT + sZ] * t0 + data[idx + sT + sY + sZ] * t1) * s0 +
+           (data[idx + sT + sX + sZ] * t0 + data[idx + sT + sX + sY + sZ] * t1) * s1) *
+              f1) *
+             g1;
+}
+
+};  // namespace Manta
+
+#endif
diff --git a/extern/mantaflow/helper/util/vectorbase.cpp b/extern/mantaflow/helper/util/vectorbase.cpp
new file mode 100644
index 00000000000..413ae086d1c
--- /dev/null
+++ b/extern/mantaflow/helper/util/vectorbase.cpp
@@ -0,0 +1,49 @@
+/******************************************************************************
+ *
+ * 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
+ *
+ * Basic vector class
+ *
+ ******************************************************************************/
+
+#include "vectorbase.h"
+
+using namespace std;
+
+namespace Manta {
+
+template<> const Vector3D<int> Vector3D<int>::Zero(0, 0, 0);
+template<> const Vector3D<float> Vector3D<float>::Zero(0.f, 0.f, 0.f);
+template<> const Vector3D<double> Vector3D<double>::Zero(0., 0., 0.);
+template<>
+const Vector3D<float> Vector3D<float>::Invalid(numeric_limits<float>::quiet_NaN(),
+                                               numeric_limits<float>::quiet_NaN(),
+                                               numeric_limits<float>::quiet_NaN());
+template<>
+const Vector3D<double> Vector3D<double>::Invalid(numeric_limits<double>::quiet_NaN(),
+                                                 numeric_limits<double>::quiet_NaN(),
+                                                 numeric_limits<double>::quiet_NaN());
+
+template<> bool Vector3D<float>::isValid() const
+{
+  return !c_isnan(x) && !c_isnan(y) && !c_isnan(z);
+}
+template<> bool Vector3D<double>::isValid() const
+{
+  return !c_isnan(x) && !c_isnan(y) && !c_isnan(z);
+}
+
+//! Specialization for readable ints
+template<> std::string Vector3D<int>::toString() const
+{
+  char buf[256];
+  snprintf(buf, 256, "[%d,%d,%d]", (*this)[0], (*this)[1], (*this)[2]);
+  return std::string(buf);
+}
+
+}  // namespace Manta
diff --git a/extern/mantaflow/helper/util/vectorbase.h b/extern/mantaflow/helper/util/vectorbase.h
new file mode 100644
index 00000000000..a3135431eb3
--- /dev/null
+++ b/extern/mantaflow/helper/util/vectorbase.h
@@ -0,0 +1,679 @@
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2011-2016 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
+ *
+ * Basic vector class
+ *
+ ******************************************************************************/
+
+#ifndef _VECTORBASE_H
+#define _VECTORBASE_H
+
+// get rid of windos min/max defines
+#if defined(WIN32) || defined(_WIN32)
+#  define NOMINMAX
+#endif
+
+#include <stdio.h>
+#include <string>
+#include <limits>
+#include <iostream>
+#include "general.h"
+
+// if min/max are still around...
+#if defined(WIN32) || defined(_WIN32)
+#  undef min
+#  undef max
+#endif
+
+// redefine usage of some windows functions
+#if defined(WIN32) || defined(_WIN32)
+#  ifndef snprintf
+#    define snprintf _snprintf
+#  endif
+#endif
+
+// use which fp-precision? 1=float, 2=double
+#ifndef FLOATINGPOINT_PRECISION
+#  define FLOATINGPOINT_PRECISION 1
+#endif
+
+// VECTOR_EPSILON is the minimal vector length
+// In order to be able to discriminate floating point values near zero, and
+// to be sure not to fail a comparison because of roundoff errors, use this
+// value as a threshold.
+#if FLOATINGPOINT_PRECISION == 1
+typedef float Real;
+#  define VECTOR_EPSILON (1e-6f)
+#else
+typedef double Real;
+#  define VECTOR_EPSILON (1e-10)
+#endif
+
+#ifndef M_PI
+#  define M_PI 3.1415926536
+#endif
+#ifndef M_E
+#  define M_E 2.7182818284
+#endif
+
+namespace Manta {
+
+//! Basic inlined vector class
+template<class S> class Vector3D {
+ public:
+  //! Constructor
+  inline Vector3D() : x(0), y(0), z(0)
+  {
+  }
+
+  //! Copy-Constructor
+  inline Vector3D(const Vector3D<S> &v) : x(v.x), y(v.y), z(v.z)
+  {
+  }
+
+  //! Copy-Constructor
+  inline Vector3D(const int *v) : x((S)v[0]), y((S)v[1]), z((S)v[2])
+  {
+  }
+
+  //! Copy-Constructor
+  inline Vector3D(const float *v) : x((S)v[0]), y((S)v[1]), z((S)v[2])
+  {
+  }
+
+  //! Copy-Constructor
+  inline Vector3D(const double *v) : x((S)v[0]), y((S)v[1]), z((S)v[2])
+  {
+  }
+
+  //! Construct a vector from one S
+  inline Vector3D(S v) : x(v), y(v), z(v)
+  {
+  }
+
+  //! Construct a vector from three Ss
+  inline Vector3D(S vx, S vy, S vz) : x(vx), y(vy), z(vz)
+  {
+  }
+
+  // Operators
+
+  //! Assignment operator
+  inline const Vector3D<S> &operator=(const Vector3D<S> &v)
+  {
+    x = v.x;
+    y = v.y;
+    z = v.z;
+    return *this;
+  }
+  //! Assignment operator
+  inline const Vector3D<S> &operator=(S s)
+  {
+    x = y = z = s;
+    return *this;
+  }
+  //! Assign and add operator
+  inline const Vector3D<S> &operator+=(const Vector3D<S> &v)
+  {
+    x += v.x;
+    y += v.y;
+    z += v.z;
+    return *this;
+  }
+  //! Assign and add operator
+  inline const Vector3D<S> &operator+=(S s)
+  {
+    x += s;
+    y += s;
+    z += s;
+    return *this;
+  }
+  //! Assign and sub operator
+  inline const Vector3D<S> &operator-=(const Vector3D<S> &v)
+  {
+    x -= v.x;
+    y -= v.y;
+    z -= v.z;
+    return *this;
+  }
+  //! Assign and sub operator
+  inline const Vector3D<S> &operator-=(S s)
+  {
+    x -= s;
+    y -= s;
+    z -= s;
+    return *this;
+  }
+  //! Assign and mult operator
+  inline const Vector3D<S> &operator*=(const Vector3D<S> &v)
+  {
+    x *= v.x;
+    y *= v.y;
+    z *= v.z;
+    return *this;
+  }
+  //! Assign and mult operator
+  inline const Vector3D<S> &operator*=(S s)
+  {
+    x *= s;
+    y *= s;
+    z *= s;
+    return *this;
+  }
+  //! Assign and div operator
+  inline const Vector3D<S> &operator/=(const Vector3D<S> &v)
+  {
+    x /= v.x;
+    y /= v.y;
+    z /= v.z;
+    return *this;
+  }
+  //! Assign and div operator
+  inline const Vector3D<S> &operator/=(S s)
+  {
+    x /= s;
+    y /= s;
+    z /= s;
+    return *this;
+  }
+  //! Negation operator
+  inline Vector3D<S> operator-() const
+  {
+    return Vector3D<S>(-x, -y, -z);
+  }
+
+  //! Get smallest component
+  inline S min() const
+  {
+    return (x < y) ? ((x < z) ? x : z) : ((y < z) ? y : z);
+  }
+  //! Get biggest component
+  inline S max() const
+  {
+    return (x > y) ? ((x > z) ? x : z) : ((y > z) ? y : z);
+  }
+
+  //! Test if all components are zero
+  inline bool empty()
+  {
+    return x == 0 && y == 0 && z == 0;
+  }
+
+  //! access operator
+  inline S &operator[](unsigned int i)
+  {
+    return value[i];
+  }
+  //! constant access operator
+  inline const S &operator[](unsigned int i) const
+  {
+    return value[i];
+  }
+
+  //! debug output vector to a string
+  std::string toString() const;
+
+  //! test if nans are present
+  bool isValid() const;
+
+  //! actual values
+  union {
+    S value[3];
+    struct {
+      S x;
+      S y;
+      S z;
+    };
+    struct {
+      S X;
+      S Y;
+      S Z;
+    };
+  };
+
+  //! zero element
+  static const Vector3D<S> Zero, Invalid;
+
+  //! For compatibility with 4d vectors (discards 4th comp)
+  inline Vector3D(S vx, S vy, S vz, S vDummy) : x(vx), y(vy), z(vz)
+  {
+  }
+
+ protected:
+};
+
+//! helper to check whether float/double value is non-zero
+inline bool notZero(Real f)
+{
+  if (std::abs(f) > VECTOR_EPSILON)
+    return true;
+  return false;
+}
+
+//************************************************************************
+// Additional operators
+//************************************************************************
+
+//! Addition operator
+template<class S> inline Vector3D<S> operator+(const Vector3D<S> &v1, const Vector3D<S> &v2)
+{
+  return Vector3D<S>(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
+}
+//! Addition operator
+template<class S, class S2> inline Vector3D<S> operator+(const Vector3D<S> &v, S2 s)
+{
+  return Vector3D<S>(v.x + s, v.y + s, v.z + s);
+}
+//! Addition operator
+template<class S, class S2> inline Vector3D<S> operator+(S2 s, const Vector3D<S> &v)
+{
+  return Vector3D<S>(v.x + s, v.y + s, v.z + s);
+}
+
+//! Subtraction operator
+template<class S> inline Vector3D<S> operator-(const Vector3D<S> &v1, const Vector3D<S> &v2)
+{
+  return Vector3D<S>(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
+}
+//! Subtraction operator
+template<class S, class S2> inline Vector3D<S> operator-(const Vector3D<S> &v, S2 s)
+{
+  return Vector3D<S>(v.x - s, v.y - s, v.z - s);
+}
+//! Subtraction operator
+template<class S, class S2> inline Vector3D<S> operator-(S2 s, const Vector3D<S> &v)
+{
+  return Vector3D<S>(s - v.x, s - v.y, s - v.z);
+}
+
+//! Multiplication operator
+template<class S> inline Vector3D<S> operator*(const Vector3D<S> &v1, const Vector3D<S> &v2)
+{
+  return Vector3D<S>(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z);
+}
+//! Multiplication operator
+template<class S, class S2> inline Vector3D<S> operator*(const Vector3D<S> &v, S2 s)
+{
+  return Vector3D<S>(v.x * s, v.y * s, v.z * s);
+}
+//! Multiplication operator
+template<class S, class S2> inline Vector3D<S> operator*(S2 s, const Vector3D<S> &v)
+{
+  return Vector3D<S>(s * v.x, s * v.y, s * v.z);
+}
+
+//! Division operator
+template<class S> inline Vector3D<S> operator/(const Vector3D<S> &v1, const Vector3D<S> &v2)
+{
+  return Vector3D<S>(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z);
+}
+//! Division operator
+template<class S, class S2> inline Vector3D<S> operator/(const Vector3D<S> &v, S2 s)
+{
+  return Vector3D<S>(v.x / s, v.y / s, v.z / s);
+}
+//! Division operator
+template<class S, class S2> inline Vector3D<S> operator/(S2 s, const Vector3D<S> &v)
+{
+  return Vector3D<S>(s / v.x, s / v.y, s / v.z);
+}
+
+//! Comparison operator
+template<class S> inline bool operator==(const Vector3D<S> &s1, const Vector3D<S> &s2)
+{
+  return s1.x == s2.x && s1.y == s2.y && s1.z == s2.z;
+}
+
+//! Comparison operator
+template<class S> inline bool operator!=(const Vector3D<S> &s1, const Vector3D<S> &s2)
+{
+  return s1.x != s2.x || s1.y != s2.y || s1.z != s2.z;
+}
+
+//************************************************************************
+// External functions
+//************************************************************************
+
+//! Min operator
+template<class S> inline Vector3D<S> vmin(const Vector3D<S> &s1, const Vector3D<S> &s2)
+{
+  return Vector3D<S>(std::min(s1.x, s2.x), std::min(s1.y, s2.y), std::min(s1.z, s2.z));
+}
+
+//! Min operator
+template<class S, class S2> inline Vector3D<S> vmin(const Vector3D<S> &s1, S2 s2)
+{
+  return Vector3D<S>(std::min(s1.x, s2), std::min(s1.y, s2), std::min(s1.z, s2));
+}
+
+//! Min operator
+template<class S1, class S> inline Vector3D<S> vmin(S1 s1, const Vector3D<S> &s2)
+{
+  return Vector3D<S>(std::min(s1, s2.x), std::min(s1, s2.y), std::min(s1, s2.z));
+}
+
+//! Max operator
+template<class S> inline Vector3D<S> vmax(const Vector3D<S> &s1, const Vector3D<S> &s2)
+{
+  return Vector3D<S>(std::max(s1.x, s2.x), std::max(s1.y, s2.y), std::max(s1.z, s2.z));
+}
+
+//! Max operator
+template<class S, class S2> inline Vector3D<S> vmax(const Vector3D<S> &s1, S2 s2)
+{
+  return Vector3D<S>(std::max(s1.x, s2), std::max(s1.y, s2), std::max(s1.z, s2));
+}
+
+//! Max operator
+template<class S1, class S> inline Vector3D<S> vmax(S1 s1, const Vector3D<S> &s2)
+{
+  return Vector3D<S>(std::max(s1, s2.x), std::max(s1, s2.y), std::max(s1, s2.z));
+}
+
+//! Dot product
+template<class S> inline S dot(const Vector3D<S> &t, const Vector3D<S> &v)
+{
+  return t.x * v.x + t.y * v.y + t.z * v.z;
+}
+
+//! Cross product
+template<class S> inline Vector3D<S> cross(const Vector3D<S> &t, const Vector3D<S> &v)
+{
+  Vector3D<S> cp(
+      ((t.y * v.z) - (t.z * v.y)), ((t.z * v.x) - (t.x * v.z)), ((t.x * v.y) - (t.y * v.x)));
+  return cp;
+}
+
+//! Project a vector into a plane, defined by its normal
+/*! Projects a vector into a plane normal to the given vector, which must
+  have unit length. Self is modified.
+  \param v The vector to project
+  \param n The plane normal
+  \return The projected vector */
+template<class S>
+inline const Vector3D<S> &projectNormalTo(const Vector3D<S> &v, const Vector3D<S> &n)
+{
+  S sprod = dot(v, n);
+  return v - n * dot(v, n);
+}
+
+//! Compute the magnitude (length) of the vector
+//! (clamps to 0 and 1 with VECTOR_EPSILON)
+template<class S> inline S norm(const Vector3D<S> &v)
+{
+  S l = v.x * v.x + v.y * v.y + v.z * v.z;
+  if (l <= VECTOR_EPSILON * VECTOR_EPSILON)
+    return (0.);
+  return (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) ? 1. : sqrt(l);
+}
+
+//! Compute squared magnitude
+template<class S> inline S normSquare(const Vector3D<S> &v)
+{
+  return v.x * v.x + v.y * v.y + v.z * v.z;
+}
+
+//! compatibility, allow use of int, Real and Vec inputs with norm/normSquare
+inline Real norm(const Real v)
+{
+  return fabs(v);
+}
+inline Real normSquare(const Real v)
+{
+  return square(v);
+}
+inline Real norm(const int v)
+{
+  return abs(v);
+}
+inline Real normSquare(const int v)
+{
+  return square(v);
+}
+
+//! Returns a normalized vector
+template<class S> inline Vector3D<S> getNormalized(const Vector3D<S> &v)
+{
+  S l = v.x * v.x + v.y * v.y + v.z * v.z;
+  if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON)
+    return v; /* normalized "enough"... */
+  else if (l > VECTOR_EPSILON * VECTOR_EPSILON) {
+    S fac = 1. / sqrt(l);
+    return Vector3D<S>(v.x * fac, v.y * fac, v.z * fac);
+  }
+  else
+    return Vector3D<S>((S)0);
+}
+
+//! Compute the norm of the vector and normalize it.
+/*! \return The value of the norm */
+template<class S> inline S normalize(Vector3D<S> &v)
+{
+  S norm;
+  S l = v.x * v.x + v.y * v.y + v.z * v.z;
+  if (fabs(l - 1.) < VECTOR_EPSILON * VECTOR_EPSILON) {
+    norm = 1.;
+  }
+  else if (l > VECTOR_EPSILON * VECTOR_EPSILON) {
+    norm = sqrt(l);
+    v *= 1. / norm;
+  }
+  else {
+    v = Vector3D<S>::Zero;
+    norm = 0.;
+  }
+  return (S)norm;
+}
+
+//! Obtain an orthogonal vector
+/*! Compute a vector that is orthonormal to the given vector.
+ *  Nothing else can be assumed for the direction of the new vector.
+ *  \return The orthonormal vector */
+template<class S> Vector3D<S> getOrthogonalVector(const Vector3D<S> &v)
+{
+  // Determine the  component with max. absolute value
+  int maxIndex = (fabs(v.x) > fabs(v.y)) ? 0 : 1;
+  maxIndex = (fabs(v[maxIndex]) > fabs(v.z)) ? maxIndex : 2;
+
+  // Choose another axis than the one with max. component and project
+  // orthogonal to self
+  Vector3D<S> o(0.0);
+  o[(maxIndex + 1) % 3] = 1;
+
+  Vector3D<S> c = cross(v, o);
+  normalize(c);
+  return c;
+}
+
+//! Convert vector to polar coordinates
+/*! Stable vector to angle conversion
+ *\param v vector to convert
+  \param phi unique angle [0,2PI]
+  \param theta unique angle [0,PI]
+ */
+template<class S> inline void vecToAngle(const Vector3D<S> &v, S &phi, S &theta)
+{
+  if (fabs(v.y) < VECTOR_EPSILON)
+    theta = M_PI / 2;
+  else if (fabs(v.x) < VECTOR_EPSILON && fabs(v.z) < VECTOR_EPSILON)
+    theta = (v.y >= 0) ? 0 : M_PI;
+  else
+    theta = atan(sqrt(v.x * v.x + v.z * v.z) / v.y);
+  if (theta < 0)
+    theta += M_PI;
+
+  if (fabs(v.x) < VECTOR_EPSILON)
+    phi = M_PI / 2;
+  else
+    phi = atan(v.z / v.x);
+  if (phi < 0)
+    phi += M_PI;
+  if (fabs(v.z) < VECTOR_EPSILON)
+    phi = (v.x >= 0) ? 0 : M_PI;
+  else if (v.z < 0)
+    phi += M_PI;
+}
+
+//! Compute vector reflected at a surface
+/*! Compute a vector, that is self (as an incoming vector)
+ * reflected at a surface with a distinct normal vector.
+ * Note that the normal is reversed, if the scalar product with it is positive.
+  \param t The incoming vector
+  \param n The surface normal
+  \return The new reflected vector
+  */
+template<class S> inline Vector3D<S> reflectVector(const Vector3D<S> &t, const Vector3D<S> &n)
+{
+  Vector3D<S> nn = (dot(t, n) > 0.0) ? (n * -1.0) : n;
+  return (t - nn * (2.0 * dot(nn, t)));
+}
+
+//! Compute vector refracted at a surface
+/*! \param t The incoming vector
+ *  \param n The surface normal
+ *  \param nt The "inside" refraction index
+ *  \param nair The "outside" refraction index
+ *  \param refRefl Set to 1 on total reflection
+ *  \return The refracted vector
+ */
+template<class S>
+inline Vector3D<S> refractVector(
+    const Vector3D<S> &t, const Vector3D<S> &normal, S nt, S nair, int &refRefl)
+{
+  // from Glassner's book, section 5.2 (Heckberts method)
+  S eta = nair / nt;
+  S n = -dot(t, normal);
+  S tt = 1.0 + eta * eta * (n * n - 1.0);
+  if (tt < 0.0) {
+    // we have total reflection!
+    refRefl = 1;
+  }
+  else {
+    // normal reflection
+    tt = eta * n - sqrt(tt);
+    return (t * eta + normal * tt);
+  }
+  return t;
+}
+
+//! Outputs the object in human readable form as string
+template<class S> std::string Vector3D<S>::toString() const
+{
+  char buf[256];
+  snprintf(buf,
+           256,
+           "[%+4.6f,%+4.6f,%+4.6f]",
+           (double)(*this)[0],
+           (double)(*this)[1],
+           (double)(*this)[2]);
+  // for debugging, optionally increase precision:
+  // snprintf ( buf,256,"[%+4.16f,%+4.16f,%+4.16f]", ( double ) ( *this ) [0], ( double ) ( *this )
+  // [1], ( double ) ( *this ) [2] );
+  return std::string(buf);
+}
+
+template<> std::string Vector3D<int>::toString() const;
+
+//! Outputs the object in human readable form to stream
+/*! Output format [x,y,z] */
+template<class S> std::ostream &operator<<(std::ostream &os, const Vector3D<S> &i)
+{
+  os << i.toString();
+  return os;
+}
+
+//! Reads the contents of the object from a stream
+/*! Input format [x,y,z] */
+template<class S> std::istream &operator>>(std::istream &is, Vector3D<S> &i)
+{
+  char c;
+  char dummy[3];
+  is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> c;
+  return is;
+}
+
+/**************************************************************************/
+// Define default vector alias
+/**************************************************************************/
+
+//! 3D vector class of type Real (typically float)
+typedef Vector3D<Real> Vec3;
+
+//! 3D vector class of type int
+typedef Vector3D<int> Vec3i;
+
+//! convert to Real Vector
+template<class T> inline Vec3 toVec3(T v)
+{
+  return Vec3(v[0], v[1], v[2]);
+}
+
+//! convert to int Vector
+template<class T> inline Vec3i toVec3i(T v)
+{
+  return Vec3i((int)v[0], (int)v[1], (int)v[2]);
+}
+
+//! convert to int Vector
+template<class T> inline Vec3i toVec3i(T v0, T v1, T v2)
+{
+  return Vec3i((int)v0, (int)v1, (int)v2);
+}
+
+//! round, and convert to int Vector
+template<class T> inline Vec3i toVec3iRound(T v)
+{
+  return Vec3i((int)round(v[0]), (int)round(v[1]), (int)round(v[2]));
+}
+
+//! convert to int Vector if values are close enough to an int
+template<class T> inline Vec3i toVec3iChecked(T v)
+{
+  Vec3i ret;
+  for (size_t i = 0; i < 3; i++) {
+    Real a = v[i];
+    if (fabs(a - floor(a + 0.5)) > 1e-5)
+      errMsg("argument is not an int, cannot convert");
+    ret[i] = (int)(a + 0.5);
+  }
+  return ret;
+}
+
+//! convert to double Vector
+template<class T> inline Vector3D<double> toVec3d(T v)
+{
+  return Vector3D<double>(v[0], v[1], v[2]);
+}
+
+//! convert to float Vector
+template<class T> inline Vector3D<float> toVec3f(T v)
+{
+  return Vector3D<float>(v[0], v[1], v[2]);
+}
+
+/**************************************************************************/
+// Specializations for common math functions
+/**************************************************************************/
+
+template<> inline Vec3 clamp<Vec3>(const Vec3 &a, const Vec3 &b, const Vec3 &c)
+{
+  return Vec3(clamp(a.x, b.x, c.x), clamp(a.y, b.y, c.y), clamp(a.z, b.z, c.z));
+}
+template<> inline Vec3 safeDivide<Vec3>(const Vec3 &a, const Vec3 &b)
+{
+  return Vec3(safeDivide(a.x, b.x), safeDivide(a.y, b.y), safeDivide(a.z, b.z));
+}
+template<> inline Vec3 nmod<Vec3>(const Vec3 &a, const Vec3 &b)
+{
+  return Vec3(nmod(a.x, b.x), nmod(a.y, b.y), nmod(a.z, b.z));
+}
+
+};  // namespace Manta
+
+#endif