// DO NOT EDIT !
// This file is generated using the MantaFlow preprocessor (prep generate).

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2017 Georg Kohl, Nils Thuerey
 *
 * This program is free software, distributed under the terms of the
 * GNU General Public License (GPL)
 * http://www.gnu.org/licenses
 *
 * Secondary particle plugin for FLIP simulations
 *
 ******************************************************************************/

#include "particle.h"
#include "commonkernels.h"

namespace Manta {

#pragma region Secondary Particles for FLIP
//----------------------------------------------------------------------------------------------------------------------------------------------------
// Secondary Particles for FLIP
//----------------------------------------------------------------------------------------------------------------------------------------------------

// helper function that clamps the value in potential to the interval [tauMin, tauMax] and
// normalizes it to [0, 1] afterwards
Real clampPotential(Real potential, Real tauMin, Real tauMax)
{
  return (std::min(potential, tauMax) - std::min(potential, tauMin)) / (tauMax - tauMin);
}

// computes all three potentials(trapped air, wave crest, kinetic energy) and the neighbor ratio
// for every fluid cell and stores it in the respective grid. Is less readable but significantly
// faster than using seperate potential computation

struct knFlipComputeSecondaryParticlePotentials : public KernelBase {
  knFlipComputeSecondaryParticlePotentials(Grid<Real> &potTA,
                                           Grid<Real> &potWC,
                                           Grid<Real> &potKE,
                                           Grid<Real> &neighborRatio,
                                           const FlagGrid &flags,
                                           const MACGrid &v,
                                           const Grid<Vec3> &normal,
                                           const int radius,
                                           const Real tauMinTA,
                                           const Real tauMaxTA,
                                           const Real tauMinWC,
                                           const Real tauMaxWC,
                                           const Real tauMinKE,
                                           const Real tauMaxKE,
                                           const Real scaleFromManta,
                                           const int itype = FlagGrid::TypeFluid,
                                           const int jtype = FlagGrid::TypeObstacle |
                                                             FlagGrid::TypeOutflow |
                                                             FlagGrid::TypeInflow)
      : KernelBase(&potTA, radius),
        potTA(potTA),
        potWC(potWC),
        potKE(potKE),
        neighborRatio(neighborRatio),
        flags(flags),
        v(v),
        normal(normal),
        radius(radius),
        tauMinTA(tauMinTA),
        tauMaxTA(tauMaxTA),
        tauMinWC(tauMinWC),
        tauMaxWC(tauMaxWC),
        tauMinKE(tauMinKE),
        tauMaxKE(tauMaxKE),
        scaleFromManta(scaleFromManta),
        itype(itype),
        jtype(jtype)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 Grid<Real> &potTA,
                 Grid<Real> &potWC,
                 Grid<Real> &potKE,
                 Grid<Real> &neighborRatio,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 const Grid<Vec3> &normal,
                 const int radius,
                 const Real tauMinTA,
                 const Real tauMaxTA,
                 const Real tauMinWC,
                 const Real tauMaxWC,
                 const Real tauMinKE,
                 const Real tauMaxKE,
                 const Real scaleFromManta,
                 const int itype = FlagGrid::TypeFluid,
                 const int jtype = FlagGrid::TypeObstacle | FlagGrid::TypeOutflow |
                                   FlagGrid::TypeInflow) const
  {

    if (!(flags(i, j, k) & itype))
      return;

    // compute trapped air potential + wave crest potential + neighbor ratio at once
    const Vec3 &xi = scaleFromManta * Vec3(i, j, k);  // scale to unit cube
    const Vec3 &vi = scaleFromManta * v.getCentered(i, j, k);
    const Vec3 &ni = getNormalized(normal(i, j, k));
    Real vdiff = 0;         // for trapped air
    Real kappa = 0;         // for wave crests
    int countFluid = 0;     // for neighbor ratio
    int countMaxFluid = 0;  // for neighbor ratio

    // iterate over neighboring cells within radius
    for (IndexInt x = i - radius; x <= i + radius; x++) {
      for (IndexInt y = j - radius; y <= j + radius; y++) {
        for (IndexInt z = k - radius; z <= k + radius; z++) {
          if ((x == i && y == j && z == k) || !flags.isInBounds(Vec3i(x, y, z)) ||
              (flags(x, y, z) & jtype))
            continue;

          if (flags(x, y, z) & itype) {
            countFluid++;
            countMaxFluid++;
          }
          else {
            countMaxFluid++;
          }

          const Vec3 &xj = scaleFromManta * Vec3(x, y, z);  // scale to unit cube
          const Vec3 &vj = scaleFromManta * v.getCentered(x, y, z);
          const Vec3 &nj = getNormalized(normal(x, y, z));
          const Vec3 xij = xi - xj;
          const Vec3 vij = vi - vj;
          Real h = !potTA.is3D() ?
                       1.414 * radius :
                       1.732 * radius;  // estimate sqrt(2)*radius resp. sqrt(3)*radius for h, due
                                        // to squared resp. cubic neighbor area
          vdiff += norm(vij) * (1 - dot(getNormalized(vij), getNormalized(xij))) *
                   (1 - norm(xij) / h);

          if (dot(getNormalized(xij), ni) < 0) {  // identifies wave crests
            kappa += (1 - dot(ni, nj)) * (1 - norm(xij) / h);
          }
        }
      }
    }

    neighborRatio(i, j, k) = float(countFluid) / float(countMaxFluid);

    potTA(i, j, k) = clampPotential(vdiff, tauMinTA, tauMaxTA);
    if (dot(getNormalized(vi), ni) >= 0.6) {  // avoid to mark boarders of the scene as wave crest
      potWC(i, j, k) = clampPotential(kappa, tauMinWC, tauMaxWC);
    }
    else {
      potWC(i, j, k) = Real(0);
    }

    // compute kinetic energy potential
    Real ek =
        Real(0.5) * 125 *
        normSquare(
            vi);  // use arbitrary constant for mass, potential adjusts with thresholds anyways
    potKE(i, j, k) = clampPotential(ek, tauMinKE, tauMaxKE);
  }
  inline Grid<Real> &getArg0()
  {
    return potTA;
  }
  typedef Grid<Real> type0;
  inline Grid<Real> &getArg1()
  {
    return potWC;
  }
  typedef Grid<Real> type1;
  inline Grid<Real> &getArg2()
  {
    return potKE;
  }
  typedef Grid<Real> type2;
  inline Grid<Real> &getArg3()
  {
    return neighborRatio;
  }
  typedef Grid<Real> type3;
  inline const FlagGrid &getArg4()
  {
    return flags;
  }
  typedef FlagGrid type4;
  inline const MACGrid &getArg5()
  {
    return v;
  }
  typedef MACGrid type5;
  inline const Grid<Vec3> &getArg6()
  {
    return normal;
  }
  typedef Grid<Vec3> type6;
  inline const int &getArg7()
  {
    return radius;
  }
  typedef int type7;
  inline const Real &getArg8()
  {
    return tauMinTA;
  }
  typedef Real type8;
  inline const Real &getArg9()
  {
    return tauMaxTA;
  }
  typedef Real type9;
  inline const Real &getArg10()
  {
    return tauMinWC;
  }
  typedef Real type10;
  inline const Real &getArg11()
  {
    return tauMaxWC;
  }
  typedef Real type11;
  inline const Real &getArg12()
  {
    return tauMinKE;
  }
  typedef Real type12;
  inline const Real &getArg13()
  {
    return tauMaxKE;
  }
  typedef Real type13;
  inline const Real &getArg14()
  {
    return scaleFromManta;
  }
  typedef Real type14;
  inline const int &getArg15()
  {
    return itype;
  }
  typedef int type15;
  inline const int &getArg16()
  {
    return jtype;
  }
  typedef int type16;
  void runMessage()
  {
    debMsg("Executing kernel knFlipComputeSecondaryParticlePotentials ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = radius; j < _maxY; j++)
          for (int i = radius; i < _maxX; i++)
            op(i,
               j,
               k,
               potTA,
               potWC,
               potKE,
               neighborRatio,
               flags,
               v,
               normal,
               radius,
               tauMinTA,
               tauMaxTA,
               tauMinWC,
               tauMaxWC,
               tauMinKE,
               tauMaxKE,
               scaleFromManta,
               itype,
               jtype);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = radius; i < _maxX; i++)
          op(i,
             j,
             k,
             potTA,
             potWC,
             potKE,
             neighborRatio,
             flags,
             v,
             normal,
             radius,
             tauMinTA,
             tauMaxTA,
             tauMinWC,
             tauMaxWC,
             tauMinKE,
             tauMaxKE,
             scaleFromManta,
             itype,
             jtype);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(radius, maxY), *this);
  }
  Grid<Real> &potTA;
  Grid<Real> &potWC;
  Grid<Real> &potKE;
  Grid<Real> &neighborRatio;
  const FlagGrid &flags;
  const MACGrid &v;
  const Grid<Vec3> &normal;
  const int radius;
  const Real tauMinTA;
  const Real tauMaxTA;
  const Real tauMinWC;
  const Real tauMaxWC;
  const Real tauMinKE;
  const Real tauMaxKE;
  const Real scaleFromManta;
  const int itype;
  const int jtype;
};

void flipComputeSecondaryParticlePotentials(Grid<Real> &potTA,
                                            Grid<Real> &potWC,
                                            Grid<Real> &potKE,
                                            Grid<Real> &neighborRatio,
                                            const FlagGrid &flags,
                                            const MACGrid &v,
                                            Grid<Vec3> &normal,
                                            const Grid<Real> &phi,
                                            const int radius,
                                            const Real tauMinTA,
                                            const Real tauMaxTA,
                                            const Real tauMinWC,
                                            const Real tauMaxWC,
                                            const Real tauMinKE,
                                            const Real tauMaxKE,
                                            const Real scaleFromManta,
                                            const int itype = FlagGrid::TypeFluid,
                                            const int jtype = FlagGrid::TypeObstacle |
                                                              FlagGrid::TypeOutflow |
                                                              FlagGrid::TypeInflow)
{
  potTA.clear();
  potWC.clear();
  potKE.clear();
  neighborRatio.clear();
  GradientOp(normal, phi);
  knFlipComputeSecondaryParticlePotentials(potTA,
                                           potWC,
                                           potKE,
                                           neighborRatio,
                                           flags,
                                           v,
                                           normal,
                                           radius,
                                           tauMinTA,
                                           tauMaxTA,
                                           tauMinWC,
                                           tauMaxWC,
                                           tauMinKE,
                                           tauMaxKE,
                                           scaleFromManta,
                                           itype,
                                           jtype);
}
static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipComputeSecondaryParticlePotentials", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &potTA = *_args.getPtr<Grid<Real>>("potTA", 0, &_lock);
      Grid<Real> &potWC = *_args.getPtr<Grid<Real>>("potWC", 1, &_lock);
      Grid<Real> &potKE = *_args.getPtr<Grid<Real>>("potKE", 2, &_lock);
      Grid<Real> &neighborRatio = *_args.getPtr<Grid<Real>>("neighborRatio", 3, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 4, &_lock);
      const MACGrid &v = *_args.getPtr<MACGrid>("v", 5, &_lock);
      Grid<Vec3> &normal = *_args.getPtr<Grid<Vec3>>("normal", 6, &_lock);
      const Grid<Real> &phi = *_args.getPtr<Grid<Real>>("phi", 7, &_lock);
      const int radius = _args.get<int>("radius", 8, &_lock);
      const Real tauMinTA = _args.get<Real>("tauMinTA", 9, &_lock);
      const Real tauMaxTA = _args.get<Real>("tauMaxTA", 10, &_lock);
      const Real tauMinWC = _args.get<Real>("tauMinWC", 11, &_lock);
      const Real tauMaxWC = _args.get<Real>("tauMaxWC", 12, &_lock);
      const Real tauMinKE = _args.get<Real>("tauMinKE", 13, &_lock);
      const Real tauMaxKE = _args.get<Real>("tauMaxKE", 14, &_lock);
      const Real scaleFromManta = _args.get<Real>("scaleFromManta", 15, &_lock);
      const int itype = _args.getOpt<int>("itype", 16, FlagGrid::TypeFluid, &_lock);
      const int jtype = _args.getOpt<int>("jtype",
                                          17,
                                          FlagGrid::TypeObstacle | FlagGrid::TypeOutflow |
                                              FlagGrid::TypeInflow,
                                          &_lock);
      _retval = getPyNone();
      flipComputeSecondaryParticlePotentials(potTA,
                                             potWC,
                                             potKE,
                                             neighborRatio,
                                             flags,
                                             v,
                                             normal,
                                             phi,
                                             radius,
                                             tauMinTA,
                                             tauMaxTA,
                                             tauMinWC,
                                             tauMaxWC,
                                             tauMinKE,
                                             tauMaxKE,
                                             scaleFromManta,
                                             itype,
                                             jtype);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipComputeSecondaryParticlePotentials", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipComputeSecondaryParticlePotentials", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipComputeSecondaryParticlePotentials(
    "", "flipComputeSecondaryParticlePotentials", _W_0);
extern "C" {
void PbRegister_flipComputeSecondaryParticlePotentials()
{
  KEEP_UNUSED(_RP_flipComputeSecondaryParticlePotentials);
}
}

// adds secondary particles to &pts_sec for every fluid cell in &flags according to the potential
// grids &potTA, &potWC and &potKE secondary particles are uniformly sampled in every fluid cell in
// a randomly offset cylinder in fluid movement direction In contrast to
// flipSampleSecondaryParticles this uses more cylinders per cell and interpolates velocity and
// potentials. To control number of cylinders in each dimension adjust radius(0.25=>2 cyl,
// 0.1666=>3 cyl, 0.125=>3cyl etc.).

struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
  knFlipSampleSecondaryParticlesMoreCylinders(const FlagGrid &flags,
                                              const MACGrid &v,
                                              BasicParticleSystem &pts_sec,
                                              ParticleDataImpl<Vec3> &v_sec,
                                              ParticleDataImpl<Real> &l_sec,
                                              const Real lMin,
                                              const Real lMax,
                                              const Grid<Real> &potTA,
                                              const Grid<Real> &potWC,
                                              const Grid<Real> &potKE,
                                              const Grid<Real> &neighborRatio,
                                              const Real c_s,
                                              const Real c_b,
                                              const Real k_ta,
                                              const Real k_wc,
                                              const Real dt,
                                              const int itype = FlagGrid::TypeFluid)
      : KernelBase(&flags, 0),
        flags(flags),
        v(v),
        pts_sec(pts_sec),
        v_sec(v_sec),
        l_sec(l_sec),
        lMin(lMin),
        lMax(lMax),
        potTA(potTA),
        potWC(potWC),
        potKE(potKE),
        neighborRatio(neighborRatio),
        c_s(c_s),
        c_b(c_b),
        k_ta(k_ta),
        k_wc(k_wc),
        dt(dt),
        itype(itype)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 BasicParticleSystem &pts_sec,
                 ParticleDataImpl<Vec3> &v_sec,
                 ParticleDataImpl<Real> &l_sec,
                 const Real lMin,
                 const Real lMax,
                 const Grid<Real> &potTA,
                 const Grid<Real> &potWC,
                 const Grid<Real> &potKE,
                 const Grid<Real> &neighborRatio,
                 const Real c_s,
                 const Real c_b,
                 const Real k_ta,
                 const Real k_wc,
                 const Real dt,
                 const int itype = FlagGrid::TypeFluid)
  {

    if (!(flags(i, j, k) & itype))
      return;

    RandomStream mRand(9832);
    Real radius =
        0.25;  // diameter=0.5 => sampling with two cylinders in each dimension since cell size=1
    for (Real x = i - radius; x <= i + radius; x += 2 * radius) {
      for (Real y = j - radius; y <= j + radius; y += 2 * radius) {
        for (Real z = k - radius; z <= k + radius; z += 2 * radius) {

          Vec3 xi = Vec3(x, y, z);
          Real KE = potKE.getInterpolated(xi);
          Real TA = potTA.getInterpolated(xi);
          Real WC = potWC.getInterpolated(xi);

          const int n = KE * (k_ta * TA + k_wc * WC) * dt;  // number of secondary particles
          if (n == 0)
            continue;
          Vec3 vi = v.getInterpolated(xi);
          Vec3 dir = dt * vi;  // direction of movement of current particle
          Vec3 e1 = getNormalized(Vec3(dir.z, 0, -dir.x));  // perpendicular to dir
          Vec3 e2 = getNormalized(
              cross(e1, dir));  // perpendicular to dir and e1, so e1 and e1 create reference plane

          for (int di = 0; di < n; di++) {
            const Real r = radius * sqrt(mRand.getReal());        // distance to cylinder axis
            const Real theta = mRand.getReal() * Real(2) * M_PI;  // azimuth
            const Real h = mRand.getReal() * norm(dt * vi);       // distance to reference plane
            Vec3 xd = xi + r * cos(theta) * e1 + r * sin(theta) * e2 + h * getNormalized(vi);
            if (!flags.is3D())
              xd.z = 0;
            pts_sec.add(xd);

            v_sec[v_sec.size() - 1] = r * cos(theta) * e1 + r * sin(theta) * e2 +
                                      vi;  // init velocity of new particle
            Real temp = (KE + TA + WC) / 3;
            l_sec[l_sec.size() - 1] = ((lMax - lMin) * temp) + lMin +
                                      mRand.getReal() * 0.1;  // init lifetime of new particle

            // init type of new particle
            if (neighborRatio(i, j, k) < c_s) {
              pts_sec[pts_sec.size() - 1].flag = ParticleBase::PSPRAY;
            }
            else if (neighborRatio(i, j, k) > c_b) {
              pts_sec[pts_sec.size() - 1].flag = ParticleBase::PBUBBLE;
            }
            else {
              pts_sec[pts_sec.size() - 1].flag = ParticleBase::PFOAM;
            }
          }
        }
      }
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return v;
  }
  typedef MACGrid type1;
  inline BasicParticleSystem &getArg2()
  {
    return pts_sec;
  }
  typedef BasicParticleSystem type2;
  inline ParticleDataImpl<Vec3> &getArg3()
  {
    return v_sec;
  }
  typedef ParticleDataImpl<Vec3> type3;
  inline ParticleDataImpl<Real> &getArg4()
  {
    return l_sec;
  }
  typedef ParticleDataImpl<Real> type4;
  inline const Real &getArg5()
  {
    return lMin;
  }
  typedef Real type5;
  inline const Real &getArg6()
  {
    return lMax;
  }
  typedef Real type6;
  inline const Grid<Real> &getArg7()
  {
    return potTA;
  }
  typedef Grid<Real> type7;
  inline const Grid<Real> &getArg8()
  {
    return potWC;
  }
  typedef Grid<Real> type8;
  inline const Grid<Real> &getArg9()
  {
    return potKE;
  }
  typedef Grid<Real> type9;
  inline const Grid<Real> &getArg10()
  {
    return neighborRatio;
  }
  typedef Grid<Real> type10;
  inline const Real &getArg11()
  {
    return c_s;
  }
  typedef Real type11;
  inline const Real &getArg12()
  {
    return c_b;
  }
  typedef Real type12;
  inline const Real &getArg13()
  {
    return k_ta;
  }
  typedef Real type13;
  inline const Real &getArg14()
  {
    return k_wc;
  }
  typedef Real type14;
  inline const Real &getArg15()
  {
    return dt;
  }
  typedef Real type15;
  inline const int &getArg16()
  {
    return itype;
  }
  typedef int type16;
  void runMessage()
  {
    debMsg("Executing kernel knFlipSampleSecondaryParticlesMoreCylinders ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    for (int k = minZ; k < maxZ; k++)
      for (int j = 0; j < _maxY; j++)
        for (int i = 0; i < _maxX; i++)
          op(i,
             j,
             k,
             flags,
             v,
             pts_sec,
             v_sec,
             l_sec,
             lMin,
             lMax,
             potTA,
             potWC,
             potKE,
             neighborRatio,
             c_s,
             c_b,
             k_ta,
             k_wc,
             dt,
             itype);
  }
  const FlagGrid &flags;
  const MACGrid &v;
  BasicParticleSystem &pts_sec;
  ParticleDataImpl<Vec3> &v_sec;
  ParticleDataImpl<Real> &l_sec;
  const Real lMin;
  const Real lMax;
  const Grid<Real> &potTA;
  const Grid<Real> &potWC;
  const Grid<Real> &potKE;
  const Grid<Real> &neighborRatio;
  const Real c_s;
  const Real c_b;
  const Real k_ta;
  const Real k_wc;
  const Real dt;
  const int itype;
};

// adds secondary particles to &pts_sec for every fluid cell in &flags according to the potential
// grids &potTA, &potWC and &potKE secondary particles are uniformly sampled in every fluid cell in
// a randomly offset cylinder in fluid movement direction

struct knFlipSampleSecondaryParticles : public KernelBase {
  knFlipSampleSecondaryParticles(const FlagGrid &flags,
                                 const MACGrid &v,
                                 BasicParticleSystem &pts_sec,
                                 ParticleDataImpl<Vec3> &v_sec,
                                 ParticleDataImpl<Real> &l_sec,
                                 const Real lMin,
                                 const Real lMax,
                                 const Grid<Real> &potTA,
                                 const Grid<Real> &potWC,
                                 const Grid<Real> &potKE,
                                 const Grid<Real> &neighborRatio,
                                 const Real c_s,
                                 const Real c_b,
                                 const Real k_ta,
                                 const Real k_wc,
                                 const Real dt,
                                 const int itype = FlagGrid::TypeFluid)
      : KernelBase(&flags, 0),
        flags(flags),
        v(v),
        pts_sec(pts_sec),
        v_sec(v_sec),
        l_sec(l_sec),
        lMin(lMin),
        lMax(lMax),
        potTA(potTA),
        potWC(potWC),
        potKE(potKE),
        neighborRatio(neighborRatio),
        c_s(c_s),
        c_b(c_b),
        k_ta(k_ta),
        k_wc(k_wc),
        dt(dt),
        itype(itype)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 BasicParticleSystem &pts_sec,
                 ParticleDataImpl<Vec3> &v_sec,
                 ParticleDataImpl<Real> &l_sec,
                 const Real lMin,
                 const Real lMax,
                 const Grid<Real> &potTA,
                 const Grid<Real> &potWC,
                 const Grid<Real> &potKE,
                 const Grid<Real> &neighborRatio,
                 const Real c_s,
                 const Real c_b,
                 const Real k_ta,
                 const Real k_wc,
                 const Real dt,
                 const int itype = FlagGrid::TypeFluid)
  {

    if (!(flags(i, j, k) & itype))
      return;

    Real KE = potKE(i, j, k);
    Real TA = potTA(i, j, k);
    Real WC = potWC(i, j, k);

    const int n = KE * (k_ta * TA + k_wc * WC) * dt;  // number of secondary particles
    if (n == 0)
      return;
    RandomStream mRand(9832);

    Vec3 xi = Vec3(i + mRand.getReal(),
                   j + mRand.getReal(),
                   k + mRand.getReal());  // randomized offset uniform in cell
    Vec3 vi = v.getInterpolated(xi);
    Vec3 dir = dt * vi;                               // direction of movement of current particle
    Vec3 e1 = getNormalized(Vec3(dir.z, 0, -dir.x));  // perpendicular to dir
    Vec3 e2 = getNormalized(
        cross(e1, dir));  // perpendicular to dir and e1, so e1 and e1 create reference plane

    for (int di = 0; di < n; di++) {
      const Real r = Real(0.5) * sqrt(mRand.getReal());     // distance to cylinder axis
      const Real theta = mRand.getReal() * Real(2) * M_PI;  // azimuth
      const Real h = mRand.getReal() * norm(dt * vi);       // distance to reference plane
      Vec3 xd = xi + r * cos(theta) * e1 + r * sin(theta) * e2 + h * getNormalized(vi);
      if (!flags.is3D())
        xd.z = 0;
      pts_sec.add(xd);

      v_sec[v_sec.size() - 1] = r * cos(theta) * e1 + r * sin(theta) * e2 +
                                vi;  // init velocity of new particle
      Real temp = (KE + TA + WC) / 3;
      l_sec[l_sec.size() - 1] = ((lMax - lMin) * temp) + lMin +
                                mRand.getReal() * 0.1;  // init lifetime of new particle

      // init type of new particle
      if (neighborRatio(i, j, k) < c_s) {
        pts_sec[pts_sec.size() - 1].flag = ParticleBase::PSPRAY;
      }
      else if (neighborRatio(i, j, k) > c_b) {
        pts_sec[pts_sec.size() - 1].flag = ParticleBase::PBUBBLE;
      }
      else {
        pts_sec[pts_sec.size() - 1].flag = ParticleBase::PFOAM;
      }
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return v;
  }
  typedef MACGrid type1;
  inline BasicParticleSystem &getArg2()
  {
    return pts_sec;
  }
  typedef BasicParticleSystem type2;
  inline ParticleDataImpl<Vec3> &getArg3()
  {
    return v_sec;
  }
  typedef ParticleDataImpl<Vec3> type3;
  inline ParticleDataImpl<Real> &getArg4()
  {
    return l_sec;
  }
  typedef ParticleDataImpl<Real> type4;
  inline const Real &getArg5()
  {
    return lMin;
  }
  typedef Real type5;
  inline const Real &getArg6()
  {
    return lMax;
  }
  typedef Real type6;
  inline const Grid<Real> &getArg7()
  {
    return potTA;
  }
  typedef Grid<Real> type7;
  inline const Grid<Real> &getArg8()
  {
    return potWC;
  }
  typedef Grid<Real> type8;
  inline const Grid<Real> &getArg9()
  {
    return potKE;
  }
  typedef Grid<Real> type9;
  inline const Grid<Real> &getArg10()
  {
    return neighborRatio;
  }
  typedef Grid<Real> type10;
  inline const Real &getArg11()
  {
    return c_s;
  }
  typedef Real type11;
  inline const Real &getArg12()
  {
    return c_b;
  }
  typedef Real type12;
  inline const Real &getArg13()
  {
    return k_ta;
  }
  typedef Real type13;
  inline const Real &getArg14()
  {
    return k_wc;
  }
  typedef Real type14;
  inline const Real &getArg15()
  {
    return dt;
  }
  typedef Real type15;
  inline const int &getArg16()
  {
    return itype;
  }
  typedef int type16;
  void runMessage()
  {
    debMsg("Executing kernel knFlipSampleSecondaryParticles ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void run()
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    for (int k = minZ; k < maxZ; k++)
      for (int j = 0; j < _maxY; j++)
        for (int i = 0; i < _maxX; i++)
          op(i,
             j,
             k,
             flags,
             v,
             pts_sec,
             v_sec,
             l_sec,
             lMin,
             lMax,
             potTA,
             potWC,
             potKE,
             neighborRatio,
             c_s,
             c_b,
             k_ta,
             k_wc,
             dt,
             itype);
  }
  const FlagGrid &flags;
  const MACGrid &v;
  BasicParticleSystem &pts_sec;
  ParticleDataImpl<Vec3> &v_sec;
  ParticleDataImpl<Real> &l_sec;
  const Real lMin;
  const Real lMax;
  const Grid<Real> &potTA;
  const Grid<Real> &potWC;
  const Grid<Real> &potKE;
  const Grid<Real> &neighborRatio;
  const Real c_s;
  const Real c_b;
  const Real k_ta;
  const Real k_wc;
  const Real dt;
  const int itype;
};

void flipSampleSecondaryParticles(const std::string mode,
                                  const FlagGrid &flags,
                                  const MACGrid &v,
                                  BasicParticleSystem &pts_sec,
                                  ParticleDataImpl<Vec3> &v_sec,
                                  ParticleDataImpl<Real> &l_sec,
                                  const Real lMin,
                                  const Real lMax,
                                  const Grid<Real> &potTA,
                                  const Grid<Real> &potWC,
                                  const Grid<Real> &potKE,
                                  const Grid<Real> &neighborRatio,
                                  const Real c_s,
                                  const Real c_b,
                                  const Real k_ta,
                                  const Real k_wc,
                                  const Real dt,
                                  const int itype = FlagGrid::TypeFluid)
{
  if (mode == "single") {
    knFlipSampleSecondaryParticles(flags,
                                   v,
                                   pts_sec,
                                   v_sec,
                                   l_sec,
                                   lMin,
                                   lMax,
                                   potTA,
                                   potWC,
                                   potKE,
                                   neighborRatio,
                                   c_s,
                                   c_b,
                                   k_ta,
                                   k_wc,
                                   dt,
                                   itype);
  }
  else if (mode == "multiple") {
    knFlipSampleSecondaryParticlesMoreCylinders(flags,
                                                v,
                                                pts_sec,
                                                v_sec,
                                                l_sec,
                                                lMin,
                                                lMax,
                                                potTA,
                                                potWC,
                                                potKE,
                                                neighborRatio,
                                                c_s,
                                                c_b,
                                                k_ta,
                                                k_wc,
                                                dt,
                                                itype);
  }
  else {
    throw std::invalid_argument("Unknown mode: use \"single\" or \"multiple\" instead!");
  }
}
static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipSampleSecondaryParticles", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const std::string mode = _args.get<std::string>("mode", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      const MACGrid &v = *_args.getPtr<MACGrid>("v", 2, &_lock);
      BasicParticleSystem &pts_sec = *_args.getPtr<BasicParticleSystem>("pts_sec", 3, &_lock);
      ParticleDataImpl<Vec3> &v_sec = *_args.getPtr<ParticleDataImpl<Vec3>>("v_sec", 4, &_lock);
      ParticleDataImpl<Real> &l_sec = *_args.getPtr<ParticleDataImpl<Real>>("l_sec", 5, &_lock);
      const Real lMin = _args.get<Real>("lMin", 6, &_lock);
      const Real lMax = _args.get<Real>("lMax", 7, &_lock);
      const Grid<Real> &potTA = *_args.getPtr<Grid<Real>>("potTA", 8, &_lock);
      const Grid<Real> &potWC = *_args.getPtr<Grid<Real>>("potWC", 9, &_lock);
      const Grid<Real> &potKE = *_args.getPtr<Grid<Real>>("potKE", 10, &_lock);
      const Grid<Real> &neighborRatio = *_args.getPtr<Grid<Real>>("neighborRatio", 11, &_lock);
      const Real c_s = _args.get<Real>("c_s", 12, &_lock);
      const Real c_b = _args.get<Real>("c_b", 13, &_lock);
      const Real k_ta = _args.get<Real>("k_ta", 14, &_lock);
      const Real k_wc = _args.get<Real>("k_wc", 15, &_lock);
      const Real dt = _args.get<Real>("dt", 16, &_lock);
      const int itype = _args.getOpt<int>("itype", 17, FlagGrid::TypeFluid, &_lock);
      _retval = getPyNone();
      flipSampleSecondaryParticles(mode,
                                   flags,
                                   v,
                                   pts_sec,
                                   v_sec,
                                   l_sec,
                                   lMin,
                                   lMax,
                                   potTA,
                                   potWC,
                                   potKE,
                                   neighborRatio,
                                   c_s,
                                   c_b,
                                   k_ta,
                                   k_wc,
                                   dt,
                                   itype);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipSampleSecondaryParticles", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipSampleSecondaryParticles", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipSampleSecondaryParticles("",
                                                           "flipSampleSecondaryParticles",
                                                           _W_1);
extern "C" {
void PbRegister_flipSampleSecondaryParticles()
{
  KEEP_UNUSED(_RP_flipSampleSecondaryParticles);
}
}

// evaluates cubic spline with radius h and distance l in dim dimensions
Real cubicSpline(const Real h, const Real l, const int dim)
{
  const Real h2 = square(h), h3 = h2 * h;
  const Real c[] = {
      Real(2e0 / (3e0 * h)), Real(10e0 / (7e0 * M_PI * h2)), Real(1e0 / (M_PI * h3))};
  const Real q = l / h;
  if (q < 1e0)
    return c[dim - 1] * (1e0 - 1.5 * square(q) + 0.75 * cubed(q));
  else if (q < 2e0)
    return c[dim - 1] * (0.25 * cubed(2e0 - q));
  return 0;
}

// updates position &pts_sec.pos and velocity &v_sec of secondary particles according to the
// particle type determined by the neighbor ratio with linear interpolation

struct knFlipUpdateSecondaryParticlesLinear : public KernelBase {
  knFlipUpdateSecondaryParticlesLinear(BasicParticleSystem &pts_sec,
                                       ParticleDataImpl<Vec3> &v_sec,
                                       ParticleDataImpl<Real> &l_sec,
                                       const ParticleDataImpl<Vec3> &f_sec,
                                       const FlagGrid &flags,
                                       const MACGrid &v,
                                       const Grid<Real> &neighborRatio,
                                       const Vec3 gravity,
                                       const Real k_b,
                                       const Real k_d,
                                       const Real c_s,
                                       const Real c_b,
                                       const Real dt,
                                       const int exclude,
                                       const int antitunneling)
      : KernelBase(pts_sec.size()),
        pts_sec(pts_sec),
        v_sec(v_sec),
        l_sec(l_sec),
        f_sec(f_sec),
        flags(flags),
        v(v),
        neighborRatio(neighborRatio),
        gravity(gravity),
        k_b(k_b),
        k_d(k_d),
        c_s(c_s),
        c_b(c_b),
        dt(dt),
        exclude(exclude),
        antitunneling(antitunneling)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 BasicParticleSystem &pts_sec,
                 ParticleDataImpl<Vec3> &v_sec,
                 ParticleDataImpl<Real> &l_sec,
                 const ParticleDataImpl<Vec3> &f_sec,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 const Grid<Real> &neighborRatio,
                 const Vec3 gravity,
                 const Real k_b,
                 const Real k_d,
                 const Real c_s,
                 const Real c_b,
                 const Real dt,
                 const int exclude,
                 const int antitunneling) const
  {

    if (!pts_sec.isActive(idx) || pts_sec[idx].flag & exclude)
      return;
    if (!flags.isInBounds(pts_sec[idx].pos)) {
      pts_sec.kill(idx);
      return;
    }

    Vec3i gridpos = toVec3i(pts_sec[idx].pos);

    // spray particle
    if (neighborRatio(gridpos) < c_s) {
      pts_sec[idx].flag |= ParticleBase::PSPRAY;
      pts_sec[idx].flag &= ~(ParticleBase::PBUBBLE | ParticleBase::PFOAM);
      v_sec[idx] += dt *
                    ((f_sec[idx] / 1) + gravity);  // TODO: if forces are added (e.g. fluid
                                                   // guiding), add parameter for mass instead of 1

      // anti tunneling for small obstacles
      for (int ct = 1; ct < antitunneling; ct++) {
        Vec3i tempPos = toVec3i(pts_sec[idx].pos +
                                ct * (1 / Real(antitunneling)) * dt * v_sec[idx]);
        if (!flags.isInBounds(tempPos) || flags(tempPos) & FlagGrid::TypeObstacle) {
          pts_sec.kill(idx);
          return;
        }
      }
      pts_sec[idx].pos += dt * v_sec[idx];
    }

    // air bubble particle
    else if (neighborRatio(gridpos) > c_b) {
      pts_sec[idx].flag |= ParticleBase::PBUBBLE;
      pts_sec[idx].flag &= ~(ParticleBase::PSPRAY | ParticleBase::PFOAM);

      const Vec3 vj = (v.getInterpolated(pts_sec[idx].pos) - v_sec[idx]) / dt;
      v_sec[idx] += dt * (k_b * -gravity + k_d * vj);

      // anti tunneling for small obstacles
      for (int ct = 1; ct < antitunneling; ct++) {
        Vec3i tempPos = toVec3i(pts_sec[idx].pos +
                                ct * (1 / Real(antitunneling)) * dt * v_sec[idx]);
        if (!flags.isInBounds(tempPos) || flags(tempPos) & FlagGrid::TypeObstacle) {
          pts_sec.kill(idx);
          return;
        }
      }
      pts_sec[idx].pos += dt * v_sec[idx];
    }

    // foam particle
    else {
      pts_sec[idx].flag |= ParticleBase::PFOAM;
      pts_sec[idx].flag &= ~(ParticleBase::PBUBBLE | ParticleBase::PSPRAY);

      const Vec3 vj = v.getInterpolated(pts_sec[idx].pos);
      // anti tunneling for small obstacles
      for (int ct = 1; ct < antitunneling; ct++) {
        Vec3i tempPos = toVec3i(pts_sec[idx].pos + ct * (1 / Real(antitunneling)) * dt * vj);
        if (!flags.isInBounds(tempPos) || flags(tempPos) & FlagGrid::TypeObstacle) {
          pts_sec.kill(idx);
          return;
        }
      }
      pts_sec[idx].pos += dt * v.getInterpolated(pts_sec[idx].pos);
    }

    // lifetime
    l_sec[idx] -= dt;
    if (l_sec[idx] <= Real(0)) {
      pts_sec.kill(idx);
    }
  }
  inline BasicParticleSystem &getArg0()
  {
    return pts_sec;
  }
  typedef BasicParticleSystem type0;
  inline ParticleDataImpl<Vec3> &getArg1()
  {
    return v_sec;
  }
  typedef ParticleDataImpl<Vec3> type1;
  inline ParticleDataImpl<Real> &getArg2()
  {
    return l_sec;
  }
  typedef ParticleDataImpl<Real> type2;
  inline const ParticleDataImpl<Vec3> &getArg3()
  {
    return f_sec;
  }
  typedef ParticleDataImpl<Vec3> type3;
  inline const FlagGrid &getArg4()
  {
    return flags;
  }
  typedef FlagGrid type4;
  inline const MACGrid &getArg5()
  {
    return v;
  }
  typedef MACGrid type5;
  inline const Grid<Real> &getArg6()
  {
    return neighborRatio;
  }
  typedef Grid<Real> type6;
  inline const Vec3 &getArg7()
  {
    return gravity;
  }
  typedef Vec3 type7;
  inline const Real &getArg8()
  {
    return k_b;
  }
  typedef Real type8;
  inline const Real &getArg9()
  {
    return k_d;
  }
  typedef Real type9;
  inline const Real &getArg10()
  {
    return c_s;
  }
  typedef Real type10;
  inline const Real &getArg11()
  {
    return c_b;
  }
  typedef Real type11;
  inline const Real &getArg12()
  {
    return dt;
  }
  typedef Real type12;
  inline const int &getArg13()
  {
    return exclude;
  }
  typedef int type13;
  inline const int &getArg14()
  {
    return antitunneling;
  }
  typedef int type14;
  void runMessage()
  {
    debMsg("Executing kernel knFlipUpdateSecondaryParticlesLinear ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx,
         pts_sec,
         v_sec,
         l_sec,
         f_sec,
         flags,
         v,
         neighborRatio,
         gravity,
         k_b,
         k_d,
         c_s,
         c_b,
         dt,
         exclude,
         antitunneling);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  BasicParticleSystem &pts_sec;
  ParticleDataImpl<Vec3> &v_sec;
  ParticleDataImpl<Real> &l_sec;
  const ParticleDataImpl<Vec3> &f_sec;
  const FlagGrid &flags;
  const MACGrid &v;
  const Grid<Real> &neighborRatio;
  const Vec3 gravity;
  const Real k_b;
  const Real k_d;
  const Real c_s;
  const Real c_b;
  const Real dt;
  const int exclude;
  const int antitunneling;
};
// updates position &pts_sec.pos and velocity &v_sec of secondary particles according to the
// particle type determined by the neighbor ratio with cubic spline interpolation

struct knFlipUpdateSecondaryParticlesCubic : public KernelBase {
  knFlipUpdateSecondaryParticlesCubic(BasicParticleSystem &pts_sec,
                                      ParticleDataImpl<Vec3> &v_sec,
                                      ParticleDataImpl<Real> &l_sec,
                                      const ParticleDataImpl<Vec3> &f_sec,
                                      const FlagGrid &flags,
                                      const MACGrid &v,
                                      const Grid<Real> &neighborRatio,
                                      const int radius,
                                      const Vec3 gravity,
                                      const Real k_b,
                                      const Real k_d,
                                      const Real c_s,
                                      const Real c_b,
                                      const Real dt,
                                      const int exclude,
                                      const int antitunneling,
                                      const int itype)
      : KernelBase(pts_sec.size()),
        pts_sec(pts_sec),
        v_sec(v_sec),
        l_sec(l_sec),
        f_sec(f_sec),
        flags(flags),
        v(v),
        neighborRatio(neighborRatio),
        radius(radius),
        gravity(gravity),
        k_b(k_b),
        k_d(k_d),
        c_s(c_s),
        c_b(c_b),
        dt(dt),
        exclude(exclude),
        antitunneling(antitunneling),
        itype(itype)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 BasicParticleSystem &pts_sec,
                 ParticleDataImpl<Vec3> &v_sec,
                 ParticleDataImpl<Real> &l_sec,
                 const ParticleDataImpl<Vec3> &f_sec,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 const Grid<Real> &neighborRatio,
                 const int radius,
                 const Vec3 gravity,
                 const Real k_b,
                 const Real k_d,
                 const Real c_s,
                 const Real c_b,
                 const Real dt,
                 const int exclude,
                 const int antitunneling,
                 const int itype) const
  {

    if (!pts_sec.isActive(idx) || pts_sec[idx].flag & exclude)
      return;
    if (!flags.isInBounds(pts_sec[idx].pos)) {
      pts_sec.kill(idx);
      return;
    }

    Vec3i gridpos = toVec3i(pts_sec[idx].pos);
    int i = gridpos.x;
    int j = gridpos.y;
    int k = gridpos.z;

    // spray particle
    if (neighborRatio(gridpos) < c_s) {
      pts_sec[idx].flag |= ParticleBase::PSPRAY;
      pts_sec[idx].flag &= ~(ParticleBase::PBUBBLE | ParticleBase::PFOAM);
      v_sec[idx] += dt *
                    ((f_sec[idx] / 1) + gravity);  // TODO: if forces are added (e.g. fluid
                                                   // guiding), add parameter for mass instead of 1

      // anti tunneling for small obstacles
      for (int ct = 1; ct < antitunneling; ct++) {
        Vec3i tempPos = toVec3i(pts_sec[idx].pos +
                                ct * (1 / Real(antitunneling)) * dt * v_sec[idx]);
        if (!flags.isInBounds(tempPos) || flags(tempPos) & FlagGrid::TypeObstacle) {
          pts_sec.kill(idx);
          return;
        }
      }
      pts_sec[idx].pos += dt * v_sec[idx];
    }

    // air bubble particle
    else if (neighborRatio(gridpos) > c_b) {
      pts_sec[idx].flag |= ParticleBase::PBUBBLE;
      pts_sec[idx].flag &= ~(ParticleBase::PSPRAY | ParticleBase::PFOAM);
      const Vec3 &xi = pts_sec[idx].pos;
      Vec3 sumNumerator = Vec3(0, 0, 0);
      Real sumDenominator = 0;
      for (IndexInt x = i - radius; x <= i + radius; x++) {
        for (IndexInt y = j - radius; y <= j + radius; y++) {
          for (IndexInt z = k - radius; z <= k + radius; z++) {
            Vec3i xj = Vec3i(x, y, z);
            if ((x == i && y == j && z == k) || !flags.isInBounds(xj))
              continue;
            if (!(flags(xj) & itype))
              continue;
            const Real len_xij = norm(xi - Vec3(x, y, z));

            int dim = flags.is3D() ? 3 : 2;
            Real dist = flags.is3D() ? 1.732 : 1.414;
            Real weight = cubicSpline(radius * dist, len_xij, dim);
            sumNumerator += v.getCentered(xj) *
                            weight;  // estimate next position by current velocity
            sumDenominator += weight;
          }
        }
      }
      const Vec3 temp = ((sumNumerator / sumDenominator) - v_sec[idx]) / dt;
      v_sec[idx] += dt * (k_b * -gravity + k_d * temp);

      // anti tunneling for small obstacles
      for (int ct = 1; ct < antitunneling; ct++) {
        Vec3i tempPos = toVec3i(pts_sec[idx].pos +
                                ct * (1 / Real(antitunneling)) * dt * v_sec[idx]);
        if (!flags.isInBounds(tempPos) || flags(tempPos) & FlagGrid::TypeObstacle) {
          pts_sec.kill(idx);
          return;
        }
      }
      pts_sec[idx].pos += dt * v_sec[idx];
    }

    // foam particle
    else {
      pts_sec[idx].flag |= ParticleBase::PFOAM;
      pts_sec[idx].flag &= ~(ParticleBase::PBUBBLE | ParticleBase::PSPRAY);
      const Vec3 &xi = pts_sec[idx].pos;
      Vec3 sumNumerator = Vec3(0, 0, 0);
      Real sumDenominator = 0;
      for (IndexInt x = i - radius; x <= i + radius; x++) {
        for (IndexInt y = j - radius; y <= j + radius; y++) {
          for (IndexInt z = k - radius; z <= k + radius; z++) {
            Vec3i xj = Vec3i(x, y, z);
            if ((x == i && y == j && z == k) || !flags.isInBounds(xj))
              continue;
            if (!(flags(xj) & itype))
              continue;
            const Real len_xij = norm(xi - Vec3(x, y, z));

            int dim = flags.is3D() ? 3 : 2;
            Real dist = flags.is3D() ? 1.732 : 1.414;
            Real weight = cubicSpline(radius * dist, len_xij, dim);
            sumNumerator += v.getCentered(xj) *
                            weight;  // estimate next position by current velocity
            sumDenominator += weight;
          }
        }
      }

      // anti tunneling for small obstacles
      for (int ct = 1; ct < antitunneling; ct++) {
        Vec3i tempPos = toVec3i(pts_sec[idx].pos + ct * (1 / Real(antitunneling)) * dt *
                                                       (sumNumerator / sumDenominator));
        if (!flags.isInBounds(tempPos) || flags(tempPos) & FlagGrid::TypeObstacle) {
          pts_sec.kill(idx);
          return;
        }
      }
      pts_sec[idx].pos += dt * (sumNumerator / sumDenominator);
    }

    // lifetime
    l_sec[idx] -= dt;
    if (l_sec[idx] <= Real(0)) {
      pts_sec.kill(idx);
    }
  }
  inline BasicParticleSystem &getArg0()
  {
    return pts_sec;
  }
  typedef BasicParticleSystem type0;
  inline ParticleDataImpl<Vec3> &getArg1()
  {
    return v_sec;
  }
  typedef ParticleDataImpl<Vec3> type1;
  inline ParticleDataImpl<Real> &getArg2()
  {
    return l_sec;
  }
  typedef ParticleDataImpl<Real> type2;
  inline const ParticleDataImpl<Vec3> &getArg3()
  {
    return f_sec;
  }
  typedef ParticleDataImpl<Vec3> type3;
  inline const FlagGrid &getArg4()
  {
    return flags;
  }
  typedef FlagGrid type4;
  inline const MACGrid &getArg5()
  {
    return v;
  }
  typedef MACGrid type5;
  inline const Grid<Real> &getArg6()
  {
    return neighborRatio;
  }
  typedef Grid<Real> type6;
  inline const int &getArg7()
  {
    return radius;
  }
  typedef int type7;
  inline const Vec3 &getArg8()
  {
    return gravity;
  }
  typedef Vec3 type8;
  inline const Real &getArg9()
  {
    return k_b;
  }
  typedef Real type9;
  inline const Real &getArg10()
  {
    return k_d;
  }
  typedef Real type10;
  inline const Real &getArg11()
  {
    return c_s;
  }
  typedef Real type11;
  inline const Real &getArg12()
  {
    return c_b;
  }
  typedef Real type12;
  inline const Real &getArg13()
  {
    return dt;
  }
  typedef Real type13;
  inline const int &getArg14()
  {
    return exclude;
  }
  typedef int type14;
  inline const int &getArg15()
  {
    return antitunneling;
  }
  typedef int type15;
  inline const int &getArg16()
  {
    return itype;
  }
  typedef int type16;
  void runMessage()
  {
    debMsg("Executing kernel knFlipUpdateSecondaryParticlesCubic ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx,
         pts_sec,
         v_sec,
         l_sec,
         f_sec,
         flags,
         v,
         neighborRatio,
         radius,
         gravity,
         k_b,
         k_d,
         c_s,
         c_b,
         dt,
         exclude,
         antitunneling,
         itype);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  BasicParticleSystem &pts_sec;
  ParticleDataImpl<Vec3> &v_sec;
  ParticleDataImpl<Real> &l_sec;
  const ParticleDataImpl<Vec3> &f_sec;
  const FlagGrid &flags;
  const MACGrid &v;
  const Grid<Real> &neighborRatio;
  const int radius;
  const Vec3 gravity;
  const Real k_b;
  const Real k_d;
  const Real c_s;
  const Real c_b;
  const Real dt;
  const int exclude;
  const int antitunneling;
  const int itype;
};

void flipUpdateSecondaryParticles(const std::string mode,
                                  BasicParticleSystem &pts_sec,
                                  ParticleDataImpl<Vec3> &v_sec,
                                  ParticleDataImpl<Real> &l_sec,
                                  const ParticleDataImpl<Vec3> &f_sec,
                                  FlagGrid &flags,
                                  const MACGrid &v,
                                  const Grid<Real> &neighborRatio,
                                  const int radius,
                                  const Vec3 gravity,
                                  const Real k_b,
                                  const Real k_d,
                                  const Real c_s,
                                  const Real c_b,
                                  const Real dt,
                                  const int exclude = ParticleBase::PTRACER,
                                  const int antitunneling = 0,
                                  const int itype = FlagGrid::TypeFluid)
{

  Vec3 g = gravity / flags.getDx();
  if (mode == "linear") {
    knFlipUpdateSecondaryParticlesLinear(pts_sec,
                                         v_sec,
                                         l_sec,
                                         f_sec,
                                         flags,
                                         v,
                                         neighborRatio,
                                         g,
                                         k_b,
                                         k_d,
                                         c_s,
                                         c_b,
                                         dt,
                                         exclude,
                                         antitunneling);
  }
  else if (mode == "cubic") {
    knFlipUpdateSecondaryParticlesCubic(pts_sec,
                                        v_sec,
                                        l_sec,
                                        f_sec,
                                        flags,
                                        v,
                                        neighborRatio,
                                        radius,
                                        g,
                                        k_b,
                                        k_d,
                                        c_s,
                                        c_b,
                                        dt,
                                        exclude,
                                        antitunneling,
                                        itype);
  }
  else {
    throw std::invalid_argument("Unknown mode: use \"linear\" or \"cubic\" instead!");
  }
  pts_sec.doCompress();
}
static PyObject *_W_2(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipUpdateSecondaryParticles", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const std::string mode = _args.get<std::string>("mode", 0, &_lock);
      BasicParticleSystem &pts_sec = *_args.getPtr<BasicParticleSystem>("pts_sec", 1, &_lock);
      ParticleDataImpl<Vec3> &v_sec = *_args.getPtr<ParticleDataImpl<Vec3>>("v_sec", 2, &_lock);
      ParticleDataImpl<Real> &l_sec = *_args.getPtr<ParticleDataImpl<Real>>("l_sec", 3, &_lock);
      const ParticleDataImpl<Vec3> &f_sec = *_args.getPtr<ParticleDataImpl<Vec3>>(
          "f_sec", 4, &_lock);
      FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 5, &_lock);
      const MACGrid &v = *_args.getPtr<MACGrid>("v", 6, &_lock);
      const Grid<Real> &neighborRatio = *_args.getPtr<Grid<Real>>("neighborRatio", 7, &_lock);
      const int radius = _args.get<int>("radius", 8, &_lock);
      const Vec3 gravity = _args.get<Vec3>("gravity", 9, &_lock);
      const Real k_b = _args.get<Real>("k_b", 10, &_lock);
      const Real k_d = _args.get<Real>("k_d", 11, &_lock);
      const Real c_s = _args.get<Real>("c_s", 12, &_lock);
      const Real c_b = _args.get<Real>("c_b", 13, &_lock);
      const Real dt = _args.get<Real>("dt", 14, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 15, ParticleBase::PTRACER, &_lock);
      const int antitunneling = _args.getOpt<int>("antitunneling", 16, 0, &_lock);
      const int itype = _args.getOpt<int>("itype", 17, FlagGrid::TypeFluid, &_lock);
      _retval = getPyNone();
      flipUpdateSecondaryParticles(mode,
                                   pts_sec,
                                   v_sec,
                                   l_sec,
                                   f_sec,
                                   flags,
                                   v,
                                   neighborRatio,
                                   radius,
                                   gravity,
                                   k_b,
                                   k_d,
                                   c_s,
                                   c_b,
                                   dt,
                                   exclude,
                                   antitunneling,
                                   itype);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipUpdateSecondaryParticles", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipUpdateSecondaryParticles", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipUpdateSecondaryParticles("",
                                                           "flipUpdateSecondaryParticles",
                                                           _W_2);
extern "C" {
void PbRegister_flipUpdateSecondaryParticles()
{
  KEEP_UNUSED(_RP_flipUpdateSecondaryParticles);
}
}

// removes secondary particles in &pts_sec that are inside boundaries (cells that are marked as
// obstacle/outflow in &flags)

struct knFlipDeleteParticlesInObstacle : public KernelBase {
  knFlipDeleteParticlesInObstacle(BasicParticleSystem &pts, const FlagGrid &flags)
      : KernelBase(pts.size()), pts(pts), flags(flags)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, BasicParticleSystem &pts, const FlagGrid &flags) const
  {

    if (!pts.isActive(idx))
      return;

    const Vec3 &xi = pts[idx].pos;
    const Vec3i xidx = toVec3i(xi);
    // remove particles that completely left the bounds
    if (!flags.isInBounds(xidx)) {
      pts.kill(idx);
      return;
    }
    int gridIndex = flags.index(xidx);
    // remove particles that penetrate obstacles
    if (flags[gridIndex] == FlagGrid::TypeObstacle || flags[gridIndex] == FlagGrid::TypeOutflow) {
      pts.kill(idx);
    }
  }
  inline BasicParticleSystem &getArg0()
  {
    return pts;
  }
  typedef BasicParticleSystem type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  void runMessage()
  {
    debMsg("Executing kernel knFlipDeleteParticlesInObstacle ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, pts, flags);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  BasicParticleSystem &pts;
  const FlagGrid &flags;
};

void flipDeleteParticlesInObstacle(BasicParticleSystem &pts, const FlagGrid &flags)
{

  knFlipDeleteParticlesInObstacle(pts, flags);
  pts.doCompress();
}
static PyObject *_W_3(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipDeleteParticlesInObstacle", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      BasicParticleSystem &pts = *_args.getPtr<BasicParticleSystem>("pts", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      _retval = getPyNone();
      flipDeleteParticlesInObstacle(pts, flags);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipDeleteParticlesInObstacle", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipDeleteParticlesInObstacle", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipDeleteParticlesInObstacle("",
                                                            "flipDeleteParticlesInObstacle",
                                                            _W_3);
extern "C" {
void PbRegister_flipDeleteParticlesInObstacle()
{
  KEEP_UNUSED(_RP_flipDeleteParticlesInObstacle);
}
}

// helper method to debug statistical data from grid

void debugGridInfo(const FlagGrid &flags,
                   Grid<Real> &grid,
                   std::string name,
                   const int itype = FlagGrid::TypeFluid)
{
  FluidSolver *s = flags.getParent();
  int countFluid = 0;
  int countLargerZero = 0;
  Real avg = 0;
  Real max = 0;
  Real sum = 0;
  Real avgLargerZero = 0;
  FOR_IJK_BND(grid, 1)
  {
    if (!(flags(i, j, k) & itype))
      continue;
    countFluid++;
    if (grid(i, j, k) > 0)
      countLargerZero++;
    sum += grid(i, j, k);
    if (grid(i, j, k) > max)
      max = grid(i, j, k);
  }
  avg = sum / std::max(Real(countFluid), Real(1));
  avgLargerZero = sum / std::max(Real(countLargerZero), Real(1));

  debMsg("Step: " << s->mFrame << " - Grid " << name << "\n\tcountFluid \t\t" << countFluid
                  << "\n\tcountLargerZero \t" << countLargerZero << "\n\tsum \t\t\t" << sum
                  << "\n\tavg \t\t\t" << avg << "\n\tavgLargerZero \t\t" << avgLargerZero
                  << "\n\tmax \t\t\t" << max,
         1);
}
static PyObject *_W_4(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "debugGridInfo", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      Grid<Real> &grid = *_args.getPtr<Grid<Real>>("grid", 1, &_lock);
      std::string name = _args.get<std::string>("name", 2, &_lock);
      const int itype = _args.getOpt<int>("itype", 3, FlagGrid::TypeFluid, &_lock);
      _retval = getPyNone();
      debugGridInfo(flags, grid, name, itype);
      _args.check();
    }
    pbFinalizePlugin(parent, "debugGridInfo", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("debugGridInfo", e.what());
    return 0;
  }
}
static const Pb::Register _RP_debugGridInfo("", "debugGridInfo", _W_4);
extern "C" {
void PbRegister_debugGridInfo()
{
  KEEP_UNUSED(_RP_debugGridInfo);
}
}

// The following methods are helper functions to recreate the velocity and flag grid from the
// underlying FLIP simulation. They cannot simply be loaded because of the upres to a higher
// resolution, instead a levelset is used.

struct knSetFlagsFromLevelset : public KernelBase {
  knSetFlagsFromLevelset(FlagGrid &flags,
                         const Grid<Real> &phi,
                         const int exclude = FlagGrid::TypeObstacle,
                         const int itype = FlagGrid::TypeFluid)
      : KernelBase(&flags, 0), flags(flags), phi(phi), exclude(exclude), itype(itype)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 FlagGrid &flags,
                 const Grid<Real> &phi,
                 const int exclude = FlagGrid::TypeObstacle,
                 const int itype = FlagGrid::TypeFluid) const
  {
    if (phi(idx) < 0 && !(flags(idx) & exclude))
      flags(idx) = itype;
  }
  inline FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const Grid<Real> &getArg1()
  {
    return phi;
  }
  typedef Grid<Real> type1;
  inline const int &getArg2()
  {
    return exclude;
  }
  typedef int type2;
  inline const int &getArg3()
  {
    return itype;
  }
  typedef int type3;
  void runMessage()
  {
    debMsg("Executing kernel knSetFlagsFromLevelset ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, flags, phi, exclude, itype);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  FlagGrid &flags;
  const Grid<Real> &phi;
  const int exclude;
  const int itype;
};

void setFlagsFromLevelset(FlagGrid &flags,
                          const Grid<Real> &phi,
                          const int exclude = FlagGrid::TypeObstacle,
                          const int itype = FlagGrid::TypeFluid)
{
  knSetFlagsFromLevelset(flags, phi, exclude, itype);
}
static PyObject *_W_5(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "setFlagsFromLevelset", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      const Grid<Real> &phi = *_args.getPtr<Grid<Real>>("phi", 1, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 2, FlagGrid::TypeObstacle, &_lock);
      const int itype = _args.getOpt<int>("itype", 3, FlagGrid::TypeFluid, &_lock);
      _retval = getPyNone();
      setFlagsFromLevelset(flags, phi, exclude, itype);
      _args.check();
    }
    pbFinalizePlugin(parent, "setFlagsFromLevelset", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("setFlagsFromLevelset", e.what());
    return 0;
  }
}
static const Pb::Register _RP_setFlagsFromLevelset("", "setFlagsFromLevelset", _W_5);
extern "C" {
void PbRegister_setFlagsFromLevelset()
{
  KEEP_UNUSED(_RP_setFlagsFromLevelset);
}
}

struct knSetMACFromLevelset : public KernelBase {
  knSetMACFromLevelset(MACGrid &v, const Grid<Real> &phi, const Vec3 c)
      : KernelBase(&v, 0), v(v), phi(phi), c(c)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, MACGrid &v, const Grid<Real> &phi, const Vec3 c) const
  {
    if (phi.getInterpolated(Vec3(i, j, k)) > 0)
      v(i, j, k) = c;
  }
  inline MACGrid &getArg0()
  {
    return v;
  }
  typedef MACGrid type0;
  inline const Grid<Real> &getArg1()
  {
    return phi;
  }
  typedef Grid<Real> type1;
  inline const Vec3 &getArg2()
  {
    return c;
  }
  typedef Vec3 type2;
  void runMessage()
  {
    debMsg("Executing kernel knSetMACFromLevelset ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 0; j < _maxY; j++)
          for (int i = 0; i < _maxX; i++)
            op(i, j, k, v, phi, c);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 0; i < _maxX; i++)
          op(i, j, k, v, phi, c);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  MACGrid &v;
  const Grid<Real> &phi;
  const Vec3 c;
};

void setMACFromLevelset(MACGrid &v, const Grid<Real> &phi, const Vec3 c)
{
  knSetMACFromLevelset(v, phi, c);
}
static PyObject *_W_6(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "setMACFromLevelset", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      MACGrid &v = *_args.getPtr<MACGrid>("v", 0, &_lock);
      const Grid<Real> &phi = *_args.getPtr<Grid<Real>>("phi", 1, &_lock);
      const Vec3 c = _args.get<Vec3>("c", 2, &_lock);
      _retval = getPyNone();
      setMACFromLevelset(v, phi, c);
      _args.check();
    }
    pbFinalizePlugin(parent, "setMACFromLevelset", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("setMACFromLevelset", e.what());
    return 0;
  }
}
static const Pb::Register _RP_setMACFromLevelset("", "setMACFromLevelset", _W_6);
extern "C" {
void PbRegister_setMACFromLevelset()
{
  KEEP_UNUSED(_RP_setMACFromLevelset);
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------
// END Secondary Particles for FLIP
//----------------------------------------------------------------------------------------------------------------------------------------------------
#pragma endregion

#pragma region Legacy Methods(still useful for debugging)
//-----------------------------------------------------------------------------------------------------------------------------------
//-----------------
// Legacy Methods (still useful for debugging)
//----------------------------------------------------------------------------------------------------------------------------------------------------

// LEGACY METHOD! Use flipComputeSecondaryParticlePotentials instead!
// computes trapped air potential for all fluid cells in &flags and saves it in &pot

struct knFlipComputePotentialTrappedAir : public KernelBase {
  knFlipComputePotentialTrappedAir(Grid<Real> &pot,
                                   const FlagGrid &flags,
                                   const MACGrid &v,
                                   const int radius,
                                   const Real tauMin,
                                   const Real tauMax,
                                   const Real scaleFromManta,
                                   const int itype = FlagGrid::TypeFluid,
                                   const int jtype = FlagGrid::TypeFluid)
      : KernelBase(&pot, 1),
        pot(pot),
        flags(flags),
        v(v),
        radius(radius),
        tauMin(tauMin),
        tauMax(tauMax),
        scaleFromManta(scaleFromManta),
        itype(itype),
        jtype(jtype)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 Grid<Real> &pot,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 const int radius,
                 const Real tauMin,
                 const Real tauMax,
                 const Real scaleFromManta,
                 const int itype = FlagGrid::TypeFluid,
                 const int jtype = FlagGrid::TypeFluid) const
  {

    if (!(flags(i, j, k) & itype))
      return;

    const Vec3 &xi = scaleFromManta * Vec3(i, j, k);  // scale to unit cube
    const Vec3 &vi = scaleFromManta * v.getCentered(i, j, k);
    Real vdiff = 0;
    for (IndexInt x = i - radius; x <= i + radius; x++) {
      for (IndexInt y = j - radius; y <= j + radius; y++) {
        for (IndexInt z = k - radius; z <= k + radius; z++) {
          if ((x == i && y == j && z == k) || !(flags(x, y, z) & jtype))
            continue;

          const Vec3 &xj = scaleFromManta * Vec3(x, y, z);  // scale to unit cube
          const Vec3 &vj = scaleFromManta * v.getCentered(x, y, z);
          const Vec3 xij = xi - xj;
          const Vec3 vij = vi - vj;
          Real h = !pot.is3D() ? 1.414 * radius :
                                 1.732 * radius;  // estimate sqrt(2)*radius resp. sqrt(3)*radius
                                                  // for h, due to squared resp. cubic neighbor area
          vdiff += norm(vij) * (1 - dot(getNormalized(vij), getNormalized(xij))) *
                   (1 - norm(xij) / h);
        }
      }
    }
    pot(i, j, k) = (std::min(vdiff, tauMax) - std::min(vdiff, tauMin)) / (tauMax - tauMin);
  }
  inline Grid<Real> &getArg0()
  {
    return pot;
  }
  typedef Grid<Real> type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const MACGrid &getArg2()
  {
    return v;
  }
  typedef MACGrid type2;
  inline const int &getArg3()
  {
    return radius;
  }
  typedef int type3;
  inline const Real &getArg4()
  {
    return tauMin;
  }
  typedef Real type4;
  inline const Real &getArg5()
  {
    return tauMax;
  }
  typedef Real type5;
  inline const Real &getArg6()
  {
    return scaleFromManta;
  }
  typedef Real type6;
  inline const int &getArg7()
  {
    return itype;
  }
  typedef int type7;
  inline const int &getArg8()
  {
    return jtype;
  }
  typedef int type8;
  void runMessage()
  {
    debMsg("Executing kernel knFlipComputePotentialTrappedAir ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, pot, flags, v, radius, tauMin, tauMax, scaleFromManta, itype, jtype);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 1; i < _maxX; i++)
          op(i, j, k, pot, flags, v, radius, tauMin, tauMax, scaleFromManta, itype, jtype);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
  }
  Grid<Real> &pot;
  const FlagGrid &flags;
  const MACGrid &v;
  const int radius;
  const Real tauMin;
  const Real tauMax;
  const Real scaleFromManta;
  const int itype;
  const int jtype;
};

void flipComputePotentialTrappedAir(Grid<Real> &pot,
                                    const FlagGrid &flags,
                                    const MACGrid &v,
                                    const int radius,
                                    const Real tauMin,
                                    const Real tauMax,
                                    const Real scaleFromManta,
                                    const int itype = FlagGrid::TypeFluid,
                                    const int jtype = FlagGrid::TypeFluid)
{
  pot.clear();
  knFlipComputePotentialTrappedAir(
      pot, flags, v, radius, tauMin, tauMax, scaleFromManta, itype, jtype);
}
static PyObject *_W_7(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipComputePotentialTrappedAir", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &pot = *_args.getPtr<Grid<Real>>("pot", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      const MACGrid &v = *_args.getPtr<MACGrid>("v", 2, &_lock);
      const int radius = _args.get<int>("radius", 3, &_lock);
      const Real tauMin = _args.get<Real>("tauMin", 4, &_lock);
      const Real tauMax = _args.get<Real>("tauMax", 5, &_lock);
      const Real scaleFromManta = _args.get<Real>("scaleFromManta", 6, &_lock);
      const int itype = _args.getOpt<int>("itype", 7, FlagGrid::TypeFluid, &_lock);
      const int jtype = _args.getOpt<int>("jtype", 8, FlagGrid::TypeFluid, &_lock);
      _retval = getPyNone();
      flipComputePotentialTrappedAir(
          pot, flags, v, radius, tauMin, tauMax, scaleFromManta, itype, jtype);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipComputePotentialTrappedAir", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipComputePotentialTrappedAir", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipComputePotentialTrappedAir("",
                                                             "flipComputePotentialTrappedAir",
                                                             _W_7);
extern "C" {
void PbRegister_flipComputePotentialTrappedAir()
{
  KEEP_UNUSED(_RP_flipComputePotentialTrappedAir);
}
}

// LEGACY METHOD! Use flipComputeSecondaryParticlePotentials instead!
// computes kinetic energy potential for all fluid cells in &flags and saves it in &pot

struct knFlipComputePotentialKineticEnergy : public KernelBase {
  knFlipComputePotentialKineticEnergy(Grid<Real> &pot,
                                      const FlagGrid &flags,
                                      const MACGrid &v,
                                      const Real tauMin,
                                      const Real tauMax,
                                      const Real scaleFromManta,
                                      const int itype = FlagGrid::TypeFluid)
      : KernelBase(&pot, 0),
        pot(pot),
        flags(flags),
        v(v),
        tauMin(tauMin),
        tauMax(tauMax),
        scaleFromManta(scaleFromManta),
        itype(itype)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 Grid<Real> &pot,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 const Real tauMin,
                 const Real tauMax,
                 const Real scaleFromManta,
                 const int itype = FlagGrid::TypeFluid) const
  {

    if (!(flags(i, j, k) & itype))
      return;

    const Vec3 &vi = scaleFromManta * v.getCentered(i, j, k);  // scale to unit cube
    Real ek =
        Real(0.5) * 125 *
        normSquare(
            vi);  // use arbitrary constant for mass, potential adjusts with thresholds anyways
    pot(i, j, k) = (std::min(ek, tauMax) - std::min(ek, tauMin)) / (tauMax - tauMin);
  }
  inline Grid<Real> &getArg0()
  {
    return pot;
  }
  typedef Grid<Real> type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const MACGrid &getArg2()
  {
    return v;
  }
  typedef MACGrid type2;
  inline const Real &getArg3()
  {
    return tauMin;
  }
  typedef Real type3;
  inline const Real &getArg4()
  {
    return tauMax;
  }
  typedef Real type4;
  inline const Real &getArg5()
  {
    return scaleFromManta;
  }
  typedef Real type5;
  inline const int &getArg6()
  {
    return itype;
  }
  typedef int type6;
  void runMessage()
  {
    debMsg("Executing kernel knFlipComputePotentialKineticEnergy ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 0; j < _maxY; j++)
          for (int i = 0; i < _maxX; i++)
            op(i, j, k, pot, flags, v, tauMin, tauMax, scaleFromManta, itype);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 0; i < _maxX; i++)
          op(i, j, k, pot, flags, v, tauMin, tauMax, scaleFromManta, itype);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  Grid<Real> &pot;
  const FlagGrid &flags;
  const MACGrid &v;
  const Real tauMin;
  const Real tauMax;
  const Real scaleFromManta;
  const int itype;
};

void flipComputePotentialKineticEnergy(Grid<Real> &pot,
                                       const FlagGrid &flags,
                                       const MACGrid &v,
                                       const Real tauMin,
                                       const Real tauMax,
                                       const Real scaleFromManta,
                                       const int itype = FlagGrid::TypeFluid)
{
  pot.clear();
  knFlipComputePotentialKineticEnergy(pot, flags, v, tauMin, tauMax, scaleFromManta, itype);
}
static PyObject *_W_8(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipComputePotentialKineticEnergy", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &pot = *_args.getPtr<Grid<Real>>("pot", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      const MACGrid &v = *_args.getPtr<MACGrid>("v", 2, &_lock);
      const Real tauMin = _args.get<Real>("tauMin", 3, &_lock);
      const Real tauMax = _args.get<Real>("tauMax", 4, &_lock);
      const Real scaleFromManta = _args.get<Real>("scaleFromManta", 5, &_lock);
      const int itype = _args.getOpt<int>("itype", 6, FlagGrid::TypeFluid, &_lock);
      _retval = getPyNone();
      flipComputePotentialKineticEnergy(pot, flags, v, tauMin, tauMax, scaleFromManta, itype);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipComputePotentialKineticEnergy", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipComputePotentialKineticEnergy", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipComputePotentialKineticEnergy(
    "", "flipComputePotentialKineticEnergy", _W_8);
extern "C" {
void PbRegister_flipComputePotentialKineticEnergy()
{
  KEEP_UNUSED(_RP_flipComputePotentialKineticEnergy);
}
}

// LEGACY METHOD! Use flipComputeSecondaryParticlePotentials instead!
// computes wave crest potential for all fluid cells in &flags and saves it in &pot

struct knFlipComputePotentialWaveCrest : public KernelBase {
  knFlipComputePotentialWaveCrest(Grid<Real> &pot,
                                  const FlagGrid &flags,
                                  const MACGrid &v,
                                  const int radius,
                                  Grid<Vec3> &normal,
                                  const Real tauMin,
                                  const Real tauMax,
                                  const Real scaleFromManta,
                                  const int itype = FlagGrid::TypeFluid,
                                  const int jtype = FlagGrid::TypeFluid)
      : KernelBase(&pot, 1),
        pot(pot),
        flags(flags),
        v(v),
        radius(radius),
        normal(normal),
        tauMin(tauMin),
        tauMax(tauMax),
        scaleFromManta(scaleFromManta),
        itype(itype),
        jtype(jtype)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 Grid<Real> &pot,
                 const FlagGrid &flags,
                 const MACGrid &v,
                 const int radius,
                 Grid<Vec3> &normal,
                 const Real tauMin,
                 const Real tauMax,
                 const Real scaleFromManta,
                 const int itype = FlagGrid::TypeFluid,
                 const int jtype = FlagGrid::TypeFluid) const
  {

    if (!(flags(i, j, k) & itype))
      return;

    const Vec3 &xi = scaleFromManta * Vec3(i, j, k);  // scale to unit cube
    const Vec3 &vi = scaleFromManta * v.getCentered(i, j, k);
    const Vec3 &ni = normal(i, j, k);
    Real kappa = 0;
    for (IndexInt x = i - radius; x <= i + radius; x++) {
      for (IndexInt y = j - radius; y <= j + radius; y++) {
        for (IndexInt z = k - radius; z <= k + radius; z++) {
          if ((x == i && y == j && z == k) || !(flags(x, y, z) & jtype))
            continue;
          const Vec3 &xj = scaleFromManta * Vec3(x, y, z);  // scale to unit cube
          const Vec3 &nj = normal(x, y, z);
          const Vec3 xij = xi - xj;
          if (dot(getNormalized(xij), ni) < 0) {  // identifies wave crests
            Real h = !pot.is3D() ?
                         1.414 * radius :
                         1.732 * radius;  // estimate sqrt(2)*radius resp. sqrt(3)*radius for h,
                                          // due to squared resp. cubic neighbor area
            kappa += (1 - dot(ni, nj)) * (1 - norm(xij) / h);
          }
        }
      }
    }

    if (dot(getNormalized(vi), ni) >= 0.6) {  // avoid to mark boarders of the scene as wave crest
      pot(i, j, k) = (std::min(kappa, tauMax) - std::min(kappa, tauMin)) / (tauMax - tauMin);
    }
    else {
      pot(i, j, k) = Real(0);
    }
  }
  inline Grid<Real> &getArg0()
  {
    return pot;
  }
  typedef Grid<Real> type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const MACGrid &getArg2()
  {
    return v;
  }
  typedef MACGrid type2;
  inline const int &getArg3()
  {
    return radius;
  }
  typedef int type3;
  inline Grid<Vec3> &getArg4()
  {
    return normal;
  }
  typedef Grid<Vec3> type4;
  inline const Real &getArg5()
  {
    return tauMin;
  }
  typedef Real type5;
  inline const Real &getArg6()
  {
    return tauMax;
  }
  typedef Real type6;
  inline const Real &getArg7()
  {
    return scaleFromManta;
  }
  typedef Real type7;
  inline const int &getArg8()
  {
    return itype;
  }
  typedef int type8;
  inline const int &getArg9()
  {
    return jtype;
  }
  typedef int type9;
  void runMessage()
  {
    debMsg("Executing kernel knFlipComputePotentialWaveCrest ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i,
               j,
               k,
               pot,
               flags,
               v,
               radius,
               normal,
               tauMin,
               tauMax,
               scaleFromManta,
               itype,
               jtype);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 1; i < _maxX; i++)
          op(i, j, k, pot, flags, v, radius, normal, tauMin, tauMax, scaleFromManta, itype, jtype);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
  }
  Grid<Real> &pot;
  const FlagGrid &flags;
  const MACGrid &v;
  const int radius;
  Grid<Vec3> &normal;
  const Real tauMin;
  const Real tauMax;
  const Real scaleFromManta;
  const int itype;
  const int jtype;
};

void flipComputePotentialWaveCrest(Grid<Real> &pot,
                                   const FlagGrid &flags,
                                   const MACGrid &v,
                                   const int radius,
                                   Grid<Vec3> &normal,
                                   const Real tauMin,
                                   const Real tauMax,
                                   const Real scaleFromManta,
                                   const int itype = FlagGrid::TypeFluid,
                                   const int jtype = FlagGrid::TypeFluid)
{

  pot.clear();
  knFlipComputePotentialWaveCrest(
      pot, flags, v, radius, normal, tauMin, tauMax, scaleFromManta, itype, jtype);
}
static PyObject *_W_9(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipComputePotentialWaveCrest", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &pot = *_args.getPtr<Grid<Real>>("pot", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      const MACGrid &v = *_args.getPtr<MACGrid>("v", 2, &_lock);
      const int radius = _args.get<int>("radius", 3, &_lock);
      Grid<Vec3> &normal = *_args.getPtr<Grid<Vec3>>("normal", 4, &_lock);
      const Real tauMin = _args.get<Real>("tauMin", 5, &_lock);
      const Real tauMax = _args.get<Real>("tauMax", 6, &_lock);
      const Real scaleFromManta = _args.get<Real>("scaleFromManta", 7, &_lock);
      const int itype = _args.getOpt<int>("itype", 8, FlagGrid::TypeFluid, &_lock);
      const int jtype = _args.getOpt<int>("jtype", 9, FlagGrid::TypeFluid, &_lock);
      _retval = getPyNone();
      flipComputePotentialWaveCrest(
          pot, flags, v, radius, normal, tauMin, tauMax, scaleFromManta, itype, jtype);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipComputePotentialWaveCrest", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipComputePotentialWaveCrest", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipComputePotentialWaveCrest("",
                                                            "flipComputePotentialWaveCrest",
                                                            _W_9);
extern "C" {
void PbRegister_flipComputePotentialWaveCrest()
{
  KEEP_UNUSED(_RP_flipComputePotentialWaveCrest);
}
}

// LEGACY METHOD! Use flipComputeSecondaryParticlePotentials instead!
// computes normal grid &normal as gradient of levelset &phi and normalizes it

struct knFlipComputeSurfaceNormals : public KernelBase {
  knFlipComputeSurfaceNormals(Grid<Vec3> &normal, const Grid<Real> &phi)
      : KernelBase(&normal, 0), normal(normal), phi(phi)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, Grid<Vec3> &normal, const Grid<Real> &phi) const
  {
    normal[idx] = getNormalized(normal[idx]);
  }
  inline Grid<Vec3> &getArg0()
  {
    return normal;
  }
  typedef Grid<Vec3> type0;
  inline const Grid<Real> &getArg1()
  {
    return phi;
  }
  typedef Grid<Real> type1;
  void runMessage()
  {
    debMsg("Executing kernel knFlipComputeSurfaceNormals ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, normal, phi);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  Grid<Vec3> &normal;
  const Grid<Real> &phi;
};

void flipComputeSurfaceNormals(Grid<Vec3> &normal, const Grid<Real> &phi)
{
  GradientOp(normal, phi);
  knFlipComputeSurfaceNormals(normal, phi);
}
static PyObject *_W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipComputeSurfaceNormals", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Vec3> &normal = *_args.getPtr<Grid<Vec3>>("normal", 0, &_lock);
      const Grid<Real> &phi = *_args.getPtr<Grid<Real>>("phi", 1, &_lock);
      _retval = getPyNone();
      flipComputeSurfaceNormals(normal, phi);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipComputeSurfaceNormals", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipComputeSurfaceNormals", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipComputeSurfaceNormals("", "flipComputeSurfaceNormals", _W_10);
extern "C" {
void PbRegister_flipComputeSurfaceNormals()
{
  KEEP_UNUSED(_RP_flipComputeSurfaceNormals);
}
}

// LEGACY METHOD! Use flipComputeSecondaryParticlePotentials instead!
// computes the neighbor ratio for every fluid cell in &flags as the number of fluid neighbors over
// the maximum possible number of fluid neighbors

struct knFlipUpdateNeighborRatio : public KernelBase {
  knFlipUpdateNeighborRatio(const FlagGrid &flags,
                            Grid<Real> &neighborRatio,
                            const int radius,
                            const int itype = FlagGrid::TypeFluid,
                            const int jtype = FlagGrid::TypeObstacle)
      : KernelBase(&flags, 1),
        flags(flags),
        neighborRatio(neighborRatio),
        radius(radius),
        itype(itype),
        jtype(jtype)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 Grid<Real> &neighborRatio,
                 const int radius,
                 const int itype = FlagGrid::TypeFluid,
                 const int jtype = FlagGrid::TypeObstacle) const
  {

    if (!(flags(i, j, k) & itype))
      return;

    int countFluid = 0;
    int countMaxFluid = 0;
    for (IndexInt x = i - radius; x <= i + radius; x++) {
      for (IndexInt y = j - radius; y <= j + radius; y++) {
        for (IndexInt z = k - radius; z <= k + radius; z++) {
          if ((x == i && y == j && z == k) || (flags(x, y, z) & jtype))
            continue;
          if (flags(x, y, z) & itype) {
            countFluid++;
            countMaxFluid++;
          }
          else {
            countMaxFluid++;
          }
        }
      }
    }
    neighborRatio(i, j, k) = float(countFluid) / float(countMaxFluid);
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<Real> &getArg1()
  {
    return neighborRatio;
  }
  typedef Grid<Real> type1;
  inline const int &getArg2()
  {
    return radius;
  }
  typedef int type2;
  inline const int &getArg3()
  {
    return itype;
  }
  typedef int type3;
  inline const int &getArg4()
  {
    return jtype;
  }
  typedef int type4;
  void runMessage()
  {
    debMsg("Executing kernel knFlipUpdateNeighborRatio ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, neighborRatio, radius, itype, jtype);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 1; i < _maxX; i++)
          op(i, j, k, flags, neighborRatio, radius, itype, jtype);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
  }
  const FlagGrid &flags;
  Grid<Real> &neighborRatio;
  const int radius;
  const int itype;
  const int jtype;
};

void flipUpdateNeighborRatio(const FlagGrid &flags,
                             Grid<Real> &neighborRatio,
                             const int radius,
                             const int itype = FlagGrid::TypeFluid,
                             const int jtype = FlagGrid::TypeObstacle)
{

  neighborRatio.clear();
  knFlipUpdateNeighborRatio(flags, neighborRatio, radius, itype, jtype);
}
static PyObject *_W_11(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
{
  try {
    PbArgs _args(_linargs, _kwds);
    FluidSolver *parent = _args.obtainParent();
    bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
    pbPreparePlugin(parent, "flipUpdateNeighborRatio", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      Grid<Real> &neighborRatio = *_args.getPtr<Grid<Real>>("neighborRatio", 1, &_lock);
      const int radius = _args.get<int>("radius", 2, &_lock);
      const int itype = _args.getOpt<int>("itype", 3, FlagGrid::TypeFluid, &_lock);
      const int jtype = _args.getOpt<int>("jtype", 4, FlagGrid::TypeObstacle, &_lock);
      _retval = getPyNone();
      flipUpdateNeighborRatio(flags, neighborRatio, radius, itype, jtype);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipUpdateNeighborRatio", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipUpdateNeighborRatio", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipUpdateNeighborRatio("", "flipUpdateNeighborRatio", _W_11);
extern "C" {
void PbRegister_flipUpdateNeighborRatio()
{
  KEEP_UNUSED(_RP_flipUpdateNeighborRatio);
}
}

//----------------------------------------------------------------------------------------------------------------------------------------------------
// Legacy Methods (still useful for debugging)
//----------------------------------------------------------------------------------------------------------------------------------------------------
#pragma endregion

}  // namespace Manta