Merge branch 'master' into codesign_error_tracker

16 files changed, 550 insertions, 444 deletions

diff --git a/extern/mantaflow/preprocessed/plugin/advection.cpp b/extern/mantaflow/preprocessed/plugin/advection.cpp
index 13f53140348..4ccf33cc362 100644
--- a/extern/mantaflow/preprocessed/plugin/advection.cpp
+++ b/extern/mantaflow/preprocessed/plugin/advection.cpp
@@ -1336,7 +1336,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "resetPhiInObs", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -1476,7 +1476,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "advectSemiLagrange", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid *flags = _args.getPtr<FlagGrid>("flags", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/apic.cpp b/extern/mantaflow/preprocessed/plugin/apic.cpp
index 6ff893014c9..9df6b597e05 100644
--- a/extern/mantaflow/preprocessed/plugin/apic.cpp
+++ b/extern/mantaflow/preprocessed/plugin/apic.cpp
@@ -6,7 +6,7 @@
 // ----------------------------------------------------------------------------
 //
 // MantaFlow fluid solver framework
-// Copyright 2016-2017 Kiwon Um, Nils Thuerey
+// Copyright 2016-2020 Kiwon Um, Nils Thuerey
 //
 // This program is free software, distributed under the terms of the
 // Apache License, Version 2.0
@@ -21,6 +21,52 @@
 
 namespace Manta {
 
+#define FOR_INT_IJK(num) \
+  for (int i = 0; i < num; ++i) \
+    for (int j = 0; j < num; ++j) \
+      for (int k = 0; k < num; ++k)
+
+static inline IndexInt indexUFace(const Vec3 &pos, const MACGrid &ref)
+{
+  const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+  const IndexInt index = f.x * ref.getStrideX() + c.y * ref.getStrideY() + c.z * ref.getStrideZ();
+  assertDeb(ref.isInBounds(index),
+            "U face index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+                                                                 << pos.z << "]");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexVFace(const Vec3 &pos, const MACGrid &ref)
+{
+  const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+  const IndexInt index = c.x * ref.getStrideX() + f.y * ref.getStrideY() + c.z * ref.getStrideZ();
+  assertDeb(ref.isInBounds(index),
+            "V face index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+                                                                 << pos.z << "]");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexWFace(const Vec3 &pos, const MACGrid &ref)
+{
+  const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+  const IndexInt index = c.x * ref.getStrideX() + c.y * ref.getStrideY() + f.z * ref.getStrideZ();
+  assertDeb(ref.isInBounds(index),
+            "W face index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+                                                                 << pos.z << "]");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexOffset(
+    const IndexInt gidx, const int i, const int j, const int k, const MACGrid &ref)
+{
+  const IndexInt dX[2] = {0, ref.getStrideX()};
+  const IndexInt dY[2] = {0, ref.getStrideY()};
+  const IndexInt dZ[2] = {0, ref.getStrideZ()};
+  const IndexInt index = gidx + dX[i] + dY[j] + dZ[k];
+  assertDeb(ref.isInBounds(index), "Offset index " << index << " is out of bounds");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
 struct knApicMapLinearVec3ToMACGrid : public KernelBase {
   knApicMapLinearVec3ToMACGrid(const BasicParticleSystem &p,
                                MACGrid &mg,
@@ -30,7 +76,8 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
                                const ParticleDataImpl<Vec3> &cpy,
                                const ParticleDataImpl<Vec3> &cpz,
                                const ParticleDataImpl<int> *ptype,
-                               const int exclude)
+                               const int exclude,
+                               const int boundaryWidth)
       : KernelBase(p.size()),
         p(p),
         mg(mg),
@@ -40,7 +87,8 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
         cpy(cpy),
         cpz(cpz),
         ptype(ptype),
-        exclude(exclude)
+        exclude(exclude),
+        boundaryWidth(boundaryWidth)
   {
     runMessage();
     run();
@@ -54,73 +102,91 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
                  const ParticleDataImpl<Vec3> &cpy,
                  const ParticleDataImpl<Vec3> &cpz,
                  const ParticleDataImpl<int> *ptype,
-                 const int exclude)
+                 const int exclude,
+                 const int boundaryWidth)
   {
     if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
       return;
-    const IndexInt dX[2] = {0, vg.getStrideX()};
-    const IndexInt dY[2] = {0, vg.getStrideY()};
-    const IndexInt dZ[2] = {0, vg.getStrideZ()};
-
-    const Vec3 &pos = p[idx].pos, &vel = vp[idx];
-    const IndexInt fi = static_cast<IndexInt>(pos.x), fj = static_cast<IndexInt>(pos.y),
-                   fk = static_cast<IndexInt>(pos.z);
-    const IndexInt ci = static_cast<IndexInt>(pos.x - 0.5),
-                   cj = static_cast<IndexInt>(pos.y - 0.5),
-                   ck = static_cast<IndexInt>(pos.z - 0.5);
-    const Real wfi = clamp(pos.x - fi, Real(0), Real(1)),
-               wfj = clamp(pos.y - fj, Real(0), Real(1)),
-               wfk = clamp(pos.z - fk, Real(0), Real(1));
-    const Real wci = clamp(Real(pos.x - ci - 0.5), Real(0), Real(1)),
-               wcj = clamp(Real(pos.y - cj - 0.5), Real(0), Real(1)),
-               wck = clamp(Real(pos.z - ck - 0.5), Real(0), Real(1));
-    // TODO: check index for safety
+    if (!vg.isInBounds(p.getPos(idx), boundaryWidth)) {
+      debMsg("Skipping particle at index " << idx
+                                           << ". Is out of bounds and cannot be applied to grid.",
+             1);
+      return;
+    }
+
+    const Vec3 &pos = p.getPos(idx), &vel = vp[idx];
+    const Vec3i f = toVec3i(pos);
+    const Vec3i c = toVec3i(pos - 0.5);
+    const Vec3 wf = clamp(pos - toVec3(f), Vec3(0.), Vec3(1.));
+    const Vec3 wc = clamp(pos - toVec3(c) - 0.5, Vec3(0.), Vec3(1.));
+
     {  // u-face
-      const IndexInt gidx = fi * dX[1] + cj * dY[1] + ck * dZ[1];
-      const Vec3 gpos(fi, cj + 0.5, ck + 0.5);
-      const Real wi[2] = {Real(1) - wfi, wfi};
-      const Real wj[2] = {Real(1) - wcj, wcj};
-      const Real wk[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const Real w = wi[i] * wj[j] * wk[k];
-            mg[gidx + dX[i] + dY[j] + dZ[k]].x += w;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * vel.x;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * dot(cpx[idx], gpos + Vec3(i, j, k) - pos);
-          }
+      const IndexInt gidx = indexUFace(pos, vg);
+      if (gidx < 0)
+        return;  // debug will fail before
+
+      const Vec3 gpos(f.x, c.y + 0.5, c.z + 0.5);
+      const Real wi[2] = {Real(1) - wf.x, wf.x};
+      const Real wj[2] = {Real(1) - wc.y, wc.y};
+      const Real wk[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const Real w = wi[i] * wj[j] * wk[k];
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;  // debug will fail before
+
+        mg[vidx].x += w;
+        vg[vidx].x += w * vel.x;
+        vg[vidx].x += w * dot(cpx[idx], gpos + Vec3(i, j, k) - pos);
+      }
     }
     {  // v-face
-      const IndexInt gidx = ci * dX[1] + fj * dY[1] + ck * dZ[1];
-      const Vec3 gpos(ci + 0.5, fj, ck + 0.5);
-      const Real wi[2] = {Real(1) - wci, wci};
-      const Real wj[2] = {Real(1) - wfj, wfj};
-      const Real wk[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const Real w = wi[i] * wj[j] * wk[k];
-            mg[gidx + dX[i] + dY[j] + dZ[k]].y += w;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * vel.y;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * dot(cpy[idx], gpos + Vec3(i, j, k) - pos);
-          }
+      const IndexInt gidx = indexVFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Vec3 gpos(c.x + 0.5, f.y, c.z + 0.5);
+      const Real wi[2] = {Real(1) - wc.x, wc.x};
+      const Real wj[2] = {Real(1) - wf.y, wf.y};
+      const Real wk[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const Real w = wi[i] * wj[j] * wk[k];
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        mg[vidx].y += w;
+        vg[vidx].y += w * vel.y;
+        vg[vidx].y += w * dot(cpy[idx], gpos + Vec3(i, j, k) - pos);
+      }
     }
     if (!vg.is3D())
       return;
     {  // w-face
-      const IndexInt gidx = ci * dX[1] + cj * dY[1] + fk * dZ[1];
-      const Vec3 gpos(ci + 0.5, cj + 0.5, fk);
-      const Real wi[2] = {Real(1) - wci, wci};
-      const Real wj[2] = {Real(1) - wcj, wcj};
-      const Real wk[2] = {Real(1) - wfk, wfk};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const Real w = wi[i] * wj[j] * wk[k];
-            mg[gidx + dX[i] + dY[j] + dZ[k]].z += w;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * vel.z;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * dot(cpz[idx], gpos + Vec3(i, j, k) - pos);
-          }
+      const IndexInt gidx = indexWFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Vec3 gpos(c.x + 0.5, c.y + 0.5, f.z);
+      const Real wi[2] = {Real(1) - wc.x, wc.x};
+      const Real wj[2] = {Real(1) - wc.y, wc.y};
+      const Real wk[2] = {Real(1) - wf.z, wf.z};
+
+      FOR_INT_IJK(2)
+      {
+        const Real w = wi[i] * wj[j] * wk[k];
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        mg[vidx].z += w;
+        vg[vidx].z += w * vel.z;
+        vg[vidx].z += w * dot(cpz[idx], gpos + Vec3(i, j, k) - pos);
+      }
     }
   }
   inline const BasicParticleSystem &getArg0()
@@ -168,6 +234,11 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
     return exclude;
   }
   typedef int type8;
+  inline const int &getArg9()
+  {
+    return boundaryWidth;
+  }
+  typedef int type9;
   void runMessage()
   {
     debMsg("Executing kernel knApicMapLinearVec3ToMACGrid ", 3);
@@ -179,7 +250,7 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
   {
     const IndexInt _sz = size;
     for (IndexInt i = 0; i < _sz; i++)
-      op(i, p, mg, vg, vp, cpx, cpy, cpz, ptype, exclude);
+      op(i, p, mg, vg, vp, cpx, cpy, cpz, ptype, exclude, boundaryWidth);
   }
   const BasicParticleSystem &p;
   MACGrid &mg;
@@ -190,6 +261,7 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
   const ParticleDataImpl<Vec3> &cpz;
   const ParticleDataImpl<int> *ptype;
   const int exclude;
+  const int boundaryWidth;
 };
 
 void apicMapPartsToMAC(const FlagGrid &flags,
@@ -199,25 +271,24 @@ void apicMapPartsToMAC(const FlagGrid &flags,
                        const ParticleDataImpl<Vec3> &cpx,
                        const ParticleDataImpl<Vec3> &cpy,
                        const ParticleDataImpl<Vec3> &cpz,
-                       MACGrid *mass = NULL,
-                       const ParticleDataImpl<int> *ptype = NULL,
-                       const int exclude = 0)
+                       MACGrid *mass = nullptr,
+                       const ParticleDataImpl<int> *ptype = nullptr,
+                       const int exclude = 0,
+                       const int boundaryWidth = 0)
 {
-  // affine map
-  // let's assume that the particle mass is constant, 1.0
-  const bool freeMass = !mass;
-  if (!mass)
-    mass = new MACGrid(flags.getParent());
-  else
-    mass->clear();
+  // affine map: let's assume that the particle mass is constant, 1.0
+  MACGrid tmpmass(vel.getParent());
 
+  tmpmass.clear();
   vel.clear();
-  knApicMapLinearVec3ToMACGrid(parts, *mass, vel, partVel, cpx, cpy, cpz, ptype, exclude);
-  mass->stomp(VECTOR_EPSILON);
-  vel.safeDivide(*mass);
 
-  if (freeMass)
-    delete mass;
+  knApicMapLinearVec3ToMACGrid(
+      parts, tmpmass, vel, partVel, cpx, cpy, cpz, ptype, exclude, boundaryWidth);
+  tmpmass.stomp(VECTOR_EPSILON);
+  vel.safeDivide(tmpmass);
+
+  if (mass)
+    (*mass).swap(tmpmass);
 }
 static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
 {
@@ -226,7 +297,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "apicMapPartsToMAC", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -237,12 +308,14 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const ParticleDataImpl<Vec3> &cpx = *_args.getPtr<ParticleDataImpl<Vec3>>("cpx", 4, &_lock);
       const ParticleDataImpl<Vec3> &cpy = *_args.getPtr<ParticleDataImpl<Vec3>>("cpy", 5, &_lock);
       const ParticleDataImpl<Vec3> &cpz = *_args.getPtr<ParticleDataImpl<Vec3>>("cpz", 6, &_lock);
-      MACGrid *mass = _args.getPtrOpt<MACGrid>("mass", 7, NULL, &_lock);
+      MACGrid *mass = _args.getPtrOpt<MACGrid>("mass", 7, nullptr, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 8, NULL, &_lock);
+          "ptype", 8, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 9, 0, &_lock);
+      const int boundaryWidth = _args.getOpt<int>("boundaryWidth", 10, 0, &_lock);
       _retval = getPyNone();
-      apicMapPartsToMAC(flags, vel, parts, partVel, cpx, cpy, cpz, mass, ptype, exclude);
+      apicMapPartsToMAC(
+          flags, vel, parts, partVel, cpx, cpy, cpz, mass, ptype, exclude, boundaryWidth);
       _args.check();
     }
     pbFinalizePlugin(parent, "apicMapPartsToMAC", !noTiming);
@@ -270,7 +343,8 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
                                const MACGrid &vg,
                                const FlagGrid &flags,
                                const ParticleDataImpl<int> *ptype,
-                               const int exclude)
+                               const int exclude,
+                               const int boundaryWidth)
       : KernelBase(vp.size()),
         vp(vp),
         cpx(cpx),
@@ -280,7 +354,8 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
         vg(vg),
         flags(flags),
         ptype(ptype),
-        exclude(exclude)
+        exclude(exclude),
+        boundaryWidth(boundaryWidth)
   {
     runMessage();
     run();
@@ -294,78 +369,94 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
                  const MACGrid &vg,
                  const FlagGrid &flags,
                  const ParticleDataImpl<int> *ptype,
-                 const int exclude) const
+                 const int exclude,
+                 const int boundaryWidth) const
   {
     if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
       return;
+    if (!vg.isInBounds(p.getPos(idx), boundaryWidth)) {
+      debMsg("Skipping particle at index " << idx
+                                           << ". Is out of bounds and cannot get value from grid.",
+             1);
+      return;
+    }
 
     vp[idx] = cpx[idx] = cpy[idx] = cpz[idx] = Vec3(Real(0));
-    const IndexInt dX[2] = {0, vg.getStrideX()}, dY[2] = {0, vg.getStrideY()},
-                   dZ[2] = {0, vg.getStrideZ()};
     const Real gw[2] = {-Real(1), Real(1)};
 
-    const Vec3 &pos = p[idx].pos;
-    const IndexInt fi = static_cast<IndexInt>(pos.x), fj = static_cast<IndexInt>(pos.y),
-                   fk = static_cast<IndexInt>(pos.z);
-    const IndexInt ci = static_cast<IndexInt>(pos.x - 0.5),
-                   cj = static_cast<IndexInt>(pos.y - 0.5),
-                   ck = static_cast<IndexInt>(pos.z - 0.5);
-    const Real wfi = clamp(pos.x - fi, Real(0), Real(1)),
-               wfj = clamp(pos.y - fj, Real(0), Real(1)),
-               wfk = clamp(pos.z - fk, Real(0), Real(1));
-    const Real wci = clamp(Real(pos.x - ci - 0.5), Real(0), Real(1)),
-               wcj = clamp(Real(pos.y - cj - 0.5), Real(0), Real(1)),
-               wck = clamp(Real(pos.z - ck - 0.5), Real(0), Real(1));
-    // TODO: check index for safety
-    {  // u
-      const IndexInt gidx = fi * dX[1] + cj * dY[1] + ck * dZ[1];
-      const Real wx[2] = {Real(1) - wfi, wfi};
-      const Real wy[2] = {Real(1) - wcj, wcj};
-      const Real wz[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
-            Real vgx = vg[vidx].x;
-            vp[idx].x += wx[i] * wy[j] * wz[k] * vgx;
-            cpx[idx].x += gw[i] * wy[j] * wz[k] * vgx;
-            cpx[idx].y += wx[i] * gw[j] * wz[k] * vgx;
-            cpx[idx].z += wx[i] * wy[j] * gw[k] * vgx;
-          }
+    const Vec3 &pos = p.getPos(idx);
+    const Vec3i f = toVec3i(pos);
+    const Vec3i c = toVec3i(pos - 0.5);
+    const Vec3 wf = clamp(pos - toVec3(f), Vec3(0.), Vec3(1.));
+    const Vec3 wc = clamp(pos - toVec3(c) - 0.5, Vec3(0.), Vec3(1.));
+
+    {  // u-face
+      const IndexInt gidx = indexUFace(pos, vg);
+      if (gidx < 0)
+        return;  // debug will fail before
+
+      const Real wx[2] = {Real(1) - wf.x, wf.x};
+      const Real wy[2] = {Real(1) - wc.y, wc.y};
+      const Real wz[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;  // debug will fail before
+
+        const Real vgx = vg[vidx].x;
+        vp[idx].x += wx[i] * wy[j] * wz[k] * vgx;
+        cpx[idx].x += gw[i] * wy[j] * wz[k] * vgx;
+        cpx[idx].y += wx[i] * gw[j] * wz[k] * vgx;
+        cpx[idx].z += wx[i] * wy[j] * gw[k] * vgx;
+      }
     }
-    {  // v
-      const IndexInt gidx = ci * dX[1] + fj * dY[1] + ck * dZ[1];
-      const Real wx[2] = {Real(1) - wci, wci};
-      const Real wy[2] = {Real(1) - wfj, wfj};
-      const Real wz[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
-            Real vgy = vg[vidx].y;
-            vp[idx].y += wx[i] * wy[j] * wz[k] * vgy;
-            cpy[idx].x += gw[i] * wy[j] * wz[k] * vgy;
-            cpy[idx].y += wx[i] * gw[j] * wz[k] * vgy;
-            cpy[idx].z += wx[i] * wy[j] * gw[k] * vgy;
-          }
+    {  // v-face
+      const IndexInt gidx = indexVFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Real wx[2] = {Real(1) - wc.x, wc.x};
+      const Real wy[2] = {Real(1) - wf.y, wf.y};
+      const Real wz[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        const Real vgy = vg[vidx].y;
+        vp[idx].y += wx[i] * wy[j] * wz[k] * vgy;
+        cpy[idx].x += gw[i] * wy[j] * wz[k] * vgy;
+        cpy[idx].y += wx[i] * gw[j] * wz[k] * vgy;
+        cpy[idx].z += wx[i] * wy[j] * gw[k] * vgy;
+      }
     }
     if (!vg.is3D())
       return;
-    {  // w
-      const IndexInt gidx = ci * dX[1] + cj * dY[1] + fk * dZ[1];
-      const Real wx[2] = {Real(1) - wci, wci};
-      const Real wy[2] = {Real(1) - wcj, wcj};
-      const Real wz[2] = {Real(1) - wfk, wfk};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
-            Real vgz = vg[vidx].z;
-            vp[idx].z += wx[i] * wy[j] * wz[k] * vgz;
-            cpz[idx].x += gw[i] * wy[j] * wz[k] * vgz;
-            cpz[idx].y += wx[i] * gw[j] * wz[k] * vgz;
-            cpz[idx].z += wx[i] * wy[j] * gw[k] * vgz;
-          }
+    {  // w-face
+      const IndexInt gidx = indexWFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Real wx[2] = {Real(1) - wc.x, wc.x};
+      const Real wy[2] = {Real(1) - wc.y, wc.y};
+      const Real wz[2] = {Real(1) - wf.z, wf.z};
+
+      FOR_INT_IJK(2)
+      {
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        const Real vgz = vg[vidx].z;
+        vp[idx].z += wx[i] * wy[j] * wz[k] * vgz;
+        cpz[idx].x += gw[i] * wy[j] * wz[k] * vgz;
+        cpz[idx].y += wx[i] * gw[j] * wz[k] * vgz;
+        cpz[idx].z += wx[i] * wy[j] * gw[k] * vgz;
+      }
     }
   }
   inline ParticleDataImpl<Vec3> &getArg0()
@@ -413,6 +504,11 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
     return exclude;
   }
   typedef int type8;
+  inline const int &getArg9()
+  {
+    return boundaryWidth;
+  }
+  typedef int type9;
   void runMessage()
   {
     debMsg("Executing kernel knApicMapLinearMACGridToVec3 ", 3);
@@ -423,7 +519,7 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
   void operator()(const tbb::blocked_range<IndexInt> &__r) const
   {
     for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
-      op(idx, vp, cpx, cpy, cpz, p, vg, flags, ptype, exclude);
+      op(idx, vp, cpx, cpy, cpz, p, vg, flags, ptype, exclude, boundaryWidth);
   }
   void run()
   {
@@ -438,6 +534,7 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
   const FlagGrid &flags;
   const ParticleDataImpl<int> *ptype;
   const int exclude;
+  const int boundaryWidth;
 };
 
 void apicMapMACGridToParts(ParticleDataImpl<Vec3> &partVel,
@@ -447,10 +544,12 @@ void apicMapMACGridToParts(ParticleDataImpl<Vec3> &partVel,
                            const BasicParticleSystem &parts,
                            const MACGrid &vel,
                            const FlagGrid &flags,
-                           const ParticleDataImpl<int> *ptype = NULL,
-                           const int exclude = 0)
+                           const ParticleDataImpl<int> *ptype = nullptr,
+                           const int exclude = 0,
+                           const int boundaryWidth = 0)
 {
-  knApicMapLinearMACGridToVec3(partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude);
+  knApicMapLinearMACGridToVec3(
+      partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude, boundaryWidth);
 }
 static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
 {
@@ -459,7 +558,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "apicMapMACGridToParts", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
@@ -471,10 +570,12 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 5, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 6, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 7, NULL, &_lock);
+          "ptype", 7, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 8, 0, &_lock);
+      const int boundaryWidth = _args.getOpt<int>("boundaryWidth", 9, 0, &_lock);
       _retval = getPyNone();
-      apicMapMACGridToParts(partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude);
+      apicMapMACGridToParts(
+          partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude, boundaryWidth);
       _args.check();
     }
     pbFinalizePlugin(parent, "apicMapMACGridToParts", !noTiming);
diff --git a/extern/mantaflow/preprocessed/plugin/extforces.cpp b/extern/mantaflow/preprocessed/plugin/extforces.cpp
index 798bb3daeee..558008afe3b 100644
--- a/extern/mantaflow/preprocessed/plugin/extforces.cpp
+++ b/extern/mantaflow/preprocessed/plugin/extforces.cpp
@@ -248,7 +248,7 @@ struct KnApplyForce : public KernelBase {
 void addGravity(const FlagGrid &flags,
                 MACGrid &vel,
                 Vec3 gravity,
-                const Grid<Real> *exclude = NULL,
+                const Grid<Real> *exclude = nullptr,
                 bool scale = true)
 {
   float gridScale = (scale) ? flags.getDx() : 1;
@@ -262,13 +262,13 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "addGravity", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
       Vec3 gravity = _args.get<Vec3>("gravity", 2, &_lock);
-      const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 3, NULL, &_lock);
+      const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 3, nullptr, &_lock);
       bool scale = _args.getOpt<bool>("scale", 4, true, &_lock);
       _retval = getPyNone();
       addGravity(flags, vel, gravity, exclude, scale);
@@ -294,7 +294,7 @@ void PbRegister_addGravity()
 void addGravityNoScale(const FlagGrid &flags,
                        MACGrid &vel,
                        const Vec3 &gravity,
-                       const Grid<Real> *exclude = NULL)
+                       const Grid<Real> *exclude = nullptr)
 {
   addGravity(flags, vel, gravity, exclude, false);
 }
@@ -305,13 +305,13 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "addGravityNoScale", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
       const Vec3 &gravity = _args.get<Vec3>("gravity", 2, &_lock);
-      const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 3, NULL, &_lock);
+      const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 3, nullptr, &_lock);
       _retval = getPyNone();
       addGravityNoScale(flags, vel, gravity, exclude);
       _args.check();
@@ -434,7 +434,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "addBuoyancy", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -537,7 +537,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setOpenBound", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -567,11 +567,11 @@ void PbRegister_setOpenBound()
 //! delete fluid and ensure empty flag in outflow cells, delete particles and density and set phi
 //! to 0.5
 void resetOutflow(FlagGrid &flags,
-                  Grid<Real> *phi = 0,
-                  BasicParticleSystem *parts = 0,
-                  Grid<Real> *real = 0,
-                  Grid<int> *index = 0,
-                  ParticleIndexSystem *indexSys = 0)
+                  Grid<Real> *phi = nullptr,
+                  BasicParticleSystem *parts = nullptr,
+                  Grid<Real> *real = nullptr,
+                  Grid<int> *index = nullptr,
+                  ParticleIndexSystem *indexSys = nullptr)
 {
   // check if phi and parts -> pindex and gpi already created -> access particles from cell index,
   // avoid extra looping over particles
@@ -623,16 +623,17 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "resetOutflow", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
-      Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 1, 0, &_lock);
-      BasicParticleSystem *parts = _args.getPtrOpt<BasicParticleSystem>("parts", 2, 0, &_lock);
-      Grid<Real> *real = _args.getPtrOpt<Grid<Real>>("real", 3, 0, &_lock);
-      Grid<int> *index = _args.getPtrOpt<Grid<int>>("index", 4, 0, &_lock);
+      Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 1, nullptr, &_lock);
+      BasicParticleSystem *parts = _args.getPtrOpt<BasicParticleSystem>(
+          "parts", 2, nullptr, &_lock);
+      Grid<Real> *real = _args.getPtrOpt<Grid<Real>>("real", 3, nullptr, &_lock);
+      Grid<int> *index = _args.getPtrOpt<Grid<int>>("index", 4, nullptr, &_lock);
       ParticleIndexSystem *indexSys = _args.getPtrOpt<ParticleIndexSystem>(
-          "indexSys", 5, 0, &_lock);
+          "indexSys", 5, nullptr, &_lock);
       _retval = getPyNone();
       resetOutflow(flags, phi, parts, real, index, indexSys);
       _args.check();
@@ -747,7 +748,7 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setInflowBcs", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
@@ -1064,13 +1065,13 @@ struct KnSetWallBcsFrac : public KernelBase {
 };
 
 //! set zero normal velocity boundary condition on walls
-// (optionally with second order accuracy using the obstacle SDF , fractions grid currentlyl not
+// (optionally with second order accuracy using the obstacle SDF , fractions grid currently not
 // needed)
 void setWallBcs(const FlagGrid &flags,
                 MACGrid &vel,
-                const MACGrid *obvel = 0,
-                const MACGrid *fractions = 0,
-                const Grid<Real> *phiObs = 0,
+                const MACGrid *obvel = nullptr,
+                const MACGrid *fractions = nullptr,
+                const Grid<Real> *phiObs = nullptr,
                 int boundaryWidth = 0)
 {
   if (!phiObs || !fractions) {
@@ -1089,14 +1090,14 @@ static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setWallBcs", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
-      const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 2, 0, &_lock);
-      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 3, 0, &_lock);
-      const Grid<Real> *phiObs = _args.getPtrOpt<Grid<Real>>("phiObs", 4, 0, &_lock);
+      const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 2, nullptr, &_lock);
+      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 3, nullptr, &_lock);
+      const Grid<Real> *phiObs = _args.getPtrOpt<Grid<Real>>("phiObs", 4, nullptr, &_lock);
       int boundaryWidth = _args.getOpt<int>("boundaryWidth", 5, 0, &_lock);
       _retval = getPyNone();
       setWallBcs(flags, vel, obvel, fractions, phiObs, boundaryWidth);
@@ -1219,7 +1220,7 @@ static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setInitialVelocity", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -1342,7 +1343,7 @@ struct KnConfForce : public KernelBase {
 void vorticityConfinement(MACGrid &vel,
                           const FlagGrid &flags,
                           Real strength = 0,
-                          const Grid<Real> *strengthCell = NULL)
+                          const Grid<Real> *strengthCell = nullptr)
 {
   Grid<Vec3> velCenter(flags.getParent()), curl(flags.getParent()), force(flags.getParent());
   Grid<Real> norm(flags.getParent());
@@ -1351,7 +1352,7 @@ void vorticityConfinement(MACGrid &vel,
   CurlOp(velCenter, curl);
   GridNorm(norm, curl);
   KnConfForce(force, norm, curl, strength, strengthCell);
-  KnApplyForceField(flags, vel, force, NULL, true, false);
+  KnApplyForceField(flags, vel, force, nullptr, true, false);
 }
 static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
 {
@@ -1360,14 +1361,14 @@ static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "vorticityConfinement", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
       Real strength = _args.getOpt<Real>("strength", 2, 0, &_lock);
       const Grid<Real> *strengthCell = _args.getPtrOpt<Grid<Real>>(
-          "strengthCell", 3, NULL, &_lock);
+          "strengthCell", 3, nullptr, &_lock);
       _retval = getPyNone();
       vorticityConfinement(vel, flags, strength, strengthCell);
       _args.check();
@@ -1391,7 +1392,7 @@ void PbRegister_vorticityConfinement()
 void addForceField(const FlagGrid &flags,
                    MACGrid &vel,
                    const Grid<Vec3> &force,
-                   const Grid<Real> *region = NULL,
+                   const Grid<Real> *region = nullptr,
                    bool isMAC = false)
 {
   KnApplyForceField(flags, vel, force, region, true, isMAC);
@@ -1403,13 +1404,13 @@ static PyObject *_W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "addForceField", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
       const Grid<Vec3> &force = *_args.getPtr<Grid<Vec3>>("force", 2, &_lock);
-      const Grid<Real> *region = _args.getPtrOpt<Grid<Real>>("region", 3, NULL, &_lock);
+      const Grid<Real> *region = _args.getPtrOpt<Grid<Real>>("region", 3, nullptr, &_lock);
       bool isMAC = _args.getOpt<bool>("isMAC", 4, false, &_lock);
       _retval = getPyNone();
       addForceField(flags, vel, force, region, isMAC);
@@ -1434,7 +1435,7 @@ void PbRegister_addForceField()
 void setForceField(const FlagGrid &flags,
                    MACGrid &vel,
                    const Grid<Vec3> &force,
-                   const Grid<Real> *region = NULL,
+                   const Grid<Real> *region = nullptr,
                    bool isMAC = false)
 {
   KnApplyForceField(flags, vel, force, region, false, isMAC);
@@ -1446,13 +1447,13 @@ static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setForceField", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
       const Grid<Vec3> &force = *_args.getPtr<Grid<Vec3>>("force", 2, &_lock);
-      const Grid<Real> *region = _args.getPtrOpt<Grid<Real>>("region", 3, NULL, &_lock);
+      const Grid<Real> *region = _args.getPtrOpt<Grid<Real>>("region", 3, nullptr, &_lock);
       bool isMAC = _args.getOpt<bool>("isMAC", 4, false, &_lock);
       _retval = getPyNone();
       setForceField(flags, vel, force, region, isMAC);
@@ -1645,10 +1646,10 @@ struct KnDissolveSmoke : public KernelBase {
 
 void dissolveSmoke(const FlagGrid &flags,
                    Grid<Real> &density,
-                   Grid<Real> *heat = NULL,
-                   Grid<Real> *red = NULL,
-                   Grid<Real> *green = NULL,
-                   Grid<Real> *blue = NULL,
+                   Grid<Real> *heat = nullptr,
+                   Grid<Real> *red = nullptr,
+                   Grid<Real> *green = nullptr,
+                   Grid<Real> *blue = nullptr,
                    int speed = 5,
                    bool logFalloff = true)
 {
@@ -1663,15 +1664,15 @@ static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "dissolveSmoke", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 1, &_lock);
-      Grid<Real> *heat = _args.getPtrOpt<Grid<Real>>("heat", 2, NULL, &_lock);
-      Grid<Real> *red = _args.getPtrOpt<Grid<Real>>("red", 3, NULL, &_lock);
-      Grid<Real> *green = _args.getPtrOpt<Grid<Real>>("green", 4, NULL, &_lock);
-      Grid<Real> *blue = _args.getPtrOpt<Grid<Real>>("blue", 5, NULL, &_lock);
+      Grid<Real> *heat = _args.getPtrOpt<Grid<Real>>("heat", 2, nullptr, &_lock);
+      Grid<Real> *red = _args.getPtrOpt<Grid<Real>>("red", 3, nullptr, &_lock);
+      Grid<Real> *green = _args.getPtrOpt<Grid<Real>>("green", 4, nullptr, &_lock);
+      Grid<Real> *blue = _args.getPtrOpt<Grid<Real>>("blue", 5, nullptr, &_lock);
       int speed = _args.getOpt<int>("speed", 6, 5, &_lock);
       bool logFalloff = _args.getOpt<bool>("logFalloff", 7, true, &_lock);
       _retval = getPyNone();
diff --git a/extern/mantaflow/preprocessed/plugin/fire.cpp b/extern/mantaflow/preprocessed/plugin/fire.cpp
index 9047d4bf8a1..f907bdf0504 100644
--- a/extern/mantaflow/preprocessed/plugin/fire.cpp
+++ b/extern/mantaflow/preprocessed/plugin/fire.cpp
@@ -258,10 +258,10 @@ struct KnProcessBurn : public KernelBase {
 void processBurn(Grid<Real> &fuel,
                  Grid<Real> &density,
                  Grid<Real> &react,
-                 Grid<Real> *red = NULL,
-                 Grid<Real> *green = NULL,
-                 Grid<Real> *blue = NULL,
-                 Grid<Real> *heat = NULL,
+                 Grid<Real> *red = nullptr,
+                 Grid<Real> *green = nullptr,
+                 Grid<Real> *blue = nullptr,
+                 Grid<Real> *heat = nullptr,
                  Real burningRate = 0.75f,
                  Real flameSmoke = 1.0f,
                  Real ignitionTemp = 1.25f,
@@ -290,16 +290,16 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "processBurn", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &fuel = *_args.getPtr<Grid<Real>>("fuel", 0, &_lock);
       Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 1, &_lock);
       Grid<Real> &react = *_args.getPtr<Grid<Real>>("react", 2, &_lock);
-      Grid<Real> *red = _args.getPtrOpt<Grid<Real>>("red", 3, NULL, &_lock);
-      Grid<Real> *green = _args.getPtrOpt<Grid<Real>>("green", 4, NULL, &_lock);
-      Grid<Real> *blue = _args.getPtrOpt<Grid<Real>>("blue", 5, NULL, &_lock);
-      Grid<Real> *heat = _args.getPtrOpt<Grid<Real>>("heat", 6, NULL, &_lock);
+      Grid<Real> *red = _args.getPtrOpt<Grid<Real>>("red", 3, nullptr, &_lock);
+      Grid<Real> *green = _args.getPtrOpt<Grid<Real>>("green", 4, nullptr, &_lock);
+      Grid<Real> *blue = _args.getPtrOpt<Grid<Real>>("blue", 5, nullptr, &_lock);
+      Grid<Real> *heat = _args.getPtrOpt<Grid<Real>>("heat", 6, nullptr, &_lock);
       Real burningRate = _args.getOpt<Real>("burningRate", 7, 0.75f, &_lock);
       Real flameSmoke = _args.getOpt<Real>("flameSmoke", 8, 1.0f, &_lock);
       Real ignitionTemp = _args.getOpt<Real>("ignitionTemp", 9, 1.25f, &_lock);
@@ -407,7 +407,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "updateFlame", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<Real> &react = *_args.getPtr<Grid<Real>>("react", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/flip.cpp b/extern/mantaflow/preprocessed/plugin/flip.cpp
index 8ac167c1166..1acdac1c094 100644
--- a/extern/mantaflow/preprocessed/plugin/flip.cpp
+++ b/extern/mantaflow/preprocessed/plugin/flip.cpp
@@ -69,7 +69,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "sampleFlagsWithParticles", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -156,7 +156,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "sampleLevelsetWithParticles", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 0, &_lock);
@@ -200,7 +200,7 @@ void sampleShapeWithParticles(const Shape &shape,
                               const Real randomness,
                               const bool reset = false,
                               const bool refillEmpty = false,
-                              const LevelsetGrid *exclude = NULL)
+                              const LevelsetGrid *exclude = nullptr)
 {
   const bool is3D = flags.is3D();
   const Real jlen = randomness / discretization;
@@ -243,7 +243,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "sampleShapeWithParticles", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Shape &shape = *_args.getPtr<Shape>("shape", 0, &_lock);
@@ -253,7 +253,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const Real randomness = _args.get<Real>("randomness", 4, &_lock);
       const bool reset = _args.getOpt<bool>("reset", 5, false, &_lock);
       const bool refillEmpty = _args.getOpt<bool>("refillEmpty", 6, false, &_lock);
-      const LevelsetGrid *exclude = _args.getPtrOpt<LevelsetGrid>("exclude", 7, NULL, &_lock);
+      const LevelsetGrid *exclude = _args.getPtrOpt<LevelsetGrid>("exclude", 7, nullptr, &_lock);
       _retval = getPyNone();
       sampleShapeWithParticles(
           shape, flags, parts, discretization, randomness, reset, refillEmpty, exclude);
@@ -418,8 +418,8 @@ struct knSetNbObstacle : public KernelBase {
 };
 void markFluidCells(const BasicParticleSystem &parts,
                     FlagGrid &flags,
-                    const Grid<Real> *phiObs = NULL,
-                    const ParticleDataImpl<int> *ptype = NULL,
+                    const Grid<Real> *phiObs = nullptr,
+                    const ParticleDataImpl<int> *ptype = nullptr,
                     const int exclude = 0)
 {
   // remove all fluid cells
@@ -448,14 +448,14 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "markFluidCells", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
       FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
-      const Grid<Real> *phiObs = _args.getPtrOpt<Grid<Real>>("phiObs", 2, NULL, &_lock);
+      const Grid<Real> *phiObs = _args.getPtrOpt<Grid<Real>>("phiObs", 2, nullptr, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 3, NULL, &_lock);
+          "ptype", 3, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 4, 0, &_lock);
       _retval = getPyNone();
       markFluidCells(parts, flags, phiObs, ptype, exclude);
@@ -492,7 +492,7 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "testInitGridWithPos", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &grid = *_args.getPtr<Grid<Real>>("grid", 0, &_lock);
@@ -534,7 +534,7 @@ void adjustNumber(BasicParticleSystem &parts,
                   const LevelsetGrid &phi,
                   Real radiusFactor = 1.,
                   Real narrowBand = -1.,
-                  const Grid<Real> *exclude = NULL)
+                  const Grid<Real> *exclude = nullptr)
 {
   // which levelset to use as threshold
   const Real SURFACE_LS = -1.0 * calculateRadiusFactor(phi, radiusFactor);
@@ -610,7 +610,7 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "adjustNumber", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -621,7 +621,7 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 5, &_lock);
       Real radiusFactor = _args.getOpt<Real>("radiusFactor", 6, 1., &_lock);
       Real narrowBand = _args.getOpt<Real>("narrowBand", 7, -1., &_lock);
-      const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 8, NULL, &_lock);
+      const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 8, nullptr, &_lock);
       _retval = getPyNone();
       adjustNumber(
           parts, vel, flags, minParticles, maxParticles, phi, radiusFactor, narrowBand, exclude);
@@ -660,7 +660,7 @@ static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "debugIntToReal", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<int> &source = *_args.getPtr<Grid<int>>("source", 0, &_lock);
@@ -694,7 +694,7 @@ void gridParticleIndex(const BasicParticleSystem &parts,
                        ParticleIndexSystem &indexSys,
                        const FlagGrid &flags,
                        Grid<int> &index,
-                       Grid<int> *counter = NULL)
+                       Grid<int> *counter = nullptr)
 {
   bool delCounter = false;
   if (!counter) {
@@ -761,14 +761,14 @@ static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "gridParticleIndex", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
       ParticleIndexSystem &indexSys = *_args.getPtr<ParticleIndexSystem>("indexSys", 1, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
       Grid<int> &index = *_args.getPtr<Grid<int>>("index", 3, &_lock);
-      Grid<int> *counter = _args.getPtrOpt<Grid<int>>("counter", 4, NULL, &_lock);
+      Grid<int> *counter = _args.getPtrOpt<Grid<int>>("counter", 4, nullptr, &_lock);
       _retval = getPyNone();
       gridParticleIndex(parts, indexSys, flags, index, counter);
       _args.check();
@@ -932,7 +932,7 @@ void unionParticleLevelset(const BasicParticleSystem &parts,
                            const Grid<int> &index,
                            LevelsetGrid &phi,
                            const Real radiusFactor = 1.,
-                           const ParticleDataImpl<int> *ptype = NULL,
+                           const ParticleDataImpl<int> *ptype = nullptr,
                            const int exclude = 0)
 {
   // use half a cell diagonal as base radius
@@ -949,7 +949,7 @@ static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "unionParticleLevelset", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -960,7 +960,7 @@ static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 4, &_lock);
       const Real radiusFactor = _args.getOpt<Real>("radiusFactor", 5, 1., &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 6, NULL, &_lock);
+          "ptype", 6, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 7, 0, &_lock);
       _retval = getPyNone();
       unionParticleLevelset(parts, indexSys, flags, index, phi, radiusFactor, ptype, exclude);
@@ -992,8 +992,8 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
                                 const Real radius,
                                 const ParticleDataImpl<int> *ptype,
                                 const int exclude,
-                                Grid<Vec3> *save_pAcc = NULL,
-                                Grid<Real> *save_rAcc = NULL)
+                                Grid<Vec3> *save_pAcc = nullptr,
+                                Grid<Real> *save_rAcc = nullptr)
       : KernelBase(&index, 0),
         parts(parts),
         index(index),
@@ -1018,8 +1018,8 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
                  const Real radius,
                  const ParticleDataImpl<int> *ptype,
                  const int exclude,
-                 Grid<Vec3> *save_pAcc = NULL,
-                 Grid<Real> *save_rAcc = NULL) const
+                 Grid<Vec3> *save_pAcc = nullptr,
+                 Grid<Real> *save_rAcc = nullptr) const
   {
     const Vec3 gridPos = Vec3(i, j, k) + Vec3(0.5);  // shifted by half cell
     Real phiv = radius * 1.0;                        // outside
@@ -1312,7 +1312,7 @@ void averagedParticleLevelset(const BasicParticleSystem &parts,
                               const Real radiusFactor = 1.,
                               const int smoothen = 1,
                               const int smoothenNeg = 1,
-                              const ParticleDataImpl<int> *ptype = NULL,
+                              const ParticleDataImpl<int> *ptype = nullptr,
                               const int exclude = 0)
 {
   // use half a cell diagonal as base radius
@@ -1340,7 +1340,7 @@ static PyObject *_W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "averagedParticleLevelset", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -1353,7 +1353,7 @@ static PyObject *_W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const int smoothen = _args.getOpt<int>("smoothen", 6, 1, &_lock);
       const int smoothenNeg = _args.getOpt<int>("smoothenNeg", 7, 1, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 8, NULL, &_lock);
+          "ptype", 8, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 9, 0, &_lock);
       _retval = getPyNone();
       averagedParticleLevelset(
@@ -1523,7 +1523,7 @@ void improvedParticleLevelset(const BasicParticleSystem &parts,
                               const int smoothenNeg = 1,
                               const Real t_low = 0.4,
                               const Real t_high = 3.5,
-                              const ParticleDataImpl<int> *ptype = NULL,
+                              const ParticleDataImpl<int> *ptype = nullptr,
                               const int exclude = 0)
 {
   // create temporary grids to store values from levelset weight computation
@@ -1557,7 +1557,7 @@ static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "improvedParticleLevelset", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -1572,7 +1572,7 @@ static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const Real t_low = _args.getOpt<Real>("t_low", 8, 0.4, &_lock);
       const Real t_high = _args.getOpt<Real>("t_high", 9, 3.5, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 10, NULL, &_lock);
+          "ptype", 10, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 11, 0, &_lock);
       _retval = getPyNone();
       improvedParticleLevelset(parts,
@@ -1714,7 +1714,7 @@ void pushOutofObs(BasicParticleSystem &parts,
                   const Grid<Real> &phiObs,
                   const Real shift = 0,
                   const Real thresh = 0,
-                  const ParticleDataImpl<int> *ptype = NULL,
+                  const ParticleDataImpl<int> *ptype = nullptr,
                   const int exclude = 0)
 {
   knPushOutofObs(parts, flags, phiObs, shift, thresh, ptype, exclude);
@@ -1726,7 +1726,7 @@ static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "pushOutofObs", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -1735,7 +1735,7 @@ static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const Real shift = _args.getOpt<Real>("shift", 3, 0, &_lock);
       const Real thresh = _args.getOpt<Real>("thresh", 4, 0, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 5, NULL, &_lock);
+          "ptype", 5, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 6, 0, &_lock);
       _retval = getPyNone();
       pushOutofObs(parts, flags, phiObs, shift, thresh, ptype, exclude);
@@ -1917,8 +1917,8 @@ void mapPartsToMAC(const FlagGrid &flags,
                    MACGrid &velOld,
                    const BasicParticleSystem &parts,
                    const ParticleDataImpl<Vec3> &partVel,
-                   Grid<Vec3> *weight = NULL,
-                   const ParticleDataImpl<int> *ptype = NULL,
+                   Grid<Vec3> *weight = nullptr,
+                   const ParticleDataImpl<int> *ptype = nullptr,
                    const int exclude = 0)
 {
   // interpol -> grid. tmpgrid for particle contribution weights
@@ -1949,7 +1949,7 @@ static PyObject *_W_12(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "mapPartsToMAC", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -1958,9 +1958,9 @@ static PyObject *_W_12(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 3, &_lock);
       const ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
           "partVel", 4, &_lock);
-      Grid<Vec3> *weight = _args.getPtrOpt<Grid<Vec3>>("weight", 5, NULL, &_lock);
+      Grid<Vec3> *weight = _args.getPtrOpt<Grid<Vec3>>("weight", 5, nullptr, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 6, NULL, &_lock);
+          "ptype", 6, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 7, 0, &_lock);
       _retval = getPyNone();
       mapPartsToMAC(flags, vel, velOld, parts, partVel, weight, ptype, exclude);
@@ -2076,7 +2076,7 @@ static PyObject *_W_13(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "mapPartsToGrid", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2118,7 +2118,7 @@ static PyObject *_W_14(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "mapPartsToGridVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2215,7 +2215,7 @@ static PyObject *_W_15(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "mapGridToParts", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<Real> &source = *_args.getPtr<Grid<Real>>("source", 0, &_lock);
@@ -2254,7 +2254,7 @@ static PyObject *_W_16(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "mapGridToPartsVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<Vec3> &source = *_args.getPtr<Grid<Vec3>>("source", 0, &_lock);
@@ -2371,7 +2371,7 @@ void mapMACToParts(const FlagGrid &flags,
                    const MACGrid &vel,
                    const BasicParticleSystem &parts,
                    ParticleDataImpl<Vec3> &partVel,
-                   const ParticleDataImpl<int> *ptype = NULL,
+                   const ParticleDataImpl<int> *ptype = nullptr,
                    const int exclude = 0)
 {
   knMapLinearMACGridToVec3_PIC(parts, flags, vel, partVel, ptype, exclude);
@@ -2383,7 +2383,7 @@ static PyObject *_W_17(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "mapMACToParts", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2392,7 +2392,7 @@ static PyObject *_W_17(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
           "partVel", 3, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 4, NULL, &_lock);
+          "ptype", 4, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 5, 0, &_lock);
       _retval = getPyNone();
       mapMACToParts(flags, vel, parts, partVel, ptype, exclude);
@@ -2526,7 +2526,7 @@ void flipVelocityUpdate(const FlagGrid &flags,
                         const BasicParticleSystem &parts,
                         ParticleDataImpl<Vec3> &partVel,
                         const Real flipRatio,
-                        const ParticleDataImpl<int> *ptype = NULL,
+                        const ParticleDataImpl<int> *ptype = nullptr,
                         const int exclude = 0)
 {
   knMapLinearMACGridToVec3_FLIP(parts, flags, vel, velOld, partVel, flipRatio, ptype, exclude);
@@ -2538,7 +2538,7 @@ static PyObject *_W_18(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipVelocityUpdate", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2549,7 +2549,7 @@ static PyObject *_W_18(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
           "partVel", 4, &_lock);
       const Real flipRatio = _args.get<Real>("flipRatio", 5, &_lock);
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
-          "ptype", 6, NULL, &_lock);
+          "ptype", 6, nullptr, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 7, 0, &_lock);
       _retval = getPyNone();
       flipVelocityUpdate(flags, vel, velOld, parts, partVel, flipRatio, ptype, exclude);
@@ -2700,7 +2700,7 @@ struct knCombineVels : public KernelBase {
 void combineGridVel(MACGrid &vel,
                     const Grid<Vec3> &weight,
                     MACGrid &combineVel,
-                    const LevelsetGrid *phi = NULL,
+                    const LevelsetGrid *phi = nullptr,
                     Real narrowBand = 0.0,
                     Real thresh = 0.0)
 {
@@ -2713,13 +2713,13 @@ static PyObject *_W_19(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "combineGridVel", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
       const Grid<Vec3> &weight = *_args.getPtr<Grid<Vec3>>("weight", 1, &_lock);
       MACGrid &combineVel = *_args.getPtr<MACGrid>("combineVel", 2, &_lock);
-      const LevelsetGrid *phi = _args.getPtrOpt<LevelsetGrid>("phi", 3, NULL, &_lock);
+      const LevelsetGrid *phi = _args.getPtrOpt<LevelsetGrid>("phi", 3, nullptr, &_lock);
       Real narrowBand = _args.getOpt<Real>("narrowBand", 4, 0.0, &_lock);
       Real thresh = _args.getOpt<Real>("thresh", 5, 0.0, &_lock);
       _retval = getPyNone();
@@ -2754,7 +2754,7 @@ static PyObject *_W_20(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "getLaplacian", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &laplacian = *_args.getPtr<Grid<Real>>("laplacian", 0, &_lock);
@@ -2790,7 +2790,7 @@ static PyObject *_W_21(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "getCurvature", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &curv = *_args.getPtr<Grid<Real>>("curv", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/fluidguiding.cpp b/extern/mantaflow/preprocessed/plugin/fluidguiding.cpp
index 18a5a37771f..de881840a2e 100644
--- a/extern/mantaflow/preprocessed/plugin/fluidguiding.cpp
+++ b/extern/mantaflow/preprocessed/plugin/fluidguiding.cpp
@@ -404,7 +404,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "getSpiralVelocity", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -454,7 +454,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setGradientYWeight", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &W = *_args.getPtr<Grid<Real>>("W", 0, &_lock);
@@ -574,10 +574,10 @@ void solvePressure(MACGrid &vel,
                    Grid<Real> &pressure,
                    const FlagGrid &flags,
                    Real cgAccuracy = 1e-3,
-                   const Grid<Real> *phi = 0,
-                   const Grid<Real> *perCellCorr = 0,
-                   const MACGrid *fractions = 0,
-                   const MACGrid *obvel = 0,
+                   const Grid<Real> *phi = nullptr,
+                   const Grid<Real> *perCellCorr = nullptr,
+                   const MACGrid *fractions = nullptr,
+                   const MACGrid *obvel = nullptr,
                    Real gfClamp = 1e-04,
                    Real cgMaxIterFac = 1.5,
                    bool precondition = true,
@@ -585,9 +585,9 @@ void solvePressure(MACGrid &vel,
                    bool enforceCompatibility = false,
                    bool useL2Norm = false,
                    bool zeroPressureFixing = false,
-                   const Grid<Real> *curv = NULL,
+                   const Grid<Real> *curv = nullptr,
                    const Real surfTens = 0.0,
-                   Grid<Real> *retRhs = NULL);
+                   Grid<Real> *retRhs = nullptr);
 
 //! Main function for fluid guiding , includes "regular" pressure solve
 
@@ -603,16 +603,16 @@ void PD_fluid_guiding(MACGrid &vel,
                       Real epsRel = 1e-3,
                       Real epsAbs = 1e-3,
                       int maxIters = 200,
-                      Grid<Real> *phi = 0,
-                      Grid<Real> *perCellCorr = 0,
-                      MACGrid *fractions = 0,
-                      MACGrid *obvel = 0,
+                      Grid<Real> *phi = nullptr,
+                      Grid<Real> *perCellCorr = nullptr,
+                      MACGrid *fractions = nullptr,
+                      MACGrid *obvel = nullptr,
                       Real gfClamp = 1e-04,
                       Real cgMaxIterFac = 1.5,
                       Real cgAccuracy = 1e-3,
                       int preconditioner = 1,
                       bool zeroPressureFixing = false,
-                      const Grid<Real> *curv = NULL,
+                      const Grid<Real> *curv = nullptr,
                       const Real surfTens = 0.)
 {
   FluidSolver *parent = vel.getParent();
@@ -693,7 +693,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "PD_fluid_guiding", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
@@ -708,16 +708,16 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       Real epsRel = _args.getOpt<Real>("epsRel", 9, 1e-3, &_lock);
       Real epsAbs = _args.getOpt<Real>("epsAbs", 10, 1e-3, &_lock);
       int maxIters = _args.getOpt<int>("maxIters", 11, 200, &_lock);
-      Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 12, 0, &_lock);
-      Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 13, 0, &_lock);
-      MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 14, 0, &_lock);
-      MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 15, 0, &_lock);
+      Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 12, nullptr, &_lock);
+      Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 13, nullptr, &_lock);
+      MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 14, nullptr, &_lock);
+      MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 15, nullptr, &_lock);
       Real gfClamp = _args.getOpt<Real>("gfClamp", 16, 1e-04, &_lock);
       Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 17, 1.5, &_lock);
       Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 18, 1e-3, &_lock);
       int preconditioner = _args.getOpt<int>("preconditioner", 19, 1, &_lock);
       bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 20, false, &_lock);
-      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 21, NULL, &_lock);
+      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 21, nullptr, &_lock);
       const Real surfTens = _args.getOpt<Real>("surfTens", 22, 0., &_lock);
       _retval = getPyNone();
       PD_fluid_guiding(vel,
@@ -775,7 +775,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "releaseBlurPrecomp", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       _retval = getPyNone();
diff --git a/extern/mantaflow/preprocessed/plugin/initplugins.cpp b/extern/mantaflow/preprocessed/plugin/initplugins.cpp
index 6ccd3afc8d1..dce7b72de6c 100644
--- a/extern/mantaflow/preprocessed/plugin/initplugins.cpp
+++ b/extern/mantaflow/preprocessed/plugin/initplugins.cpp
@@ -154,7 +154,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "densityInflow", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -273,7 +273,7 @@ struct KnAddNoise : public KernelBase {
 void addNoise(const FlagGrid &flags,
               Grid<Real> &density,
               const WaveletNoiseField &noise,
-              const Grid<Real> *sdf = NULL,
+              const Grid<Real> *sdf = nullptr,
               Real scale = 1.0)
 {
   KnAddNoise(flags, density, noise, sdf, scale);
@@ -285,13 +285,13 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "addNoise", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 1, &_lock);
       const WaveletNoiseField &noise = *_args.getPtr<WaveletNoiseField>("noise", 2, &_lock);
-      const Grid<Real> *sdf = _args.getPtrOpt<Grid<Real>>("sdf", 3, NULL, &_lock);
+      const Grid<Real> *sdf = _args.getPtrOpt<Grid<Real>>("sdf", 3, nullptr, &_lock);
       Real scale = _args.getOpt<Real>("scale", 4, 1.0, &_lock);
       _retval = getPyNone();
       addNoise(flags, density, noise, sdf, scale);
@@ -448,7 +448,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setNoisePdata", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -489,7 +489,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setNoisePdataVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -530,7 +530,7 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setNoisePdataInt", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -592,7 +592,7 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "obstacleGradient", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -636,7 +636,7 @@ static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "obstacleLevelset", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -777,7 +777,7 @@ struct KnApplyEmission : public KernelBase {
 void applyEmission(FlagGrid &flags,
                    Grid<Real> &target,
                    Grid<Real> &source,
-                   Grid<Real> *emissionTexture = NULL,
+                   Grid<Real> *emissionTexture = nullptr,
                    bool isAbsolute = true,
                    int type = 0)
 {
@@ -790,14 +790,14 @@ static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "applyEmission", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       Grid<Real> &target = *_args.getPtr<Grid<Real>>("target", 1, &_lock);
       Grid<Real> &source = *_args.getPtr<Grid<Real>>("source", 2, &_lock);
       Grid<Real> *emissionTexture = _args.getPtrOpt<Grid<Real>>(
-          "emissionTexture", 3, NULL, &_lock);
+          "emissionTexture", 3, nullptr, &_lock);
       bool isAbsolute = _args.getOpt<bool>("isAbsolute", 4, true, &_lock);
       int type = _args.getOpt<int>("type", 5, 0, &_lock);
       _retval = getPyNone();
@@ -925,7 +925,7 @@ static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "densityInflowMeshNoise", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -974,7 +974,7 @@ static PyObject *_W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "densityInflowMesh", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -1161,12 +1161,12 @@ struct KnResetInObstacle : public KernelBase {
 void resetInObstacle(FlagGrid &flags,
                      MACGrid &vel,
                      Grid<Real> *density,
-                     Grid<Real> *heat = NULL,
-                     Grid<Real> *fuel = NULL,
-                     Grid<Real> *flame = NULL,
-                     Grid<Real> *red = NULL,
-                     Grid<Real> *green = NULL,
-                     Grid<Real> *blue = NULL,
+                     Grid<Real> *heat = nullptr,
+                     Grid<Real> *fuel = nullptr,
+                     Grid<Real> *flame = nullptr,
+                     Grid<Real> *red = nullptr,
+                     Grid<Real> *green = nullptr,
+                     Grid<Real> *blue = nullptr,
                      Real resetValue = 0)
 {
   KnResetInObstacle(flags, vel, density, heat, fuel, flame, red, green, blue, resetValue);
@@ -1178,18 +1178,18 @@ static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "resetInObstacle", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
       Grid<Real> *density = _args.getPtr<Grid<Real>>("density", 2, &_lock);
-      Grid<Real> *heat = _args.getPtrOpt<Grid<Real>>("heat", 3, NULL, &_lock);
-      Grid<Real> *fuel = _args.getPtrOpt<Grid<Real>>("fuel", 4, NULL, &_lock);
-      Grid<Real> *flame = _args.getPtrOpt<Grid<Real>>("flame", 5, NULL, &_lock);
-      Grid<Real> *red = _args.getPtrOpt<Grid<Real>>("red", 6, NULL, &_lock);
-      Grid<Real> *green = _args.getPtrOpt<Grid<Real>>("green", 7, NULL, &_lock);
-      Grid<Real> *blue = _args.getPtrOpt<Grid<Real>>("blue", 8, NULL, &_lock);
+      Grid<Real> *heat = _args.getPtrOpt<Grid<Real>>("heat", 3, nullptr, &_lock);
+      Grid<Real> *fuel = _args.getPtrOpt<Grid<Real>>("fuel", 4, nullptr, &_lock);
+      Grid<Real> *flame = _args.getPtrOpt<Grid<Real>>("flame", 5, nullptr, &_lock);
+      Grid<Real> *red = _args.getPtrOpt<Grid<Real>>("red", 6, nullptr, &_lock);
+      Grid<Real> *green = _args.getPtrOpt<Grid<Real>>("green", 7, nullptr, &_lock);
+      Grid<Real> *blue = _args.getPtrOpt<Grid<Real>>("blue", 8, nullptr, &_lock);
       Real resetValue = _args.getOpt<Real>("resetValue", 9, 0, &_lock);
       _retval = getPyNone();
       resetInObstacle(flags, vel, density, heat, fuel, flame, red, green, blue, resetValue);
@@ -1216,7 +1216,7 @@ void PbRegister_resetInObstacle()
 //! check for symmetry , optionally enfore by copying
 
 void checkSymmetry(
-    Grid<Real> &a, Grid<Real> *err = NULL, bool symmetrize = false, int axis = 0, int bound = 0)
+    Grid<Real> &a, Grid<Real> *err = nullptr, bool symmetrize = false, int axis = 0, int bound = 0)
 {
   const int c = axis;
   const int s = a.getSize()[c];
@@ -1241,11 +1241,11 @@ static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "checkSymmetry", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &a = *_args.getPtr<Grid<Real>>("a", 0, &_lock);
-      Grid<Real> *err = _args.getPtrOpt<Grid<Real>>("err", 1, NULL, &_lock);
+      Grid<Real> *err = _args.getPtrOpt<Grid<Real>>("err", 1, nullptr, &_lock);
       bool symmetrize = _args.getOpt<bool>("symmetrize", 2, false, &_lock);
       int axis = _args.getOpt<int>("axis", 3, 0, &_lock);
       int bound = _args.getOpt<int>("bound", 4, 0, &_lock);
@@ -1272,7 +1272,7 @@ void PbRegister_checkSymmetry()
 //! check for symmetry , mac grid version
 
 void checkSymmetryVec3(Grid<Vec3> &a,
-                       Grid<Real> *err = NULL,
+                       Grid<Real> *err = nullptr,
                        bool symmetrize = false,
                        int axis = 0,
                        int bound = 0,
@@ -1359,11 +1359,11 @@ static PyObject *_W_12(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "checkSymmetryVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Vec3> &a = *_args.getPtr<Grid<Vec3>>("a", 0, &_lock);
-      Grid<Real> *err = _args.getPtrOpt<Grid<Real>>("err", 1, NULL, &_lock);
+      Grid<Real> *err = _args.getPtrOpt<Grid<Real>>("err", 1, nullptr, &_lock);
       bool symmetrize = _args.getOpt<bool>("symmetrize", 2, false, &_lock);
       int axis = _args.getOpt<int>("axis", 3, 0, &_lock);
       int bound = _args.getOpt<int>("bound", 4, 0, &_lock);
@@ -1407,7 +1407,7 @@ static PyObject *_W_13(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "projectPpmFull", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<Real> &val = *_args.getPtr<Grid<Real>>("val", 0, &_lock);
@@ -1453,7 +1453,7 @@ static PyObject *_W_14(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "addTestParts", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -1504,7 +1504,7 @@ static PyObject *_W_15(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "pdataMaxDiff", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const ParticleDataImpl<Real> *a = _args.getPtr<ParticleDataImpl<Real>>("a", 0, &_lock);
@@ -1539,7 +1539,7 @@ static PyObject *_W_16(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "pdataMaxDiffInt", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const ParticleDataImpl<int> *a = _args.getPtr<ParticleDataImpl<int>>("a", 0, &_lock);
@@ -1574,7 +1574,7 @@ static PyObject *_W_17(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "pdataMaxDiffVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const ParticleDataImpl<Vec3> *a = _args.getPtr<ParticleDataImpl<Vec3>>("a", 0, &_lock);
@@ -1620,7 +1620,7 @@ static PyObject *_W_18(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "calcCenterOfMass", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 0, &_lock);
@@ -1842,7 +1842,7 @@ static PyObject *_W_19(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "updateFractions", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2007,9 +2007,9 @@ struct KnUpdateFlagsObs : public KernelBase {
 //! optionally uses fill fractions for obstacle
 void setObstacleFlags(FlagGrid &flags,
                       const Grid<Real> &phiObs,
-                      const MACGrid *fractions = NULL,
-                      const Grid<Real> *phiOut = NULL,
-                      const Grid<Real> *phiIn = NULL,
+                      const MACGrid *fractions = nullptr,
+                      const Grid<Real> *phiOut = nullptr,
+                      const Grid<Real> *phiIn = nullptr,
                       int boundaryWidth = 1)
 {
   KnUpdateFlagsObs(flags, fractions, phiObs, phiOut, phiIn, boundaryWidth);
@@ -2021,14 +2021,14 @@ static PyObject *_W_20(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setObstacleFlags", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       const Grid<Real> &phiObs = *_args.getPtr<Grid<Real>>("phiObs", 1, &_lock);
-      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 2, NULL, &_lock);
-      const Grid<Real> *phiOut = _args.getPtrOpt<Grid<Real>>("phiOut", 3, NULL, &_lock);
-      const Grid<Real> *phiIn = _args.getPtrOpt<Grid<Real>>("phiIn", 4, NULL, &_lock);
+      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 2, nullptr, &_lock);
+      const Grid<Real> *phiOut = _args.getPtrOpt<Grid<Real>>("phiOut", 3, nullptr, &_lock);
+      const Grid<Real> *phiIn = _args.getPtrOpt<Grid<Real>>("phiIn", 4, nullptr, &_lock);
       int boundaryWidth = _args.getOpt<int>("boundaryWidth", 5, 1, &_lock);
       _retval = getPyNone();
       setObstacleFlags(flags, phiObs, fractions, phiOut, phiIn, boundaryWidth);
@@ -2166,7 +2166,7 @@ static PyObject *_W_21(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "initVortexVelocity", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<Real> &phiObs = *_args.getPtr<Grid<Real>>("phiObs", 0, &_lock);
@@ -2203,10 +2203,10 @@ struct GaussianKernelCreator {
   int mDim;
   float *mMat1D;
 
-  GaussianKernelCreator() : mSigma(0.0f), mDim(0), mMat1D(NULL)
+  GaussianKernelCreator() : mSigma(0.0f), mDim(0), mMat1D(nullptr)
   {
   }
-  GaussianKernelCreator(float sigma, int dim = 0) : mSigma(0.0f), mDim(0), mMat1D(NULL)
+  GaussianKernelCreator(float sigma, int dim = 0) : mSigma(0.0f), mDim(0), mMat1D(nullptr)
   {
     setGaussianSigma(sigma, dim);
   }
@@ -2518,7 +2518,7 @@ static PyObject *_W_22(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "blurMacGrid", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &oG = *_args.getPtr<MACGrid>("oG", 0, &_lock);
@@ -2554,7 +2554,7 @@ static PyObject *_W_23(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "blurRealGrid", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &oG = *_args.getPtr<Grid<Real>>("oG", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/kepsilon.cpp b/extern/mantaflow/preprocessed/plugin/kepsilon.cpp
index 306db9e20cc..32425a5756e 100644
--- a/extern/mantaflow/preprocessed/plugin/kepsilon.cpp
+++ b/extern/mantaflow/preprocessed/plugin/kepsilon.cpp
@@ -305,7 +305,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "KEpsilonComputeProduction", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
@@ -420,7 +420,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "KEpsilonSources", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &k = *_args.getPtr<Grid<Real>>("k", 0, &_lock);
@@ -469,7 +469,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "KEpsilonBcs", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -547,7 +547,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "KEpsilonGradientDiffusion", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &k = *_args.getPtr<Grid<Real>>("k", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/meshplugins.cpp b/extern/mantaflow/preprocessed/plugin/meshplugins.cpp
index 415bca153d0..043660db20b 100644
--- a/extern/mantaflow/preprocessed/plugin/meshplugins.cpp
+++ b/extern/mantaflow/preprocessed/plugin/meshplugins.cpp
@@ -118,7 +118,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "smoothMesh", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
@@ -657,7 +657,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "subdivideMesh", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
@@ -752,7 +752,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "killSmallComponents", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/pressure.cpp b/extern/mantaflow/preprocessed/plugin/pressure.cpp
index 780ba44a2b5..dfba8e0082b 100644
--- a/extern/mantaflow/preprocessed/plugin/pressure.cpp
+++ b/extern/mantaflow/preprocessed/plugin/pressure.cpp
@@ -912,7 +912,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "releaseMG", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       FluidSolver *solver = _args.getPtrOpt<FluidSolver>("solver", 0, nullptr, &_lock);
@@ -949,10 +949,10 @@ void computePressureRhs(Grid<Real> &rhs,
                         const Grid<Real> &pressure,
                         const FlagGrid &flags,
                         Real cgAccuracy = 1e-3,
-                        const Grid<Real> *phi = 0,
-                        const Grid<Real> *perCellCorr = 0,
-                        const MACGrid *fractions = 0,
-                        const MACGrid *obvel = 0,
+                        const Grid<Real> *phi = nullptr,
+                        const Grid<Real> *perCellCorr = nullptr,
+                        const MACGrid *fractions = nullptr,
+                        const MACGrid *obvel = nullptr,
                         Real gfClamp = 1e-04,
                         Real cgMaxIterFac = 1.5,
                         bool precondition = true,
@@ -960,7 +960,7 @@ void computePressureRhs(Grid<Real> &rhs,
                         bool enforceCompatibility = false,
                         bool useL2Norm = false,
                         bool zeroPressureFixing = false,
-                        const Grid<Real> *curv = NULL,
+                        const Grid<Real> *curv = nullptr,
                         const Real surfTens = 0.)
 {
   // compute divergence and init right hand side
@@ -977,7 +977,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "computePressureRhs", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &rhs = *_args.getPtr<Grid<Real>>("rhs", 0, &_lock);
@@ -985,10 +985,11 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 2, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock);
       Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 4, 1e-3, &_lock);
-      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 5, 0, &_lock);
-      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 6, 0, &_lock);
-      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 7, 0, &_lock);
-      const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 8, 0, &_lock);
+      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 5, nullptr, &_lock);
+      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>(
+          "perCellCorr", 6, nullptr, &_lock);
+      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 7, nullptr, &_lock);
+      const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 8, nullptr, &_lock);
       Real gfClamp = _args.getOpt<Real>("gfClamp", 9, 1e-04, &_lock);
       Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 10, 1.5, &_lock);
       bool precondition = _args.getOpt<bool>("precondition", 11, true, &_lock);
@@ -996,7 +997,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       bool enforceCompatibility = _args.getOpt<bool>("enforceCompatibility", 13, false, &_lock);
       bool useL2Norm = _args.getOpt<bool>("useL2Norm", 14, false, &_lock);
       bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 15, false, &_lock);
-      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 16, NULL, &_lock);
+      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 16, nullptr, &_lock);
       const Real surfTens = _args.getOpt<Real>("surfTens", 17, 0., &_lock);
       _retval = getPyNone();
       computePressureRhs(rhs,
@@ -1050,9 +1051,9 @@ void solvePressureSystem(Grid<Real> &rhs,
                          Grid<Real> &pressure,
                          const FlagGrid &flags,
                          Real cgAccuracy = 1e-3,
-                         const Grid<Real> *phi = 0,
-                         const Grid<Real> *perCellCorr = 0,
-                         const MACGrid *fractions = 0,
+                         const Grid<Real> *phi = nullptr,
+                         const Grid<Real> *perCellCorr = nullptr,
+                         const MACGrid *fractions = nullptr,
                          Real gfClamp = 1e-04,
                          Real cgMaxIterFac = 1.5,
                          bool precondition = true,
@@ -1060,7 +1061,7 @@ void solvePressureSystem(Grid<Real> &rhs,
                          const bool enforceCompatibility = false,
                          const bool useL2Norm = false,
                          const bool zeroPressureFixing = false,
-                         const Grid<Real> *curv = NULL,
+                         const Grid<Real> *curv = nullptr,
                          const Real surfTens = 0.)
 {
   if (precondition == false)
@@ -1221,7 +1222,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "solvePressureSystem", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &rhs = *_args.getPtr<Grid<Real>>("rhs", 0, &_lock);
@@ -1229,9 +1230,10 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 2, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock);
       Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 4, 1e-3, &_lock);
-      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 5, 0, &_lock);
-      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 6, 0, &_lock);
-      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 7, 0, &_lock);
+      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 5, nullptr, &_lock);
+      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>(
+          "perCellCorr", 6, nullptr, &_lock);
+      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 7, nullptr, &_lock);
       Real gfClamp = _args.getOpt<Real>("gfClamp", 8, 1e-04, &_lock);
       Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 9, 1.5, &_lock);
       bool precondition = _args.getOpt<bool>("precondition", 10, true, &_lock);
@@ -1240,7 +1242,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
           "enforceCompatibility", 12, false, &_lock);
       const bool useL2Norm = _args.getOpt<bool>("useL2Norm", 13, false, &_lock);
       const bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 14, false, &_lock);
-      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 15, NULL, &_lock);
+      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 15, nullptr, &_lock);
       const Real surfTens = _args.getOpt<Real>("surfTens", 16, 0., &_lock);
       _retval = getPyNone();
       solvePressureSystem(rhs,
@@ -1284,9 +1286,9 @@ void correctVelocity(MACGrid &vel,
                      Grid<Real> &pressure,
                      const FlagGrid &flags,
                      Real cgAccuracy = 1e-3,
-                     const Grid<Real> *phi = 0,
-                     const Grid<Real> *perCellCorr = 0,
-                     const MACGrid *fractions = 0,
+                     const Grid<Real> *phi = nullptr,
+                     const Grid<Real> *perCellCorr = nullptr,
+                     const MACGrid *fractions = nullptr,
                      Real gfClamp = 1e-04,
                      Real cgMaxIterFac = 1.5,
                      bool precondition = true,
@@ -1294,7 +1296,7 @@ void correctVelocity(MACGrid &vel,
                      bool enforceCompatibility = false,
                      bool useL2Norm = false,
                      bool zeroPressureFixing = false,
-                     const Grid<Real> *curv = NULL,
+                     const Grid<Real> *curv = nullptr,
                      const Real surfTens = 0.)
 {
   knCorrectVelocity(flags, vel, pressure);
@@ -1311,16 +1313,17 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "correctVelocity", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
       Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 1, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
       Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 3, 1e-3, &_lock);
-      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 4, 0, &_lock);
-      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 5, 0, &_lock);
-      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 6, 0, &_lock);
+      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 4, nullptr, &_lock);
+      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>(
+          "perCellCorr", 5, nullptr, &_lock);
+      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 6, nullptr, &_lock);
       Real gfClamp = _args.getOpt<Real>("gfClamp", 7, 1e-04, &_lock);
       Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 8, 1.5, &_lock);
       bool precondition = _args.getOpt<bool>("precondition", 9, true, &_lock);
@@ -1328,7 +1331,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       bool enforceCompatibility = _args.getOpt<bool>("enforceCompatibility", 11, false, &_lock);
       bool useL2Norm = _args.getOpt<bool>("useL2Norm", 12, false, &_lock);
       bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 13, false, &_lock);
-      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 14, NULL, &_lock);
+      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 14, nullptr, &_lock);
       const Real surfTens = _args.getOpt<Real>("surfTens", 15, 0., &_lock);
       _retval = getPyNone();
       correctVelocity(vel,
@@ -1372,10 +1375,10 @@ void solvePressure(MACGrid &vel,
                    Grid<Real> &pressure,
                    const FlagGrid &flags,
                    Real cgAccuracy = 1e-3,
-                   const Grid<Real> *phi = 0,
-                   const Grid<Real> *perCellCorr = 0,
-                   const MACGrid *fractions = 0,
-                   const MACGrid *obvel = 0,
+                   const Grid<Real> *phi = nullptr,
+                   const Grid<Real> *perCellCorr = nullptr,
+                   const MACGrid *fractions = nullptr,
+                   const MACGrid *obvel = nullptr,
                    Real gfClamp = 1e-04,
                    Real cgMaxIterFac = 1.5,
                    bool precondition = true,
@@ -1383,9 +1386,9 @@ void solvePressure(MACGrid &vel,
                    bool enforceCompatibility = false,
                    bool useL2Norm = false,
                    bool zeroPressureFixing = false,
-                   const Grid<Real> *curv = NULL,
+                   const Grid<Real> *curv = nullptr,
                    const Real surfTens = 0.,
-                   Grid<Real> *retRhs = NULL)
+                   Grid<Real> *retRhs = nullptr)
 {
   Grid<Real> rhs(vel.getParent());
 
@@ -1455,17 +1458,18 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "solvePressure", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
       Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 1, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
       Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 3, 1e-3, &_lock);
-      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 4, 0, &_lock);
-      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 5, 0, &_lock);
-      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 6, 0, &_lock);
-      const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 7, 0, &_lock);
+      const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 4, nullptr, &_lock);
+      const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>(
+          "perCellCorr", 5, nullptr, &_lock);
+      const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 6, nullptr, &_lock);
+      const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 7, nullptr, &_lock);
       Real gfClamp = _args.getOpt<Real>("gfClamp", 8, 1e-04, &_lock);
       Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 9, 1.5, &_lock);
       bool precondition = _args.getOpt<bool>("precondition", 10, true, &_lock);
@@ -1473,9 +1477,9 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       bool enforceCompatibility = _args.getOpt<bool>("enforceCompatibility", 12, false, &_lock);
       bool useL2Norm = _args.getOpt<bool>("useL2Norm", 13, false, &_lock);
       bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 14, false, &_lock);
-      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 15, NULL, &_lock);
+      const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 15, nullptr, &_lock);
       const Real surfTens = _args.getOpt<Real>("surfTens", 16, 0., &_lock);
-      Grid<Real> *retRhs = _args.getPtrOpt<Grid<Real>>("retRhs", 17, NULL, &_lock);
+      Grid<Real> *retRhs = _args.getPtrOpt<Grid<Real>>("retRhs", 17, nullptr, &_lock);
       _retval = getPyNone();
       solvePressure(vel,
                     pressure,
diff --git a/extern/mantaflow/preprocessed/plugin/ptsplugins.cpp b/extern/mantaflow/preprocessed/plugin/ptsplugins.cpp
index a6bbccc5966..7b2aedb694e 100644
--- a/extern/mantaflow/preprocessed/plugin/ptsplugins.cpp
+++ b/extern/mantaflow/preprocessed/plugin/ptsplugins.cpp
@@ -98,7 +98,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "addForcePvel", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       ParticleDataImpl<Vec3> &vel = *_args.getPtr<ParticleDataImpl<Vec3>>("vel", 0, &_lock);
@@ -227,7 +227,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "updateVelocityFromDeltaPos", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -342,7 +342,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "eulerStep", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
@@ -470,7 +470,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setPartType", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp b/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp
index 5aae7285f1a..82afe2d7ab2 100644
--- a/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp
+++ b/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp
@@ -392,7 +392,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipComputeSecondaryParticlePotentials", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &potTA = *_args.getPtr<Grid<Real>>("potTA", 0, &_lock);
@@ -1041,7 +1041,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipSampleSecondaryParticles", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const std::string mode = _args.get<std::string>("mode", 0, &_lock);
@@ -1760,7 +1760,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipUpdateSecondaryParticles", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const std::string mode = _args.get<std::string>("mode", 0, &_lock);
@@ -1895,7 +1895,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipDeleteParticlesInObstacle", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       BasicParticleSystem &pts = *_args.getPtr<BasicParticleSystem>("pts", 0, &_lock);
@@ -1963,7 +1963,7 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "debugGridInfo", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2069,7 +2069,7 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setFlagsFromLevelset", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2170,7 +2170,7 @@ static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "setMACFromLevelset", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &v = *_args.getPtr<MACGrid>("v", 0, &_lock);
@@ -2382,7 +2382,7 @@ static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipComputePotentialTrappedAir", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &pot = *_args.getPtr<Grid<Real>>("pot", 0, &_lock);
@@ -2555,7 +2555,7 @@ static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipComputePotentialKineticEnergy", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &pot = *_args.getPtr<Grid<Real>>("pot", 0, &_lock);
@@ -2791,7 +2791,7 @@ static PyObject *_W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipComputePotentialWaveCrest", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &pot = *_args.getPtr<Grid<Real>>("pot", 0, &_lock);
@@ -2883,7 +2883,7 @@ static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipComputeSurfaceNormals", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Vec3> &normal = *_args.getPtr<Grid<Vec3>>("normal", 0, &_lock);
@@ -3040,7 +3040,7 @@ static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "flipUpdateNeighborRatio", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/surfaceturbulence.cpp b/extern/mantaflow/preprocessed/plugin/surfaceturbulence.cpp
index 465314f51ed..aa04b551e1e 100644
--- a/extern/mantaflow/preprocessed/plugin/surfaceturbulence.cpp
+++ b/extern/mantaflow/preprocessed/plugin/surfaceturbulence.cpp
@@ -2055,7 +2055,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "particleSurfaceTurbulence", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -2159,7 +2159,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "debugCheckParts", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp b/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp
index c2a21d82689..d5d7d597a7c 100644
--- a/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp
+++ b/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp
@@ -48,7 +48,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "markAsFixed", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
@@ -111,7 +111,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "texcoordInflow", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
@@ -155,7 +155,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "meshSmokeInflow", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
@@ -239,8 +239,8 @@ struct KnAcceleration : public KernelBase {
 
 void vorticitySource(VortexSheetMesh &mesh,
                      Vec3 gravity,
-                     const MACGrid *vel = NULL,
-                     const MACGrid *velOld = NULL,
+                     const MACGrid *vel = nullptr,
+                     const MACGrid *velOld = nullptr,
                      Real scale = 0.1,
                      Real maxAmount = 0,
                      Real mult = 1.0)
@@ -289,13 +289,13 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "vorticitySource", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
       Vec3 gravity = _args.get<Vec3>("gravity", 1, &_lock);
-      const MACGrid *vel = _args.getPtrOpt<MACGrid>("vel", 2, NULL, &_lock);
-      const MACGrid *velOld = _args.getPtrOpt<MACGrid>("velOld", 3, NULL, &_lock);
+      const MACGrid *vel = _args.getPtrOpt<MACGrid>("vel", 2, nullptr, &_lock);
+      const MACGrid *velOld = _args.getPtrOpt<MACGrid>("velOld", 3, nullptr, &_lock);
       Real scale = _args.getOpt<Real>("scale", 4, 0.1, &_lock);
       Real maxAmount = _args.getOpt<Real>("maxAmount", 5, 0, &_lock);
       Real mult = _args.getOpt<Real>("mult", 6, 1.0, &_lock);
@@ -373,7 +373,7 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "smoothVorticity", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
@@ -437,7 +437,7 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "VPseedK41", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       VortexParticleSystem &system = *_args.getPtr<VortexParticleSystem>("system", 0, &_lock);
@@ -473,7 +473,7 @@ void VICintegration(VortexSheetMesh &mesh,
                     Real sigma,
                     Grid<Vec3> &vel,
                     const FlagGrid &flags,
-                    Grid<Vec3> *vorticity = NULL,
+                    Grid<Vec3> *vorticity = nullptr,
                     Real cgMaxIterFac = 1.5,
                     Real cgAccuracy = 1e-3,
                     Real scale = 0.01,
@@ -604,14 +604,14 @@ static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "VICintegration", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
       Real sigma = _args.get<Real>("sigma", 1, &_lock);
       Grid<Vec3> &vel = *_args.getPtr<Grid<Vec3>>("vel", 2, &_lock);
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock);
-      Grid<Vec3> *vorticity = _args.getPtrOpt<Grid<Vec3>>("vorticity", 4, NULL, &_lock);
+      Grid<Vec3> *vorticity = _args.getPtrOpt<Grid<Vec3>>("vorticity", 4, nullptr, &_lock);
       Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 5, 1.5, &_lock);
       Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 6, 1e-3, &_lock);
       Real scale = _args.getOpt<Real>("scale", 7, 0.01, &_lock);
@@ -665,7 +665,7 @@ static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "densityFromLevelset", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/waveletturbulence.cpp b/extern/mantaflow/preprocessed/plugin/waveletturbulence.cpp
index 9d3bdaa3f21..7d867542132 100644
--- a/extern/mantaflow/preprocessed/plugin/waveletturbulence.cpp
+++ b/extern/mantaflow/preprocessed/plugin/waveletturbulence.cpp
@@ -69,7 +69,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "interpolateGrid", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &target = *_args.getPtr<Grid<Real>>("target", 0, &_lock);
@@ -116,7 +116,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "interpolateGridVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Vec3> &target = *_args.getPtr<Grid<Vec3>>("target", 0, &_lock);
@@ -263,7 +263,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "interpolateMACGrid", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       MACGrid &target = *_args.getPtr<MACGrid>("target", 0, &_lock);
@@ -395,7 +395,7 @@ void applySimpleNoiseVec3(const FlagGrid &flags,
                           Grid<Vec3> &target,
                           const WaveletNoiseField &noise,
                           Real scale = 1.0,
-                          const Grid<Real> *weight = NULL)
+                          const Grid<Real> *weight = nullptr)
 {
   // note - passing a MAC grid here is slightly inaccurate, we should evaluate each component
   // separately
@@ -408,14 +408,14 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "applySimpleNoiseVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       Grid<Vec3> &target = *_args.getPtr<Grid<Vec3>>("target", 1, &_lock);
       const WaveletNoiseField &noise = *_args.getPtr<WaveletNoiseField>("noise", 2, &_lock);
       Real scale = _args.getOpt<Real>("scale", 3, 1.0, &_lock);
-      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 4, NULL, &_lock);
+      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 4, nullptr, &_lock);
       _retval = getPyNone();
       applySimpleNoiseVec3(flags, target, noise, scale, weight);
       _args.check();
@@ -537,7 +537,7 @@ void applySimpleNoiseReal(const FlagGrid &flags,
                           Grid<Real> &target,
                           const WaveletNoiseField &noise,
                           Real scale = 1.0,
-                          const Grid<Real> *weight = NULL)
+                          const Grid<Real> *weight = nullptr)
 {
   knApplySimpleNoiseReal(flags, target, noise, scale, weight);
 }
@@ -548,14 +548,14 @@ static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "applySimpleNoiseReal", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
       Grid<Real> &target = *_args.getPtr<Grid<Real>>("target", 1, &_lock);
       const WaveletNoiseField &noise = *_args.getPtr<WaveletNoiseField>("noise", 2, &_lock);
       Real scale = _args.getOpt<Real>("scale", 3, 1.0, &_lock);
-      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 4, NULL, &_lock);
+      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 4, nullptr, &_lock);
       _retval = getPyNone();
       applySimpleNoiseReal(flags, target, noise, scale, weight);
       _args.check();
@@ -763,8 +763,8 @@ void applyNoiseVec3(const FlagGrid &flags,
                     const WaveletNoiseField &noise,
                     Real scale = 1.0,
                     Real scaleSpatial = 1.0,
-                    const Grid<Real> *weight = NULL,
-                    const Grid<Vec3> *uv = NULL)
+                    const Grid<Real> *weight = nullptr,
+                    const Grid<Vec3> *uv = nullptr)
 {
   // check whether the uv grid has a different resolution
   bool uvInterpol = false;
@@ -794,7 +794,7 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "applyNoiseVec3", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -802,8 +802,8 @@ static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const WaveletNoiseField &noise = *_args.getPtr<WaveletNoiseField>("noise", 2, &_lock);
       Real scale = _args.getOpt<Real>("scale", 3, 1.0, &_lock);
       Real scaleSpatial = _args.getOpt<Real>("scaleSpatial", 4, 1.0, &_lock);
-      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 5, NULL, &_lock);
-      const Grid<Vec3> *uv = _args.getPtrOpt<Grid<Vec3>>("uv", 6, NULL, &_lock);
+      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 5, nullptr, &_lock);
+      const Grid<Vec3> *uv = _args.getPtrOpt<Grid<Vec3>>("uv", 6, nullptr, &_lock);
       _retval = getPyNone();
       applyNoiseVec3(flags, target, noise, scale, scaleSpatial, weight, uv);
       _args.check();
@@ -905,7 +905,7 @@ static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "computeEnergy", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -943,7 +943,7 @@ static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "computeWaveletCoeffs", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &input = *_args.getPtr<Grid<Real>>("input", 0, &_lock);
@@ -968,7 +968,7 @@ void PbRegister_computeWaveletCoeffs()
 }
 
 // note - alomst the same as for vorticity confinement
-void computeVorticity(const MACGrid &vel, Grid<Vec3> &vorticity, Grid<Real> *norm = NULL)
+void computeVorticity(const MACGrid &vel, Grid<Vec3> &vorticity, Grid<Real> *norm = nullptr)
 {
   Grid<Vec3> velCenter(vel.getParent());
   GetCentered(velCenter, vel);
@@ -983,12 +983,12 @@ static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "computeVorticity", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
       Grid<Vec3> &vorticity = *_args.getPtr<Grid<Vec3>>("vorticity", 1, &_lock);
-      Grid<Real> *norm = _args.getPtrOpt<Grid<Real>>("norm", 2, NULL, &_lock);
+      Grid<Real> *norm = _args.getPtrOpt<Grid<Real>>("norm", 2, nullptr, &_lock);
       _retval = getPyNone();
       computeVorticity(vel, vorticity, norm);
       _args.check();
@@ -1104,7 +1104,7 @@ static PyObject *_W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "computeStrainRateMag", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
@@ -1219,7 +1219,7 @@ static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "extrapolateSimpleFlags", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
@@ -1263,7 +1263,7 @@ static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "getCurl", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
diff --git a/extern/mantaflow/preprocessed/plugin/waves.cpp b/extern/mantaflow/preprocessed/plugin/waves.cpp
index 7745dce4711..6fb3b16c742 100644
--- a/extern/mantaflow/preprocessed/plugin/waves.cpp
+++ b/extern/mantaflow/preprocessed/plugin/waves.cpp
@@ -101,7 +101,7 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "calcSecDeriv2d", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const Grid<Real> &v = *_args.getPtr<Grid<Real>>("v", 0, &_lock);
@@ -206,7 +206,7 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "totalSum", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &height = *_args.getPtr<Grid<Real>>("height", 0, &_lock);
@@ -243,7 +243,7 @@ static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "normalizeSumTo", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       Grid<Real> &height = *_args.getPtr<Grid<Real>>("height", 0, &_lock);
@@ -449,7 +449,7 @@ static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
     FluidSolver *parent = _args.obtainParent();
     bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
     pbPreparePlugin(parent, "cgSolveWE", !noTiming);
-    PyObject *_retval = 0;
+    PyObject *_retval = nullptr;
     {
       ArgLocker _lock;
       const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);