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

1 files changed, 2260 insertions, 0 deletions

diff --git a/extern/mantaflow/preprocessed/grid.h b/extern/mantaflow/preprocessed/grid.h
new file mode 100644
index 00000000000..bd4e0f99f85
--- /dev/null
+++ b/extern/mantaflow/preprocessed/grid.h
@@ -0,0 +1,2260 @@
+
+
+// DO NOT EDIT !
+// This file is generated using the MantaFlow preprocessor (prep generate).
+
+/******************************************************************************
+ *
+ * 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
+ *
+ * Grid representation
+ *
+ ******************************************************************************/
+
+#ifndef _GRID_H
+#define _GRID_H
+
+#include "manta.h"
+#include "vectorbase.h"
+#include "interpol.h"
+#include "interpolHigh.h"
+#include "kernel.h"
+
+namespace Manta {
+class LevelsetGrid;
+
+//! Base class for all grids
+class GridBase : public PbClass {
+ public:
+  enum GridType {
+    TypeNone = 0,
+    TypeReal = 1,
+    TypeInt = 2,
+    TypeVec3 = 4,
+    TypeMAC = 8,
+    TypeLevelset = 16,
+    TypeFlags = 32
+  };
+
+  GridBase(FluidSolver *parent);
+  static int _W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    PbClass *obj = Pb::objFromPy(_self);
+    if (obj)
+      delete obj;
+    try {
+      PbArgs _args(_linargs, _kwds);
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(0, "GridBase::GridBase", !noTiming);
+      {
+        ArgLocker _lock;
+        FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
+        obj = new GridBase(parent);
+        obj->registerObject(_self, &_args);
+        _args.check();
+      }
+      pbFinalizePlugin(obj->getParent(), "GridBase::GridBase", !noTiming);
+      return 0;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::GridBase", e.what());
+      return -1;
+    }
+  }
+
+  //! Get the grids X dimension
+  inline int getSizeX() const
+  {
+    return mSize.x;
+  }
+  static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::getSizeX", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getSizeX());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::getSizeX", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::getSizeX", e.what());
+      return 0;
+    }
+  }
+
+  //! Get the grids Y dimension
+  inline int getSizeY() const
+  {
+    return mSize.y;
+  }
+  static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::getSizeY", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getSizeY());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::getSizeY", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::getSizeY", e.what());
+      return 0;
+    }
+  }
+
+  //! Get the grids Z dimension
+  inline int getSizeZ() const
+  {
+    return mSize.z;
+  }
+  static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::getSizeZ", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getSizeZ());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::getSizeZ", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::getSizeZ", e.what());
+      return 0;
+    }
+  }
+
+  //! Get the grids dimensions
+  inline Vec3i getSize() const
+  {
+    return mSize;
+  }
+  static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::getSize", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getSize());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::getSize", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::getSize", e.what());
+      return 0;
+    }
+  }
+
+  //! Get Stride in X dimension
+  inline IndexInt getStrideX() const
+  {
+    return 1;
+  }
+  //! Get Stride in Y dimension
+  inline IndexInt getStrideY() const
+  {
+    return mSize.x;
+  }
+  //! Get Stride in Z dimension
+  inline IndexInt getStrideZ() const
+  {
+    return mStrideZ;
+  }
+
+  inline Real getDx() const
+  {
+    return mDx;
+  }
+
+  //! Check if indices are within bounds, otherwise error (should only be called when debugging)
+  inline void checkIndex(int i, int j, int k) const;
+  //! Check if indices are within bounds, otherwise error (should only be called when debugging)
+  inline void checkIndex(IndexInt idx) const;
+  //! Check if index is within given boundaries
+  inline bool isInBounds(const Vec3i &p, int bnd) const;
+  //! Check if index is within given boundaries
+  inline bool isInBounds(const Vec3i &p) const;
+  //! Check if index is within given boundaries
+  inline bool isInBounds(const Vec3 &p, int bnd = 0) const
+  {
+    return isInBounds(toVec3i(p), bnd);
+  }
+  //! Check if linear index is in the range of the array
+  inline bool isInBounds(IndexInt idx) const;
+
+  //! Get the type of grid
+  inline GridType getType() const
+  {
+    return mType;
+  }
+  //! Check dimensionality
+  inline bool is3D() const
+  {
+    return m3D;
+  }
+  static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::is3D", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->is3D());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::is3D", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::is3D", e.what());
+      return 0;
+    }
+  }
+
+  //! Get index into the data
+  inline IndexInt index(int i, int j, int k) const
+  {
+    DEBUG_ONLY(checkIndex(i, j, k));
+    return (IndexInt)i + (IndexInt)mSize.x * j + (IndexInt)mStrideZ * k;
+  }
+  //! Get index into the data
+  inline IndexInt index(const Vec3i &pos) const
+  {
+    DEBUG_ONLY(checkIndex(pos.x, pos.y, pos.z));
+    return (IndexInt)pos.x + (IndexInt)mSize.x * pos.y + (IndexInt)mStrideZ * pos.z;
+  }
+
+  //! grid4d compatibility functions
+  inline bool is4D() const
+  {
+    return false;
+  }
+  static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::is4D", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->is4D());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::is4D", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::is4D", e.what());
+      return 0;
+    }
+  }
+
+  inline int getSizeT() const
+  {
+    return 1;
+  }
+  static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::getSizeT", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getSizeT());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::getSizeT", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::getSizeT", e.what());
+      return 0;
+    }
+  }
+
+  inline int getStrideT() const
+  {
+    return 0;
+  }
+  static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      GridBase *pbo = dynamic_cast<GridBase *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "GridBase::getStrideT", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getStrideT());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "GridBase::getStrideT", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("GridBase::getStrideT", e.what());
+      return 0;
+    }
+  }
+
+  inline int index(int i, int j, int k, int unused) const
+  {
+    return index(i, j, k);
+  }
+  inline bool isInBounds(int i, int j, int k, int t, int bnd) const
+  {
+    if (t != 0)
+      return false;
+    return isInBounds(Vec3i(i, j, k), bnd);
+  }
+
+ protected:
+  GridType mType;
+  Vec3i mSize;
+  Real mDx;
+  bool m3D;  // precomputed Z shift: to ensure 2D compatibility, always use this instead of sx*sy !
+  IndexInt mStrideZ;
+ public:
+  PbArgs _args;
+}
+#define _C_GridBase
+;
+
+//! Grid class
+
+template<class T> class Grid : public GridBase {
+ public:
+  //! init new grid, values are set to zero
+  Grid(FluidSolver *parent, bool show = true);
+  static int _W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    PbClass *obj = Pb::objFromPy(_self);
+    if (obj)
+      delete obj;
+    try {
+      PbArgs _args(_linargs, _kwds);
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(0, "Grid::Grid", !noTiming);
+      {
+        ArgLocker _lock;
+        FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
+        bool show = _args.getOpt<bool>("show", 1, true, &_lock);
+        obj = new Grid(parent, show);
+        obj->registerObject(_self, &_args);
+        _args.check();
+      }
+      pbFinalizePlugin(obj->getParent(), "Grid::Grid", !noTiming);
+      return 0;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::Grid", e.what());
+      return -1;
+    }
+  }
+
+  //! init new grid with an existing array
+  Grid(FluidSolver *parent, T *data, bool show = true);
+  //! create new & copy content from another grid
+  Grid(const Grid<T> &a);
+  //! return memory to solver
+  virtual ~Grid();
+
+  typedef T BASETYPE;
+  typedef GridBase BASETYPE_GRID;
+
+  void save(std::string name);
+  static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::save", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        std::string name = _args.get<std::string>("name", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->save(name);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::save", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::save", e.what());
+      return 0;
+    }
+  }
+
+  void load(std::string name);
+  static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::load", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        std::string name = _args.get<std::string>("name", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->load(name);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::load", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::load", e.what());
+      return 0;
+    }
+  }
+
+  //! set all cells to zero
+  void clear();
+  static PyObject *_W_12(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::clear", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->clear();
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::clear", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::clear", e.what());
+      return 0;
+    }
+  }
+
+  //! all kinds of access functions, use grid(), grid[] or grid.get()
+  //! access data
+  inline T get(int i, int j, int k) const
+  {
+    return mData[index(i, j, k)];
+  }
+  //! access data
+  inline T &get(int i, int j, int k)
+  {
+    return mData[index(i, j, k)];
+  }
+  //! access data
+  inline T get(IndexInt idx) const
+  {
+    DEBUG_ONLY(checkIndex(idx));
+    return mData[idx];
+  }
+  //! access data
+  inline T get(const Vec3i &pos) const
+  {
+    return mData[index(pos)];
+  }
+  //! access data
+  inline T &operator()(int i, int j, int k)
+  {
+    return mData[index(i, j, k)];
+  }
+  //! access data
+  inline T operator()(int i, int j, int k) const
+  {
+    return mData[index(i, j, k)];
+  }
+  //! access data
+  inline T &operator()(IndexInt idx)
+  {
+    DEBUG_ONLY(checkIndex(idx));
+    return mData[idx];
+  }
+  //! access data
+  inline T operator()(IndexInt idx) const
+  {
+    DEBUG_ONLY(checkIndex(idx));
+    return mData[idx];
+  }
+  //! access data
+  inline T &operator()(const Vec3i &pos)
+  {
+    return mData[index(pos)];
+  }
+  //! access data
+  inline T operator()(const Vec3i &pos) const
+  {
+    return mData[index(pos)];
+  }
+  //! access data
+  inline T &operator[](IndexInt idx)
+  {
+    DEBUG_ONLY(checkIndex(idx));
+    return mData[idx];
+  }
+  //! access data
+  inline const T operator[](IndexInt idx) const
+  {
+    DEBUG_ONLY(checkIndex(idx));
+    return mData[idx];
+  }
+
+  // interpolated access
+  inline T getInterpolated(const Vec3 &pos) const
+  {
+    return interpol<T>(mData, mSize, mStrideZ, pos);
+  }
+  inline void setInterpolated(const Vec3 &pos, const T &val, Grid<Real> &sumBuffer) const
+  {
+    setInterpol<T>(mData, mSize, mStrideZ, pos, val, &sumBuffer[0]);
+  }
+  // higher order interpolation (1=linear, 2=cubic)
+  inline T getInterpolatedHi(const Vec3 &pos, int order) const
+  {
+    switch (order) {
+      case 1:
+        return interpol<T>(mData, mSize, mStrideZ, pos);
+      case 2:
+        return interpolCubic<T>(mData, mSize, mStrideZ, pos);
+      default:
+        assertMsg(false, "Unknown interpolation order " << order);
+    }
+    return T(0.);  // should never be reached, just to prevent compiler warnings
+  }
+
+  // assignment / copy
+
+  //! warning - do not use "=" for grids in python, this copies the reference! not the grid
+  //! content...
+  // Grid<T>& operator=(const Grid<T>& a);
+  //! copy content from other grid (use this one instead of operator= !)
+  Grid<T> &copyFrom(const Grid<T> &a, bool copyType = true);
+  static PyObject *_W_13(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::copyFrom", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        const Grid<T> &a = *_args.getPtr<Grid<T>>("a", 0, &_lock);
+        bool copyType = _args.getOpt<bool>("copyType", 1, true, &_lock);
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->copyFrom(a, copyType));
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::copyFrom", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::copyFrom", e.what());
+      return 0;
+    }
+  }
+  // old: { *this = a; }
+
+  // helper functions to work with grids in scene files
+
+  //! get grid type
+  int getGridType();
+  static PyObject *_W_14(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::getGridType", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getGridType());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::getGridType", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::getGridType", e.what());
+      return 0;
+    }
+  }
+
+  //! add/subtract other grid
+  void add(const Grid<T> &a);
+  static PyObject *_W_15(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::add", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        const Grid<T> &a = *_args.getPtr<Grid<T>>("a", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->add(a);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::add", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::add", e.what());
+      return 0;
+    }
+  }
+
+  void sub(const Grid<T> &a);
+  static PyObject *_W_16(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::sub", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        const Grid<T> &a = *_args.getPtr<Grid<T>>("a", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->sub(a);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::sub", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::sub", e.what());
+      return 0;
+    }
+  }
+
+  //! set all cells to constant value
+  void setConst(T s);
+  static PyObject *_W_17(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::setConst", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        T s = _args.get<T>("s", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->setConst(s);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::setConst", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::setConst", e.what());
+      return 0;
+    }
+  }
+
+  //! add constant to all grid cells
+  void addConst(T s);
+  static PyObject *_W_18(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::addConst", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        T s = _args.get<T>("s", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->addConst(s);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::addConst", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::addConst", e.what());
+      return 0;
+    }
+  }
+
+  //! add scaled other grid to current one (note, only "Real" factor, "T" type not supported here!)
+  void addScaled(const Grid<T> &a, const T &factor);
+  static PyObject *_W_19(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::addScaled", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        const Grid<T> &a = *_args.getPtr<Grid<T>>("a", 0, &_lock);
+        const T &factor = *_args.getPtr<T>("factor", 1, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->addScaled(a, factor);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::addScaled", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::addScaled", e.what());
+      return 0;
+    }
+  }
+
+  //! multiply contents of grid
+  void mult(const Grid<T> &a);
+  static PyObject *_W_20(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::mult", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        const Grid<T> &a = *_args.getPtr<Grid<T>>("a", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->mult(a);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::mult", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::mult", e.what());
+      return 0;
+    }
+  }
+
+  //! multiply each cell by a constant scalar value
+  void multConst(T s);
+  static PyObject *_W_21(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::multConst", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        T s = _args.get<T>("s", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->multConst(s);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::multConst", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::multConst", e.what());
+      return 0;
+    }
+  }
+
+  //! clamp content to range (for vec3, clamps each component separately)
+  void clamp(Real min, Real max);
+  static PyObject *_W_22(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::clamp", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        Real min = _args.get<Real>("min", 0, &_lock);
+        Real max = _args.get<Real>("max", 1, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->clamp(min, max);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::clamp", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::clamp", e.what());
+      return 0;
+    }
+  }
+
+  //! reduce small values to zero
+  void stomp(const T &threshold);
+  static PyObject *_W_23(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::stomp", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        const T &threshold = *_args.getPtr<T>("threshold", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->stomp(threshold);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::stomp", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::stomp", e.what());
+      return 0;
+    }
+  }
+
+  //! permute grid axes, e.g. switch y with z (0,2,1)
+  void permuteAxes(int axis0, int axis1, int axis2);
+  static PyObject *_W_24(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::permuteAxes", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int axis0 = _args.get<int>("axis0", 0, &_lock);
+        int axis1 = _args.get<int>("axis1", 1, &_lock);
+        int axis2 = _args.get<int>("axis2", 2, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->permuteAxes(axis0, axis1, axis2);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::permuteAxes", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::permuteAxes", e.what());
+      return 0;
+    }
+  }
+
+  //! permute grid axes, e.g. switch y with z (0,2,1)
+  void permuteAxesCopyToGrid(int axis0, int axis1, int axis2, Grid<T> &out);
+  static PyObject *_W_25(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::permuteAxesCopyToGrid", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int axis0 = _args.get<int>("axis0", 0, &_lock);
+        int axis1 = _args.get<int>("axis1", 1, &_lock);
+        int axis2 = _args.get<int>("axis2", 2, &_lock);
+        Grid<T> &out = *_args.getPtr<Grid<T>>("out", 3, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->permuteAxesCopyToGrid(axis0, axis1, axis2, out);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::permuteAxesCopyToGrid", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::permuteAxesCopyToGrid", e.what());
+      return 0;
+    }
+  }
+
+  // common compound operators
+  //! get absolute max value in grid
+  Real getMaxAbs() const;
+  static PyObject *_W_26(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::getMaxAbs", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getMaxAbs());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::getMaxAbs", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::getMaxAbs", e.what());
+      return 0;
+    }
+  }
+
+  //! get max value in grid
+  Real getMax() const;
+  static PyObject *_W_27(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::getMax", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getMax());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::getMax", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::getMax", e.what());
+      return 0;
+    }
+  }
+
+  //! get min value in grid
+  Real getMin() const;
+  static PyObject *_W_28(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::getMin", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getMin());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::getMin", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::getMin", e.what());
+      return 0;
+    }
+  }
+
+  //! calculate L1 norm of grid content
+  Real getL1(int bnd = 0);
+  static PyObject *_W_29(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::getL1", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int bnd = _args.getOpt<int>("bnd", 0, 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getL1(bnd));
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::getL1", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::getL1", e.what());
+      return 0;
+    }
+  }
+
+  //! calculate L2 norm of grid content
+  Real getL2(int bnd = 0);
+  static PyObject *_W_30(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::getL2", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int bnd = _args.getOpt<int>("bnd", 0, 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getL2(bnd));
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::getL2", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::getL2", e.what());
+      return 0;
+    }
+  }
+
+  //! set all boundary cells to constant value (Dirichlet)
+  void setBound(T value, int boundaryWidth = 1);
+  static PyObject *_W_31(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::setBound", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        T value = _args.get<T>("value", 0, &_lock);
+        int boundaryWidth = _args.getOpt<int>("boundaryWidth", 1, 1, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->setBound(value, boundaryWidth);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::setBound", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::setBound", e.what());
+      return 0;
+    }
+  }
+
+  //! set all boundary cells to last inner value (Neumann)
+  void setBoundNeumann(int boundaryWidth = 1);
+  static PyObject *_W_32(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::setBoundNeumann", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int boundaryWidth = _args.getOpt<int>("boundaryWidth", 0, 1, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->setBoundNeumann(boundaryWidth);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::setBoundNeumann", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::setBoundNeumann", e.what());
+      return 0;
+    }
+  }
+
+  //! get data pointer of grid
+  std::string getDataPointer();
+  static PyObject *_W_33(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::getDataPointer", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->getDataPointer());
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::getDataPointer", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::getDataPointer", e.what());
+      return 0;
+    }
+  }
+
+  //! debugging helper, print grid from python. skip boundary of width bnd
+  void printGrid(int zSlice = -1, bool printIndex = false, int bnd = 1);
+  static PyObject *_W_34(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      Grid *pbo = dynamic_cast<Grid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "Grid::printGrid", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int zSlice = _args.getOpt<int>("zSlice", 0, -1, &_lock);
+        bool printIndex = _args.getOpt<bool>("printIndex", 1, false, &_lock);
+        int bnd = _args.getOpt<int>("bnd", 2, 1, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->printGrid(zSlice, printIndex, bnd);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "Grid::printGrid", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("Grid::printGrid", e.what());
+      return 0;
+    }
+  }
+
+  // c++ only operators
+  template<class S> Grid<T> &operator+=(const Grid<S> &a);
+  template<class S> Grid<T> &operator+=(const S &a);
+  template<class S> Grid<T> &operator-=(const Grid<S> &a);
+  template<class S> Grid<T> &operator-=(const S &a);
+  template<class S> Grid<T> &operator*=(const Grid<S> &a);
+  template<class S> Grid<T> &operator*=(const S &a);
+  template<class S> Grid<T> &operator/=(const Grid<S> &a);
+  template<class S> Grid<T> &operator/=(const S &a);
+  Grid<T> &safeDivide(const Grid<T> &a);
+
+  //! Swap data with another grid (no actual data is moved)
+  void swap(Grid<T> &other);
+
+  //! grid4d compatibility functions
+  inline T &operator()(int i, int j, int k, int unused)
+  {
+    return mData[index(i, j, k)];
+  }
+  inline T operator()(int i, int j, int k, int unused) const
+  {
+    return mData[index(i, j, k)];
+  }
+
+ protected:
+  T *mData;
+  bool externalData;  // True if mData is managed outside of the Fluidsolver
+ public:
+  PbArgs _args;
+}
+#define _C_Grid
+;
+
+// Python doesn't know about templates: explicit aliases needed
+
+//! Special function for staggered grids
+class MACGrid : public Grid<Vec3> {
+ public:
+  MACGrid(FluidSolver *parent, bool show = true) : Grid<Vec3>(parent, show)
+  {
+    mType = (GridType)(TypeMAC | TypeVec3);
+  }
+  static int _W_35(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    PbClass *obj = Pb::objFromPy(_self);
+    if (obj)
+      delete obj;
+    try {
+      PbArgs _args(_linargs, _kwds);
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(0, "MACGrid::MACGrid", !noTiming);
+      {
+        ArgLocker _lock;
+        FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
+        bool show = _args.getOpt<bool>("show", 1, true, &_lock);
+        obj = new MACGrid(parent, show);
+        obj->registerObject(_self, &_args);
+        _args.check();
+      }
+      pbFinalizePlugin(obj->getParent(), "MACGrid::MACGrid", !noTiming);
+      return 0;
+    }
+    catch (std::exception &e) {
+      pbSetError("MACGrid::MACGrid", e.what());
+      return -1;
+    }
+  }
+
+  MACGrid(FluidSolver *parent, Vec3 *data, bool show = true) : Grid<Vec3>(parent, data, show)
+  {
+    mType = (GridType)(TypeMAC | TypeVec3);
+  }
+
+  // specialized functions for interpolating MAC information
+  inline Vec3 getCentered(int i, int j, int k) const;
+  inline Vec3 getCentered(const Vec3i &pos) const
+  {
+    return getCentered(pos.x, pos.y, pos.z);
+  }
+  inline Vec3 getAtMACX(int i, int j, int k) const;
+  inline Vec3 getAtMACY(int i, int j, int k) const;
+  inline Vec3 getAtMACZ(int i, int j, int k) const;
+  // interpolation
+  inline Vec3 getInterpolated(const Vec3 &pos) const
+  {
+    return interpolMAC(mData, mSize, mStrideZ, pos);
+  }
+  inline void setInterpolated(const Vec3 &pos, const Vec3 &val, Vec3 *tmp) const
+  {
+    return setInterpolMAC(mData, mSize, mStrideZ, pos, val, tmp);
+  }
+  inline Vec3 getInterpolatedHi(const Vec3 &pos, int order) const
+  {
+    switch (order) {
+      case 1:
+        return interpolMAC(mData, mSize, mStrideZ, pos);
+      case 2:
+        return interpolCubicMAC(mData, mSize, mStrideZ, pos);
+      default:
+        assertMsg(false, "Unknown interpolation order " << order);
+    }
+    return Vec3(0.);  // should never be reached, just to prevent compiler warnings
+  }
+  // specials for mac grid:
+  template<int comp> inline Real getInterpolatedComponent(Vec3 pos) const
+  {
+    return interpolComponent<comp>(mData, mSize, mStrideZ, pos);
+  }
+  template<int comp> inline Real getInterpolatedComponentHi(const Vec3 &pos, int order) const
+  {
+    switch (order) {
+      case 1:
+        return interpolComponent<comp>(mData, mSize, mStrideZ, pos);
+      case 2:
+        return interpolCubicMAC(mData, mSize, mStrideZ, pos)[comp];  // warning - not yet optimized
+      default:
+        assertMsg(false, "Unknown interpolation order " << order);
+    }
+    return 0.;  // should never be reached, just to prevent compiler warnings
+  }
+
+  //! set all boundary cells of a MAC grid to certain value (Dirchlet). Respects staggered grid
+  //! locations optionally, only set normal components
+  void setBoundMAC(Vec3 value, int boundaryWidth, bool normalOnly = false);
+  static PyObject *_W_36(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      MACGrid *pbo = dynamic_cast<MACGrid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "MACGrid::setBoundMAC", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        Vec3 value = _args.get<Vec3>("value", 0, &_lock);
+        int boundaryWidth = _args.get<int>("boundaryWidth", 1, &_lock);
+        bool normalOnly = _args.getOpt<bool>("normalOnly", 2, false, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->setBoundMAC(value, boundaryWidth, normalOnly);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "MACGrid::setBoundMAC", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("MACGrid::setBoundMAC", e.what());
+      return 0;
+    }
+  }
+
+ protected:
+ public:
+  PbArgs _args;
+}
+#define _C_MACGrid
+;
+
+//! Special functions for FlagGrid
+class FlagGrid : public Grid<int> {
+ public:
+  FlagGrid(FluidSolver *parent, int dim = 3, bool show = true) : Grid<int>(parent, show)
+  {
+    mType = (GridType)(TypeFlags | TypeInt);
+  }
+  static int _W_37(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    PbClass *obj = Pb::objFromPy(_self);
+    if (obj)
+      delete obj;
+    try {
+      PbArgs _args(_linargs, _kwds);
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(0, "FlagGrid::FlagGrid", !noTiming);
+      {
+        ArgLocker _lock;
+        FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
+        int dim = _args.getOpt<int>("dim", 1, 3, &_lock);
+        bool show = _args.getOpt<bool>("show", 2, true, &_lock);
+        obj = new FlagGrid(parent, dim, show);
+        obj->registerObject(_self, &_args);
+        _args.check();
+      }
+      pbFinalizePlugin(obj->getParent(), "FlagGrid::FlagGrid", !noTiming);
+      return 0;
+    }
+    catch (std::exception &e) {
+      pbSetError("FlagGrid::FlagGrid", e.what());
+      return -1;
+    }
+  }
+
+  FlagGrid(FluidSolver *parent, int *data, int dim = 3, bool show = true)
+      : Grid<int>(parent, data, show)
+  {
+    mType = (GridType)(TypeFlags | TypeInt);
+  }
+
+  //! types of cells, in/outflow can be combined, e.g., TypeFluid|TypeInflow
+  enum CellType {
+    TypeNone = 0,
+    TypeFluid = 1,
+    TypeObstacle = 2,
+    TypeEmpty = 4,
+    TypeInflow = 8,
+    TypeOutflow = 16,
+    TypeOpen = 32,
+    TypeStick = 64,
+    TypeReserved = 256
+  };
+
+  //! access for particles
+  inline int getAt(const Vec3 &pos) const
+  {
+    return mData[index((int)pos.x, (int)pos.y, (int)pos.z)];
+  }
+
+  //! check for different flag types
+  inline bool isObstacle(IndexInt idx) const
+  {
+    return get(idx) & TypeObstacle;
+  }
+  inline bool isObstacle(int i, int j, int k) const
+  {
+    return get(i, j, k) & TypeObstacle;
+  }
+  inline bool isObstacle(const Vec3i &pos) const
+  {
+    return get(pos) & TypeObstacle;
+  }
+  inline bool isObstacle(const Vec3 &pos) const
+  {
+    return getAt(pos) & TypeObstacle;
+  }
+  inline bool isFluid(IndexInt idx) const
+  {
+    return get(idx) & TypeFluid;
+  }
+  inline bool isFluid(int i, int j, int k) const
+  {
+    return get(i, j, k) & TypeFluid;
+  }
+  inline bool isFluid(const Vec3i &pos) const
+  {
+    return get(pos) & TypeFluid;
+  }
+  inline bool isFluid(const Vec3 &pos) const
+  {
+    return getAt(pos) & TypeFluid;
+  }
+  inline bool isInflow(IndexInt idx) const
+  {
+    return get(idx) & TypeInflow;
+  }
+  inline bool isInflow(int i, int j, int k) const
+  {
+    return get(i, j, k) & TypeInflow;
+  }
+  inline bool isInflow(const Vec3i &pos) const
+  {
+    return get(pos) & TypeInflow;
+  }
+  inline bool isInflow(const Vec3 &pos) const
+  {
+    return getAt(pos) & TypeInflow;
+  }
+  inline bool isEmpty(IndexInt idx) const
+  {
+    return get(idx) & TypeEmpty;
+  }
+  inline bool isEmpty(int i, int j, int k) const
+  {
+    return get(i, j, k) & TypeEmpty;
+  }
+  inline bool isEmpty(const Vec3i &pos) const
+  {
+    return get(pos) & TypeEmpty;
+  }
+  inline bool isEmpty(const Vec3 &pos) const
+  {
+    return getAt(pos) & TypeEmpty;
+  }
+  inline bool isOutflow(IndexInt idx) const
+  {
+    return get(idx) & TypeOutflow;
+  }
+  inline bool isOutflow(int i, int j, int k) const
+  {
+    return get(i, j, k) & TypeOutflow;
+  }
+  inline bool isOutflow(const Vec3i &pos) const
+  {
+    return get(pos) & TypeOutflow;
+  }
+  inline bool isOutflow(const Vec3 &pos) const
+  {
+    return getAt(pos) & TypeOutflow;
+  }
+  inline bool isOpen(IndexInt idx) const
+  {
+    return get(idx) & TypeOpen;
+  }
+  inline bool isOpen(int i, int j, int k) const
+  {
+    return get(i, j, k) & TypeOpen;
+  }
+  inline bool isOpen(const Vec3i &pos) const
+  {
+    return get(pos) & TypeOpen;
+  }
+  inline bool isOpen(const Vec3 &pos) const
+  {
+    return getAt(pos) & TypeOpen;
+  }
+  inline bool isStick(IndexInt idx) const
+  {
+    return get(idx) & TypeStick;
+  }
+  inline bool isStick(int i, int j, int k) const
+  {
+    return get(i, j, k) & TypeStick;
+  }
+  inline bool isStick(const Vec3i &pos) const
+  {
+    return get(pos) & TypeStick;
+  }
+  inline bool isStick(const Vec3 &pos) const
+  {
+    return getAt(pos) & TypeStick;
+  }
+
+  void InitMinXWall(const int &boundaryWidth, Grid<Real> &phiWalls);
+  void InitMaxXWall(const int &boundaryWidth, Grid<Real> &phiWalls);
+  void InitMinYWall(const int &boundaryWidth, Grid<Real> &phiWalls);
+  void InitMaxYWall(const int &boundaryWidth, Grid<Real> &phiWalls);
+  void InitMinZWall(const int &boundaryWidth, Grid<Real> &phiWalls);
+  void InitMaxZWall(const int &boundaryWidth, Grid<Real> &phiWalls);
+  // Python callables
+
+  void initDomain(const int &boundaryWidth = 0,
+                  const std::string &wall = "xXyYzZ",
+                  const std::string &open = "      ",
+                  const std::string &inflow = "      ",
+                  const std::string &outflow = "      ",
+                  Grid<Real> *phiWalls = 0x00);
+  static PyObject *_W_38(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      FlagGrid *pbo = dynamic_cast<FlagGrid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "FlagGrid::initDomain", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        const int &boundaryWidth = _args.getOpt<int>("boundaryWidth", 0, 0, &_lock);
+        const std::string &wall = _args.getOpt<std::string>("wall", 1, "xXyYzZ", &_lock);
+        const std::string &open = _args.getOpt<std::string>("open", 2, "      ", &_lock);
+        const std::string &inflow = _args.getOpt<std::string>("inflow", 3, "      ", &_lock);
+        const std::string &outflow = _args.getOpt<std::string>("outflow", 4, "      ", &_lock);
+        Grid<Real> *phiWalls = _args.getPtrOpt<Grid<Real>>("phiWalls", 5, 0x00, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->initDomain(boundaryWidth, wall, open, inflow, outflow, phiWalls);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "FlagGrid::initDomain", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("FlagGrid::initDomain", e.what());
+      return 0;
+    }
+  }
+
+  void initBoundaries(const int &boundaryWidth, const int *types);
+
+  //! set fluid flags inside levelset (liquids)
+  void updateFromLevelset(LevelsetGrid &levelset);
+  static PyObject *_W_39(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      FlagGrid *pbo = dynamic_cast<FlagGrid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "FlagGrid::updateFromLevelset", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        LevelsetGrid &levelset = *_args.getPtr<LevelsetGrid>("levelset", 0, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->updateFromLevelset(levelset);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "FlagGrid::updateFromLevelset", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("FlagGrid::updateFromLevelset", e.what());
+      return 0;
+    }
+  }
+
+  //! set all cells (except obs/in/outflow) to type (fluid by default)
+  void fillGrid(int type = TypeFluid);
+  static PyObject *_W_40(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      FlagGrid *pbo = dynamic_cast<FlagGrid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "FlagGrid::fillGrid", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int type = _args.getOpt<int>("type", 0, TypeFluid, &_lock);
+        pbo->_args.copy(_args);
+        _retval = getPyNone();
+        pbo->fillGrid(type);
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "FlagGrid::fillGrid", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("FlagGrid::fillGrid", e.what());
+      return 0;
+    }
+  }
+
+  //! count no. of cells matching flags via "AND"
+  //! warning for large grids! only regular int returned (due to python interface)
+  //! optionally creates mask in RealGrid (1 where flag matches, 0 otherwise)
+  int countCells(int flag, int bnd = 0, Grid<Real> *mask = NULL);
+  static PyObject *_W_41(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+  {
+    try {
+      PbArgs _args(_linargs, _kwds);
+      FlagGrid *pbo = dynamic_cast<FlagGrid *>(Pb::objFromPy(_self));
+      bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+      pbPreparePlugin(pbo->getParent(), "FlagGrid::countCells", !noTiming);
+      PyObject *_retval = 0;
+      {
+        ArgLocker _lock;
+        int flag = _args.get<int>("flag", 0, &_lock);
+        int bnd = _args.getOpt<int>("bnd", 1, 0, &_lock);
+        Grid<Real> *mask = _args.getPtrOpt<Grid<Real>>("mask", 2, NULL, &_lock);
+        pbo->_args.copy(_args);
+        _retval = toPy(pbo->countCells(flag, bnd, mask));
+        pbo->_args.check();
+      }
+      pbFinalizePlugin(pbo->getParent(), "FlagGrid::countCells", !noTiming);
+      return _retval;
+    }
+    catch (std::exception &e) {
+      pbSetError("FlagGrid::countCells", e.what());
+      return 0;
+    }
+  }
+
+ public:
+  PbArgs _args;
+}
+#define _C_FlagGrid
+;
+
+//! helper to compute grid conversion factor between local coordinates of two grids
+inline Vec3 calcGridSizeFactor(Vec3i s1, Vec3i s2)
+{
+  return Vec3(Real(s1[0]) / s2[0], Real(s1[1]) / s2[1], Real(s1[2]) / s2[2]);
+}
+
+// prototypes for grid plugins
+void copyMacToVec3(MACGrid &source, Grid<Vec3> &target);
+void convertMacToVec3(MACGrid &source, Grid<Vec3> &target);
+void resampleVec3ToMac(Grid<Vec3> &source, MACGrid &target);
+void resampleMacToVec3(MACGrid &source, Grid<Vec3> &target);
+
+void getComponent(const Grid<Vec3> &source, Grid<Real> &target, int component);
+void setComponent(const Grid<Real> &source, Grid<Vec3> &target, int component);
+
+//******************************************************************************
+// Implementation of inline functions
+
+inline void GridBase::checkIndex(int i, int j, int k) const
+{
+  if (i < 0 || j < 0 || k < 0 || i >= mSize.x || j >= mSize.y || k >= mSize.z) {
+    std::ostringstream s;
+    s << "Grid " << mName << " dim " << mSize << " : index " << i << "," << j << "," << k
+      << " out of bound ";
+    errMsg(s.str());
+  }
+}
+
+inline void GridBase::checkIndex(IndexInt idx) const
+{
+  if (idx < 0 || idx >= mSize.x * mSize.y * mSize.z) {
+    std::ostringstream s;
+    s << "Grid " << mName << " dim " << mSize << " : index " << idx << " out of bound ";
+    errMsg(s.str());
+  }
+}
+
+bool GridBase::isInBounds(const Vec3i &p) const
+{
+  return (p.x >= 0 && p.y >= 0 && p.z >= 0 && p.x < mSize.x && p.y < mSize.y && p.z < mSize.z);
+}
+
+bool GridBase::isInBounds(const Vec3i &p, int bnd) const
+{
+  bool ret = (p.x >= bnd && p.y >= bnd && p.x < mSize.x - bnd && p.y < mSize.y - bnd);
+  if (this->is3D()) {
+    ret &= (p.z >= bnd && p.z < mSize.z - bnd);
+  }
+  else {
+    ret &= (p.z == 0);
+  }
+  return ret;
+}
+//! Check if linear index is in the range of the array
+bool GridBase::isInBounds(IndexInt idx) const
+{
+  if (idx < 0 || idx >= mSize.x * mSize.y * mSize.z) {
+    return false;
+  }
+  return true;
+}
+
+inline Vec3 MACGrid::getCentered(int i, int j, int k) const
+{
+  DEBUG_ONLY(checkIndex(i + 1, j + 1, k));
+  const IndexInt idx = index(i, j, k);
+  Vec3 v = Vec3(
+      0.5 * (mData[idx].x + mData[idx + 1].x), 0.5 * (mData[idx].y + mData[idx + mSize.x].y), 0.);
+  if (this->is3D()) {
+    DEBUG_ONLY(checkIndex(idx + mStrideZ));
+    v[2] = 0.5 * (mData[idx].z + mData[idx + mStrideZ].z);
+  }
+  return v;
+}
+
+inline Vec3 MACGrid::getAtMACX(int i, int j, int k) const
+{
+  DEBUG_ONLY(checkIndex(i - 1, j + 1, k));
+  const IndexInt idx = index(i, j, k);
+  Vec3 v = Vec3((mData[idx].x),
+                0.25 * (mData[idx].y + mData[idx - 1].y + mData[idx + mSize.x].y +
+                        mData[idx + mSize.x - 1].y),
+                0.);
+  if (this->is3D()) {
+    DEBUG_ONLY(checkIndex(idx + mStrideZ - 1));
+    v[2] = 0.25 * (mData[idx].z + mData[idx - 1].z + mData[idx + mStrideZ].z +
+                   mData[idx + mStrideZ - 1].z);
+  }
+  return v;
+}
+
+inline Vec3 MACGrid::getAtMACY(int i, int j, int k) const
+{
+  DEBUG_ONLY(checkIndex(i + 1, j - 1, k));
+  const IndexInt idx = index(i, j, k);
+  Vec3 v = Vec3(0.25 * (mData[idx].x + mData[idx - mSize.x].x + mData[idx + 1].x +
+                        mData[idx + 1 - mSize.x].x),
+                (mData[idx].y),
+                0.);
+  if (this->is3D()) {
+    DEBUG_ONLY(checkIndex(idx + mStrideZ - mSize.x));
+    v[2] = 0.25 * (mData[idx].z + mData[idx - mSize.x].z + mData[idx + mStrideZ].z +
+                   mData[idx + mStrideZ - mSize.x].z);
+  }
+  return v;
+}
+
+inline Vec3 MACGrid::getAtMACZ(int i, int j, int k) const
+{
+  const IndexInt idx = index(i, j, k);
+  DEBUG_ONLY(checkIndex(idx - mStrideZ));
+  DEBUG_ONLY(checkIndex(idx + mSize.x - mStrideZ));
+  Vec3 v = Vec3(0.25 * (mData[idx].x + mData[idx - mStrideZ].x + mData[idx + 1].x +
+                        mData[idx + 1 - mStrideZ].x),
+                0.25 * (mData[idx].y + mData[idx - mStrideZ].y + mData[idx + mSize.x].y +
+                        mData[idx + mSize.x - mStrideZ].y),
+                (mData[idx].z));
+  return v;
+}
+
+template<class T, class S> struct gridAdd : public KernelBase {
+  gridAdd(Grid<T> &me, const Grid<S> &other) : KernelBase(&me, 0), me(me), other(other)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &me, const Grid<S> &other) const
+  {
+    me[idx] += other[idx];
+  }
+  inline Grid<T> &getArg0()
+  {
+    return me;
+  }
+  typedef Grid<T> type0;
+  inline const Grid<S> &getArg1()
+  {
+    return other;
+  }
+  typedef Grid<S> type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridAdd ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, me, other);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &me;
+  const Grid<S> &other;
+};
+template<class T, class S> struct gridSub : public KernelBase {
+  gridSub(Grid<T> &me, const Grid<S> &other) : KernelBase(&me, 0), me(me), other(other)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &me, const Grid<S> &other) const
+  {
+    me[idx] -= other[idx];
+  }
+  inline Grid<T> &getArg0()
+  {
+    return me;
+  }
+  typedef Grid<T> type0;
+  inline const Grid<S> &getArg1()
+  {
+    return other;
+  }
+  typedef Grid<S> type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridSub ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, me, other);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &me;
+  const Grid<S> &other;
+};
+template<class T, class S> struct gridMult : public KernelBase {
+  gridMult(Grid<T> &me, const Grid<S> &other) : KernelBase(&me, 0), me(me), other(other)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &me, const Grid<S> &other) const
+  {
+    me[idx] *= other[idx];
+  }
+  inline Grid<T> &getArg0()
+  {
+    return me;
+  }
+  typedef Grid<T> type0;
+  inline const Grid<S> &getArg1()
+  {
+    return other;
+  }
+  typedef Grid<S> type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridMult ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, me, other);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &me;
+  const Grid<S> &other;
+};
+template<class T, class S> struct gridDiv : public KernelBase {
+  gridDiv(Grid<T> &me, const Grid<S> &other) : KernelBase(&me, 0), me(me), other(other)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &me, const Grid<S> &other) const
+  {
+    me[idx] /= other[idx];
+  }
+  inline Grid<T> &getArg0()
+  {
+    return me;
+  }
+  typedef Grid<T> type0;
+  inline const Grid<S> &getArg1()
+  {
+    return other;
+  }
+  typedef Grid<S> type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridDiv ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, me, other);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &me;
+  const Grid<S> &other;
+};
+template<class T, class S> struct gridAddScalar : public KernelBase {
+  gridAddScalar(Grid<T> &me, const S &other) : KernelBase(&me, 0), me(me), other(other)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &me, const S &other) const
+  {
+    me[idx] += other;
+  }
+  inline Grid<T> &getArg0()
+  {
+    return me;
+  }
+  typedef Grid<T> type0;
+  inline const S &getArg1()
+  {
+    return other;
+  }
+  typedef S type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridAddScalar ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, me, other);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &me;
+  const S &other;
+};
+template<class T, class S> struct gridMultScalar : public KernelBase {
+  gridMultScalar(Grid<T> &me, const S &other) : KernelBase(&me, 0), me(me), other(other)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &me, const S &other) const
+  {
+    me[idx] *= other;
+  }
+  inline Grid<T> &getArg0()
+  {
+    return me;
+  }
+  typedef Grid<T> type0;
+  inline const S &getArg1()
+  {
+    return other;
+  }
+  typedef S type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridMultScalar ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, me, other);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &me;
+  const S &other;
+};
+template<class T, class S> struct gridScaledAdd : public KernelBase {
+  gridScaledAdd(Grid<T> &me, const Grid<T> &other, const S &factor)
+      : KernelBase(&me, 0), me(me), other(other), factor(factor)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &me, const Grid<T> &other, const S &factor) const
+  {
+    me[idx] += factor * other[idx];
+  }
+  inline Grid<T> &getArg0()
+  {
+    return me;
+  }
+  typedef Grid<T> type0;
+  inline const Grid<T> &getArg1()
+  {
+    return other;
+  }
+  typedef Grid<T> type1;
+  inline const S &getArg2()
+  {
+    return factor;
+  }
+  typedef S type2;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridScaledAdd ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, me, other, factor);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &me;
+  const Grid<T> &other;
+  const S &factor;
+};
+
+template<class T> struct gridSetConst : public KernelBase {
+  gridSetConst(Grid<T> &grid, T value) : KernelBase(&grid, 0), grid(grid), value(value)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(IndexInt idx, Grid<T> &grid, T value) const
+  {
+    grid[idx] = value;
+  }
+  inline Grid<T> &getArg0()
+  {
+    return grid;
+  }
+  typedef Grid<T> type0;
+  inline T &getArg1()
+  {
+    return value;
+  }
+  typedef T type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel gridSetConst ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, grid, value);
+  }
+  void run()
+  {
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  Grid<T> &grid;
+  T value;
+};
+
+template<class T> template<class S> Grid<T> &Grid<T>::operator+=(const Grid<S> &a)
+{
+  gridAdd<T, S>(*this, a);
+  return *this;
+}
+template<class T> template<class S> Grid<T> &Grid<T>::operator+=(const S &a)
+{
+  gridAddScalar<T, S>(*this, a);
+  return *this;
+}
+template<class T> template<class S> Grid<T> &Grid<T>::operator-=(const Grid<S> &a)
+{
+  gridSub<T, S>(*this, a);
+  return *this;
+}
+template<class T> template<class S> Grid<T> &Grid<T>::operator-=(const S &a)
+{
+  gridAddScalar<T, S>(*this, -a);
+  return *this;
+}
+template<class T> template<class S> Grid<T> &Grid<T>::operator*=(const Grid<S> &a)
+{
+  gridMult<T, S>(*this, a);
+  return *this;
+}
+template<class T> template<class S> Grid<T> &Grid<T>::operator*=(const S &a)
+{
+  gridMultScalar<T, S>(*this, a);
+  return *this;
+}
+template<class T> template<class S> Grid<T> &Grid<T>::operator/=(const Grid<S> &a)
+{
+  gridDiv<T, S>(*this, a);
+  return *this;
+}
+template<class T> template<class S> Grid<T> &Grid<T>::operator/=(const S &a)
+{
+  S rez((S)1.0 / a);
+  gridMultScalar<T, S>(*this, rez);
+  return *this;
+}
+
+//******************************************************************************
+// Other helper functions
+
+// compute gradient of a scalar grid
+inline Vec3 getGradient(const Grid<Real> &data, int i, int j, int k)
+{
+  Vec3 v;
+
+  if (i > data.getSizeX() - 2)
+    i = data.getSizeX() - 2;
+  if (j > data.getSizeY() - 2)
+    j = data.getSizeY() - 2;
+  if (i < 1)
+    i = 1;
+  if (j < 1)
+    j = 1;
+  v = Vec3(data(i + 1, j, k) - data(i - 1, j, k), data(i, j + 1, k) - data(i, j - 1, k), 0.);
+
+  if (data.is3D()) {
+    if (k > data.getSizeZ() - 2)
+      k = data.getSizeZ() - 2;
+    if (k < 1)
+      k = 1;
+    v[2] = data(i, j, k + 1) - data(i, j, k - 1);
+  }
+
+  return v;
+}
+
+// interpolate grid from one size to another size
+
+template<class S> struct knInterpolateGridTempl : public KernelBase {
+  knInterpolateGridTempl(Grid<S> &target,
+                         const Grid<S> &source,
+                         const Vec3 &sourceFactor,
+                         Vec3 offset,
+                         int orderSpace = 1)
+      : KernelBase(&target, 0),
+        target(target),
+        source(source),
+        sourceFactor(sourceFactor),
+        offset(offset),
+        orderSpace(orderSpace)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(int i,
+                 int j,
+                 int k,
+                 Grid<S> &target,
+                 const Grid<S> &source,
+                 const Vec3 &sourceFactor,
+                 Vec3 offset,
+                 int orderSpace = 1) const
+  {
+    Vec3 pos = Vec3(i, j, k) * sourceFactor + offset;
+    if (!source.is3D())
+      pos[2] = 0;  // allow 2d -> 3d
+    target(i, j, k) = source.getInterpolatedHi(pos, orderSpace);
+  }
+  inline Grid<S> &getArg0()
+  {
+    return target;
+  }
+  typedef Grid<S> type0;
+  inline const Grid<S> &getArg1()
+  {
+    return source;
+  }
+  typedef Grid<S> type1;
+  inline const Vec3 &getArg2()
+  {
+    return sourceFactor;
+  }
+  typedef Vec3 type2;
+  inline Vec3 &getArg3()
+  {
+    return offset;
+  }
+  typedef Vec3 type3;
+  inline int &getArg4()
+  {
+    return orderSpace;
+  }
+  typedef int type4;
+  void runMessage()
+  {
+    debMsg("Executing kernel knInterpolateGridTempl ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    const int _maxX = maxX;
+    const int _maxY = maxY;
+    if (maxZ > 1) {
+      for (int k = __r.begin(); k != (int)__r.end(); k++)
+        for (int j = 0; j < _maxY; j++)
+          for (int i = 0; i < _maxX; i++)
+            op(i, j, k, target, source, sourceFactor, offset, orderSpace);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 0; i < _maxX; i++)
+          op(i, j, k, target, source, sourceFactor, offset, orderSpace);
+    }
+  }
+  void run()
+  {
+    if (maxZ > 1)
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
+  }
+  Grid<S> &target;
+  const Grid<S> &source;
+  const Vec3 &sourceFactor;
+  Vec3 offset;
+  int orderSpace;
+};
+// template glue code - choose interpolation based on template arguments
+template<class GRID> void interpolGridTempl(GRID &target, GRID &source)
+{
+  errMsg("interpolGridTempl - Only valid for specific instantiations");
+}
+
+}  // namespace Manta
+#endif