blender-archive/extern/mantaflow/preprocessed/plugin/flip.cpp

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

/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey
 *
 * This program is free software, distributed under the terms of the
 * Apache License, Version 2.0
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * FLIP (fluid implicit particles)
 * for use with particle data fields
 *
 ******************************************************************************/

#include "particle.h"
#include "grid.h"
#include "commonkernels.h"
#include "randomstream.h"
#include "levelset.h"
#include "shapes.h"
#include "matrixbase.h"

using namespace std;
namespace Manta {

// init

//! note - this is a simplified version , sampleLevelsetWithParticles has more functionality

void sampleFlagsWithParticles(const FlagGrid &flags,
                              BasicParticleSystem &parts,
                              const int discretization,
                              const Real randomness)
{
  const bool is3D = flags.is3D();
  const Real jlen = randomness / discretization;
  const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization);
  RandomStream mRand(9832);

  FOR_IJK_BND(flags, 0)
  {
    if (flags.isObstacle(i, j, k))
      continue;
    if (flags.isFluid(i, j, k)) {
      const Vec3 pos(i, j, k);
      for (int dk = 0; dk < (is3D ? discretization : 1); dk++)
        for (int dj = 0; dj < discretization; dj++)
          for (int di = 0; di < discretization; di++) {
            Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk);
            subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3());
            if (!is3D)
              subpos[2] = 0.5;
            parts.addBuffered(subpos);
          }
    }
  }
  parts.insertBufferedParticles();
}
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, "sampleFlagsWithParticles", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 1, &_lock);
      const int discretization = _args.get<int>("discretization", 2, &_lock);
      const Real randomness = _args.get<Real>("randomness", 3, &_lock);
      _retval = getPyNone();
      sampleFlagsWithParticles(flags, parts, discretization, randomness);
      _args.check();
    }
    pbFinalizePlugin(parent, "sampleFlagsWithParticles", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("sampleFlagsWithParticles", e.what());
    return 0;
  }
}
static const Pb::Register _RP_sampleFlagsWithParticles("", "sampleFlagsWithParticles", _W_0);
extern "C" {
void PbRegister_sampleFlagsWithParticles()
{
  KEEP_UNUSED(_RP_sampleFlagsWithParticles);
}
}

//! sample a level set with particles, use reset to clear the particle buffer,
//! and skipEmpty for a continuous inflow (in the latter case, only empty cells will
//! be re-filled once they empty when calling sampleLevelsetWithParticles during
//! the main loop).

void sampleLevelsetWithParticles(const LevelsetGrid &phi,
                                 const FlagGrid &flags,
                                 BasicParticleSystem &parts,
                                 const int discretization,
                                 const Real randomness,
                                 const bool reset = false,
                                 const bool refillEmpty = false,
                                 const int particleFlag = -1)
{
  const bool is3D = phi.is3D();
  const Real jlen = randomness / discretization;
  const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization);
  RandomStream mRand(9832);

  if (reset) {
    parts.clear();
    parts.doCompress();
  }

  FOR_IJK_BND(phi, 0)
  {
    if (flags.isObstacle(i, j, k))
      continue;
    if (refillEmpty && flags.isFluid(i, j, k))
      continue;
    if (phi(i, j, k) < 1.733) {
      const Vec3 pos(i, j, k);
      for (int dk = 0; dk < (is3D ? discretization : 1); dk++)
        for (int dj = 0; dj < discretization; dj++)
          for (int di = 0; di < discretization; di++) {
            Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk);
            subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3());
            if (!is3D)
              subpos[2] = 0.5;
            if (phi.getInterpolated(subpos) > 0.)
              continue;
            if (particleFlag < 0) {
              parts.addBuffered(subpos);
            }
            else {
              parts.addBuffered(subpos, particleFlag);
            }
          }
    }
  }

  parts.insertBufferedParticles();
}
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, "sampleLevelsetWithParticles", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 2, &_lock);
      const int discretization = _args.get<int>("discretization", 3, &_lock);
      const Real randomness = _args.get<Real>("randomness", 4, &_lock);
      const bool reset = _args.getOpt<bool>("reset", 5, false, &_lock);
      const bool refillEmpty = _args.getOpt<bool>("refillEmpty", 6, false, &_lock);
      const int particleFlag = _args.getOpt<int>("particleFlag", 7, -1, &_lock);
      _retval = getPyNone();
      sampleLevelsetWithParticles(
          phi, flags, parts, discretization, randomness, reset, refillEmpty, particleFlag);
      _args.check();
    }
    pbFinalizePlugin(parent, "sampleLevelsetWithParticles", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("sampleLevelsetWithParticles", e.what());
    return 0;
  }
}
static const Pb::Register _RP_sampleLevelsetWithParticles("", "sampleLevelsetWithParticles", _W_1);
extern "C" {
void PbRegister_sampleLevelsetWithParticles()
{
  KEEP_UNUSED(_RP_sampleLevelsetWithParticles);
}
}

//! sample a shape with particles, use reset to clear the particle buffer,
//! and skipEmpty for a continuous inflow (in the latter case, only empty cells will
//! be re-filled once they empty when calling sampleShapeWithParticles during
//! the main loop).

void sampleShapeWithParticles(const Shape &shape,
                              const FlagGrid &flags,
                              BasicParticleSystem &parts,
                              const int discretization,
                              const Real randomness,
                              const bool reset = false,
                              const bool refillEmpty = false,
                              const LevelsetGrid *exclude = NULL)
{
  const bool is3D = flags.is3D();
  const Real jlen = randomness / discretization;
  const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization);
  RandomStream mRand(9832);

  if (reset) {
    parts.clear();
    parts.doCompress();
  }

  FOR_IJK_BND(flags, 0)
  {
    if (flags.isObstacle(i, j, k))
      continue;
    if (refillEmpty && flags.isFluid(i, j, k))
      continue;
    const Vec3 pos(i, j, k);
    for (int dk = 0; dk < (is3D ? discretization : 1); dk++)
      for (int dj = 0; dj < discretization; dj++)
        for (int di = 0; di < discretization; di++) {
          Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk);
          subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3());
          if (!is3D)
            subpos[2] = 0.5;
          if (exclude && exclude->getInterpolated(subpos) <= 0.)
            continue;
          if (!shape.isInside(subpos))
            continue;
          parts.addBuffered(subpos);
        }
  }

  parts.insertBufferedParticles();
}
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, "sampleShapeWithParticles", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const Shape &shape = *_args.getPtr<Shape>("shape", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 2, &_lock);
      const int discretization = _args.get<int>("discretization", 3, &_lock);
      const Real randomness = _args.get<Real>("randomness", 4, &_lock);
      const bool reset = _args.getOpt<bool>("reset", 5, false, &_lock);
      const bool refillEmpty = _args.getOpt<bool>("refillEmpty", 6, false, &_lock);
      const LevelsetGrid *exclude = _args.getPtrOpt<LevelsetGrid>("exclude", 7, NULL, &_lock);
      _retval = getPyNone();
      sampleShapeWithParticles(
          shape, flags, parts, discretization, randomness, reset, refillEmpty, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "sampleShapeWithParticles", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("sampleShapeWithParticles", e.what());
    return 0;
  }
}
static const Pb::Register _RP_sampleShapeWithParticles("", "sampleShapeWithParticles", _W_2);
extern "C" {
void PbRegister_sampleShapeWithParticles()
{
  KEEP_UNUSED(_RP_sampleShapeWithParticles);
}
}

//! mark fluid cells and helpers
struct knClearFluidFlags : public KernelBase {
  knClearFluidFlags(FlagGrid &flags, int dummy = 0)
      : KernelBase(&flags, 0), flags(flags), dummy(dummy)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, FlagGrid &flags, int dummy = 0) const
  {
    if (flags.isFluid(i, j, k)) {
      flags(i, j, k) = (flags(i, j, k) | FlagGrid::TypeEmpty) & ~FlagGrid::TypeFluid;
    }
  }
  inline FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline int &getArg1()
  {
    return dummy;
  }
  typedef int type1;
  void runMessage()
  {
    debMsg("Executing kernel knClearFluidFlags ", 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, flags, dummy);
    }
    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, flags, dummy);
    }
  }
  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);
  }
  FlagGrid &flags;
  int dummy;
};

struct knSetNbObstacle : public KernelBase {
  knSetNbObstacle(FlagGrid &nflags, const FlagGrid &flags, const Grid<Real> &phiObs)
      : KernelBase(&nflags, 1), nflags(nflags), flags(flags), phiObs(phiObs)
  {
    runMessage();
    run();
  }
  inline void op(
      int i, int j, int k, FlagGrid &nflags, const FlagGrid &flags, const Grid<Real> &phiObs) const
  {
    if (phiObs(i, j, k) > 0.)
      return;
    if (flags.isEmpty(i, j, k)) {
      bool set = false;
      if ((flags.isFluid(i - 1, j, k)) && (phiObs(i + 1, j, k) <= 0.))
        set = true;
      if ((flags.isFluid(i + 1, j, k)) && (phiObs(i - 1, j, k) <= 0.))
        set = true;
      if ((flags.isFluid(i, j - 1, k)) && (phiObs(i, j + 1, k) <= 0.))
        set = true;
      if ((flags.isFluid(i, j + 1, k)) && (phiObs(i, j - 1, k) <= 0.))
        set = true;
      if (flags.is3D()) {
        if ((flags.isFluid(i, j, k - 1)) && (phiObs(i, j, k + 1) <= 0.))
          set = true;
        if ((flags.isFluid(i, j, k + 1)) && (phiObs(i, j, k - 1) <= 0.))
          set = true;
      }
      if (set)
        nflags(i, j, k) = (flags(i, j, k) | FlagGrid::TypeFluid) & ~FlagGrid::TypeEmpty;
    }
  }
  inline FlagGrid &getArg0()
  {
    return nflags;
  }
  typedef FlagGrid type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const Grid<Real> &getArg2()
  {
    return phiObs;
  }
  typedef Grid<Real> type2;
  void runMessage()
  {
    debMsg("Executing kernel knSetNbObstacle ", 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, nflags, flags, phiObs);
    }
    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, nflags, flags, phiObs);
    }
  }
  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);
  }
  FlagGrid &nflags;
  const FlagGrid &flags;
  const Grid<Real> &phiObs;
};
void markFluidCells(const BasicParticleSystem &parts,
                    FlagGrid &flags,
                    const Grid<Real> *phiObs = NULL,
                    const ParticleDataImpl<int> *ptype = NULL,
                    const int exclude = 0)
{
  // remove all fluid cells
  knClearFluidFlags(flags, 0);

  // mark all particles in flaggrid as fluid
  for (IndexInt idx = 0; idx < parts.size(); idx++) {
    if (!parts.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
      continue;
    Vec3i p = toVec3i(parts.getPos(idx));
    if (flags.isInBounds(p) && flags.isEmpty(p))
      flags(p) = (flags(p) | FlagGrid::TypeFluid) & ~FlagGrid::TypeEmpty;
  }

  // special for second order obstacle BCs, check empty cells in boundary region
  if (phiObs) {
    FlagGrid tmp(flags);
    knSetNbObstacle(tmp, flags, *phiObs);
    flags.swap(tmp);
  }
}
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, "markFluidCells", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
      FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      const Grid<Real> *phiObs = _args.getPtrOpt<Grid<Real>>("phiObs", 2, NULL, &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 3, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 4, 0, &_lock);
      _retval = getPyNone();
      markFluidCells(parts, flags, phiObs, ptype, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "markFluidCells", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("markFluidCells", e.what());
    return 0;
  }
}
static const Pb::Register _RP_markFluidCells("", "markFluidCells", _W_3);
extern "C" {
void PbRegister_markFluidCells()
{
  KEEP_UNUSED(_RP_markFluidCells);
}
}

// for testing purposes only...
void testInitGridWithPos(Grid<Real> &grid)
{
  FOR_IJK(grid)
  {
    grid(i, j, k) = norm(Vec3(i, j, k));
  }
}
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, "testInitGridWithPos", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &grid = *_args.getPtr<Grid<Real>>("grid", 0, &_lock);
      _retval = getPyNone();
      testInitGridWithPos(grid);
      _args.check();
    }
    pbFinalizePlugin(parent, "testInitGridWithPos", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("testInitGridWithPos", e.what());
    return 0;
  }
}
static const Pb::Register _RP_testInitGridWithPos("", "testInitGridWithPos", _W_4);
extern "C" {
void PbRegister_testInitGridWithPos()
{
  KEEP_UNUSED(_RP_testInitGridWithPos);
}
}

//! helper to calculate particle radius factor to cover the diagonal of a cell in 2d/3d
inline Real calculateRadiusFactor(const Grid<Real> &grid, Real factor)
{
  return (grid.is3D() ? sqrt(3.) : sqrt(2.)) *
         (factor + .01);  // note, a 1% safety factor is added here
}

//! re-sample particles based on an input levelset
// optionally skip seeding new particles in "exclude" SDF

void adjustNumber(BasicParticleSystem &parts,
                  const MACGrid &vel,
                  const FlagGrid &flags,
                  int minParticles,
                  int maxParticles,
                  const LevelsetGrid &phi,
                  Real radiusFactor = 1.,
                  Real narrowBand = -1.,
                  const Grid<Real> *exclude = NULL)
{
  // which levelset to use as threshold
  const Real SURFACE_LS = -1.0 * calculateRadiusFactor(phi, radiusFactor);
  Grid<int> tmp(vel.getParent());
  std::ostringstream out;

  // count particles in cells, and delete excess particles
  for (IndexInt idx = 0; idx < (int)parts.size(); idx++) {
    if (parts.isActive(idx)) {
      Vec3i p = toVec3i(parts.getPos(idx));
      if (!tmp.isInBounds(p)) {
        parts.kill(idx);  // out of domain, remove
        continue;
      }

      Real phiv = phi.getInterpolated(parts.getPos(idx));
      if (phiv > 0) {
        parts.kill(idx);
        continue;
      }
      if (narrowBand > 0. && phiv < -narrowBand) {
        parts.kill(idx);
        continue;
      }

      bool atSurface = false;
      if (phiv > SURFACE_LS)
        atSurface = true;
      int num = tmp(p);

      // dont delete particles in non fluid cells here, the particles are "always right"
      if (num > maxParticles && (!atSurface)) {
        parts.kill(idx);
      }
      else {
        tmp(p) = num + 1;
      }
    }
  }

  // seed new particles
  RandomStream mRand(9832);
  FOR_IJK(tmp)
  {
    int cnt = tmp(i, j, k);

    // skip cells near surface
    if (phi(i, j, k) > SURFACE_LS)
      continue;
    if (narrowBand > 0. && phi(i, j, k) < -narrowBand) {
      continue;
    }
    if (exclude && ((*exclude)(i, j, k) < 0.)) {
      continue;
    }

    if (flags.isFluid(i, j, k) && cnt < minParticles) {
      for (int m = cnt; m < minParticles; m++) {
        Vec3 pos = Vec3(i, j, k) + mRand.getVec3();
        // Vec3 pos (i + 0.5, j + 0.5, k + 0.5); // cell center
        parts.addBuffered(pos);
      }
    }
  }

  parts.doCompress();
  parts.insertBufferedParticles();
}
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, "adjustNumber", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
      int minParticles = _args.get<int>("minParticles", 3, &_lock);
      int maxParticles = _args.get<int>("maxParticles", 4, &_lock);
      const LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 5, &_lock);
      Real radiusFactor = _args.getOpt<Real>("radiusFactor", 6, 1., &_lock);
      Real narrowBand = _args.getOpt<Real>("narrowBand", 7, -1., &_lock);
      const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 8, NULL, &_lock);
      _retval = getPyNone();
      adjustNumber(
          parts, vel, flags, minParticles, maxParticles, phi, radiusFactor, narrowBand, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "adjustNumber", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("adjustNumber", e.what());
    return 0;
  }
}
static const Pb::Register _RP_adjustNumber("", "adjustNumber", _W_5);
extern "C" {
void PbRegister_adjustNumber()
{
  KEEP_UNUSED(_RP_adjustNumber);
}
}

// simple and slow helper conversion to show contents of int grids like a real grid in the ui
// (use eg to quickly display contents of the particle-index grid)

void debugIntToReal(const Grid<int> &source, Grid<Real> &dest, Real factor = 1.)
{
  FOR_IJK(source)
  {
    dest(i, j, k) = (Real)source(i, j, k) * factor;
  }
}
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, "debugIntToReal", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const Grid<int> &source = *_args.getPtr<Grid<int>>("source", 0, &_lock);
      Grid<Real> &dest = *_args.getPtr<Grid<Real>>("dest", 1, &_lock);
      Real factor = _args.getOpt<Real>("factor", 2, 1., &_lock);
      _retval = getPyNone();
      debugIntToReal(source, dest, factor);
      _args.check();
    }
    pbFinalizePlugin(parent, "debugIntToReal", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("debugIntToReal", e.what());
    return 0;
  }
}
static const Pb::Register _RP_debugIntToReal("", "debugIntToReal", _W_6);
extern "C" {
void PbRegister_debugIntToReal()
{
  KEEP_UNUSED(_RP_debugIntToReal);
}
}

// build a grid that contains indices for a particle system
// the particles in a cell i,j,k are particles[index(i,j,k)] to particles[index(i+1,j,k)-1]
// (ie,  particles[index(i+1,j,k)] already belongs to cell i+1,j,k)

void gridParticleIndex(const BasicParticleSystem &parts,
                       ParticleIndexSystem &indexSys,
                       const FlagGrid &flags,
                       Grid<int> &index,
                       Grid<int> *counter = NULL)
{
  bool delCounter = false;
  if (!counter) {
    counter = new Grid<int>(flags.getParent());
    delCounter = true;
  }
  else {
    counter->clear();
  }

  // count particles in cells, and delete excess particles
  index.clear();
  int inactive = 0;
  for (IndexInt idx = 0; idx < (IndexInt)parts.size(); idx++) {
    if (parts.isActive(idx)) {
      // check index for validity...
      Vec3i p = toVec3i(parts.getPos(idx));
      if (!index.isInBounds(p)) {
        inactive++;
        continue;
      }

      index(p)++;
    }
    else {
      inactive++;
    }
  }

  // note - this one might be smaller...
  indexSys.resize(parts.size() - inactive);

  // convert per cell number to continuous index
  IndexInt idx = 0;
  FOR_IJK(index)
  {
    int num = index(i, j, k);
    index(i, j, k) = idx;
    idx += num;
  }

  // add particles to indexed array, we still need a per cell particle counter
  for (IndexInt idx = 0; idx < (IndexInt)parts.size(); idx++) {
    if (!parts.isActive(idx))
      continue;
    Vec3i p = toVec3i(parts.getPos(idx));
    if (!index.isInBounds(p)) {
      continue;
    }

    // initialize position and index into original array
    // indexSys[ index(p)+(*counter)(p) ].pos        = parts[idx].pos;
    indexSys[index(p) + (*counter)(p)].sourceIndex = idx;
    (*counter)(p)++;
  }

  if (delCounter)
    delete counter;
}
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, "gridParticleIndex", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
      ParticleIndexSystem &indexSys = *_args.getPtr<ParticleIndexSystem>("indexSys", 1, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
      Grid<int> &index = *_args.getPtr<Grid<int>>("index", 3, &_lock);
      Grid<int> *counter = _args.getPtrOpt<Grid<int>>("counter", 4, NULL, &_lock);
      _retval = getPyNone();
      gridParticleIndex(parts, indexSys, flags, index, counter);
      _args.check();
    }
    pbFinalizePlugin(parent, "gridParticleIndex", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("gridParticleIndex", e.what());
    return 0;
  }
}
static const Pb::Register _RP_gridParticleIndex("", "gridParticleIndex", _W_7);
extern "C" {
void PbRegister_gridParticleIndex()
{
  KEEP_UNUSED(_RP_gridParticleIndex);
}
}

struct ComputeUnionLevelsetPindex : public KernelBase {
  ComputeUnionLevelsetPindex(const Grid<int> &index,
                             const BasicParticleSystem &parts,
                             const ParticleIndexSystem &indexSys,
                             LevelsetGrid &phi,
                             const Real radius,
                             const ParticleDataImpl<int> *ptype,
                             const int exclude)
      : KernelBase(&index, 0),
        index(index),
        parts(parts),
        indexSys(indexSys),
        phi(phi),
        radius(radius),
        ptype(ptype),
        exclude(exclude)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const Grid<int> &index,
                 const BasicParticleSystem &parts,
                 const ParticleIndexSystem &indexSys,
                 LevelsetGrid &phi,
                 const Real radius,
                 const ParticleDataImpl<int> *ptype,
                 const int exclude) const
  {
    const Vec3 gridPos = Vec3(i, j, k) + Vec3(0.5);  // shifted by half cell
    Real phiv = radius * 1.0;                        // outside

    int r = int(radius) + 1;
    int rZ = phi.is3D() ? r : 0;
    for (int zj = k - rZ; zj <= k + rZ; zj++)
      for (int yj = j - r; yj <= j + r; yj++)
        for (int xj = i - r; xj <= i + r; xj++) {
          if (!phi.isInBounds(Vec3i(xj, yj, zj)))
            continue;

          // note, for the particle indices in indexSys the access is periodic (ie, dont skip for
          // eg inBounds(sx,10,10)
          IndexInt isysIdxS = index.index(xj, yj, zj);
          IndexInt pStart = index(isysIdxS), pEnd = 0;
          if (phi.isInBounds(isysIdxS + 1))
            pEnd = index(isysIdxS + 1);
          else
            pEnd = indexSys.size();

          // now loop over particles in cell
          for (IndexInt p = pStart; p < pEnd; ++p) {
            const int psrc = indexSys[p].sourceIndex;
            if (ptype && ((*ptype)[psrc] & exclude))
              continue;
            const Vec3 pos = parts[psrc].pos;
            phiv = std::min(phiv, fabs(norm(gridPos - pos)) - radius);
          }
        }
    phi(i, j, k) = phiv;
  }
  inline const Grid<int> &getArg0()
  {
    return index;
  }
  typedef Grid<int> type0;
  inline const BasicParticleSystem &getArg1()
  {
    return parts;
  }
  typedef BasicParticleSystem type1;
  inline const ParticleIndexSystem &getArg2()
  {
    return indexSys;
  }
  typedef ParticleIndexSystem type2;
  inline LevelsetGrid &getArg3()
  {
    return phi;
  }
  typedef LevelsetGrid type3;
  inline const Real &getArg4()
  {
    return radius;
  }
  typedef Real type4;
  inline const ParticleDataImpl<int> *getArg5()
  {
    return ptype;
  }
  typedef ParticleDataImpl<int> type5;
  inline const int &getArg6()
  {
    return exclude;
  }
  typedef int type6;
  void runMessage()
  {
    debMsg("Executing kernel ComputeUnionLevelsetPindex ", 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, index, parts, indexSys, phi, radius, ptype, exclude);
    }
    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, index, parts, indexSys, phi, radius, ptype, exclude);
    }
  }
  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);
  }
  const Grid<int> &index;
  const BasicParticleSystem &parts;
  const ParticleIndexSystem &indexSys;
  LevelsetGrid &phi;
  const Real radius;
  const ParticleDataImpl<int> *ptype;
  const int exclude;
};

void unionParticleLevelset(const BasicParticleSystem &parts,
                           const ParticleIndexSystem &indexSys,
                           const FlagGrid &flags,
                           const Grid<int> &index,
                           LevelsetGrid &phi,
                           const Real radiusFactor = 1.,
                           const ParticleDataImpl<int> *ptype = NULL,
                           const int exclude = 0)
{
  // use half a cell diagonal as base radius
  const Real radius = 0.5 * calculateRadiusFactor(phi, radiusFactor);
  // no reset of phi necessary here
  ComputeUnionLevelsetPindex(index, parts, indexSys, phi, radius, ptype, exclude);

  phi.setBound(0.5, 0);
}
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, "unionParticleLevelset", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
      const ParticleIndexSystem &indexSys = *_args.getPtr<ParticleIndexSystem>(
          "indexSys", 1, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
      const Grid<int> &index = *_args.getPtr<Grid<int>>("index", 3, &_lock);
      LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 4, &_lock);
      const Real radiusFactor = _args.getOpt<Real>("radiusFactor", 5, 1., &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 6, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 7, 0, &_lock);
      _retval = getPyNone();
      unionParticleLevelset(parts, indexSys, flags, index, phi, radiusFactor, ptype, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "unionParticleLevelset", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("unionParticleLevelset", e.what());
    return 0;
  }
}
static const Pb::Register _RP_unionParticleLevelset("", "unionParticleLevelset", _W_8);
extern "C" {
void PbRegister_unionParticleLevelset()
{
  KEEP_UNUSED(_RP_unionParticleLevelset);
}
}

//! kernel for computing averaged particle level set weights

struct ComputeAveragedLevelsetWeight : public KernelBase {
  ComputeAveragedLevelsetWeight(const BasicParticleSystem &parts,
                                const Grid<int> &index,
                                const ParticleIndexSystem &indexSys,
                                LevelsetGrid &phi,
                                const Real radius,
                                const ParticleDataImpl<int> *ptype,
                                const int exclude,
                                Grid<Vec3> *save_pAcc = NULL,
                                Grid<Real> *save_rAcc = NULL)
      : KernelBase(&index, 0),
        parts(parts),
        index(index),
        indexSys(indexSys),
        phi(phi),
        radius(radius),
        ptype(ptype),
        exclude(exclude),
        save_pAcc(save_pAcc),
        save_rAcc(save_rAcc)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const BasicParticleSystem &parts,
                 const Grid<int> &index,
                 const ParticleIndexSystem &indexSys,
                 LevelsetGrid &phi,
                 const Real radius,
                 const ParticleDataImpl<int> *ptype,
                 const int exclude,
                 Grid<Vec3> *save_pAcc = NULL,
                 Grid<Real> *save_rAcc = NULL) const
  {
    const Vec3 gridPos = Vec3(i, j, k) + Vec3(0.5);  // shifted by half cell
    Real phiv = radius * 1.0;                        // outside

    // loop over neighborhood, similar to ComputeUnionLevelsetPindex
    const Real sradiusInv = 1. / (4. * radius * radius);
    int r = int(1. * radius) + 1;
    int rZ = phi.is3D() ? r : 0;
    // accumulators
    Real wacc = 0.;
    Vec3 pacc = Vec3(0.);
    Real racc = 0.;

    for (int zj = k - rZ; zj <= k + rZ; zj++)
      for (int yj = j - r; yj <= j + r; yj++)
        for (int xj = i - r; xj <= i + r; xj++) {
          if (!phi.isInBounds(Vec3i(xj, yj, zj)))
            continue;

          IndexInt isysIdxS = index.index(xj, yj, zj);
          IndexInt pStart = index(isysIdxS), pEnd = 0;
          if (phi.isInBounds(isysIdxS + 1))
            pEnd = index(isysIdxS + 1);
          else
            pEnd = indexSys.size();
          for (IndexInt p = pStart; p < pEnd; ++p) {
            IndexInt psrc = indexSys[p].sourceIndex;
            if (ptype && ((*ptype)[psrc] & exclude))
              continue;

            Vec3 pos = parts[psrc].pos;
            Real s = normSquare(gridPos - pos) * sradiusInv;
            // Real  w = std::max(0., cubed(1.-s) );
            Real w = std::max(0., (1. - s));  // a bit smoother
            wacc += w;
            racc += radius * w;
            pacc += pos * w;
          }
        }

    if (wacc > VECTOR_EPSILON) {
      racc /= wacc;
      pacc /= wacc;
      phiv = fabs(norm(gridPos - pacc)) - racc;

      if (save_pAcc)
        (*save_pAcc)(i, j, k) = pacc;
      if (save_rAcc)
        (*save_rAcc)(i, j, k) = racc;
    }
    phi(i, j, k) = phiv;
  }
  inline const BasicParticleSystem &getArg0()
  {
    return parts;
  }
  typedef BasicParticleSystem type0;
  inline const Grid<int> &getArg1()
  {
    return index;
  }
  typedef Grid<int> type1;
  inline const ParticleIndexSystem &getArg2()
  {
    return indexSys;
  }
  typedef ParticleIndexSystem type2;
  inline LevelsetGrid &getArg3()
  {
    return phi;
  }
  typedef LevelsetGrid type3;
  inline const Real &getArg4()
  {
    return radius;
  }
  typedef Real type4;
  inline const ParticleDataImpl<int> *getArg5()
  {
    return ptype;
  }
  typedef ParticleDataImpl<int> type5;
  inline const int &getArg6()
  {
    return exclude;
  }
  typedef int type6;
  inline Grid<Vec3> *getArg7()
  {
    return save_pAcc;
  }
  typedef Grid<Vec3> type7;
  inline Grid<Real> *getArg8()
  {
    return save_rAcc;
  }
  typedef Grid<Real> type8;
  void runMessage()
  {
    debMsg("Executing kernel ComputeAveragedLevelsetWeight ", 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, parts, index, indexSys, phi, radius, ptype, exclude, save_pAcc, save_rAcc);
    }
    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, parts, index, indexSys, phi, radius, ptype, exclude, save_pAcc, save_rAcc);
    }
  }
  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);
  }
  const BasicParticleSystem &parts;
  const Grid<int> &index;
  const ParticleIndexSystem &indexSys;
  LevelsetGrid &phi;
  const Real radius;
  const ParticleDataImpl<int> *ptype;
  const int exclude;
  Grid<Vec3> *save_pAcc;
  Grid<Real> *save_rAcc;
};

template<class T> T smoothingValue(const Grid<T> val, int i, int j, int k, T center)
{
  return val(i, j, k);
}

// smoothing, and

template<class T> struct knSmoothGrid : public KernelBase {
  knSmoothGrid(const Grid<T> &me, Grid<T> &tmp, Real factor)
      : KernelBase(&me, 1), me(me), tmp(tmp), factor(factor)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, const Grid<T> &me, Grid<T> &tmp, Real factor) const
  {
    T val = me(i, j, k) + me(i + 1, j, k) + me(i - 1, j, k) + me(i, j + 1, k) + me(i, j - 1, k);
    if (me.is3D()) {
      val += me(i, j, k + 1) + me(i, j, k - 1);
    }
    tmp(i, j, k) = val * factor;
  }
  inline const Grid<T> &getArg0()
  {
    return me;
  }
  typedef Grid<T> type0;
  inline Grid<T> &getArg1()
  {
    return tmp;
  }
  typedef Grid<T> type1;
  inline Real &getArg2()
  {
    return factor;
  }
  typedef Real type2;
  void runMessage()
  {
    debMsg("Executing kernel knSmoothGrid ", 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, me, tmp, factor);
    }
    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, me, tmp, factor);
    }
  }
  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 Grid<T> &me;
  Grid<T> &tmp;
  Real factor;
};

template<class T> struct knSmoothGridNeg : public KernelBase {
  knSmoothGridNeg(const Grid<T> &me, Grid<T> &tmp, Real factor)
      : KernelBase(&me, 1), me(me), tmp(tmp), factor(factor)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, const Grid<T> &me, Grid<T> &tmp, Real factor) const
  {
    T val = me(i, j, k) + me(i + 1, j, k) + me(i - 1, j, k) + me(i, j + 1, k) + me(i, j - 1, k);
    if (me.is3D()) {
      val += me(i, j, k + 1) + me(i, j, k - 1);
    }
    val *= factor;
    if (val < tmp(i, j, k))
      tmp(i, j, k) = val;
    else
      tmp(i, j, k) = me(i, j, k);
  }
  inline const Grid<T> &getArg0()
  {
    return me;
  }
  typedef Grid<T> type0;
  inline Grid<T> &getArg1()
  {
    return tmp;
  }
  typedef Grid<T> type1;
  inline Real &getArg2()
  {
    return factor;
  }
  typedef Real type2;
  void runMessage()
  {
    debMsg("Executing kernel knSmoothGridNeg ", 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, me, tmp, factor);
    }
    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, me, tmp, factor);
    }
  }
  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 Grid<T> &me;
  Grid<T> &tmp;
  Real factor;
};

//! Zhu & Bridson particle level set creation

void averagedParticleLevelset(const BasicParticleSystem &parts,
                              const ParticleIndexSystem &indexSys,
                              const FlagGrid &flags,
                              const Grid<int> &index,
                              LevelsetGrid &phi,
                              const Real radiusFactor = 1.,
                              const int smoothen = 1,
                              const int smoothenNeg = 1,
                              const ParticleDataImpl<int> *ptype = NULL,
                              const int exclude = 0)
{
  // use half a cell diagonal as base radius
  const Real radius = 0.5 * calculateRadiusFactor(phi, radiusFactor);
  ComputeAveragedLevelsetWeight(parts, index, indexSys, phi, radius, ptype, exclude);

  // post-process level-set
  for (int i = 0; i < std::max(smoothen, smoothenNeg); ++i) {
    LevelsetGrid tmp(flags.getParent());
    if (i < smoothen) {
      knSmoothGrid<Real>(phi, tmp, 1. / (phi.is3D() ? 7. : 5.));
      phi.swap(tmp);
    }
    if (i < smoothenNeg) {
      knSmoothGridNeg<Real>(phi, tmp, 1. / (phi.is3D() ? 7. : 5.));
      phi.swap(tmp);
    }
  }
  phi.setBound(0.5, 0);
}
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, "averagedParticleLevelset", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
      const ParticleIndexSystem &indexSys = *_args.getPtr<ParticleIndexSystem>(
          "indexSys", 1, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
      const Grid<int> &index = *_args.getPtr<Grid<int>>("index", 3, &_lock);
      LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 4, &_lock);
      const Real radiusFactor = _args.getOpt<Real>("radiusFactor", 5, 1., &_lock);
      const int smoothen = _args.getOpt<int>("smoothen", 6, 1, &_lock);
      const int smoothenNeg = _args.getOpt<int>("smoothenNeg", 7, 1, &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 8, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 9, 0, &_lock);
      _retval = getPyNone();
      averagedParticleLevelset(
          parts, indexSys, flags, index, phi, radiusFactor, smoothen, smoothenNeg, ptype, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "averagedParticleLevelset", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("averagedParticleLevelset", e.what());
    return 0;
  }
}
static const Pb::Register _RP_averagedParticleLevelset("", "averagedParticleLevelset", _W_9);
extern "C" {
void PbRegister_averagedParticleLevelset()
{
  KEEP_UNUSED(_RP_averagedParticleLevelset);
}
}

//! kernel for improvedParticleLevelset

struct correctLevelset : public KernelBase {
  correctLevelset(LevelsetGrid &phi,
                  const Grid<Vec3> &pAcc,
                  const Grid<Real> &rAcc,
                  const Real radius,
                  const Real t_low,
                  const Real t_high)
      : KernelBase(&phi, 1),
        phi(phi),
        pAcc(pAcc),
        rAcc(rAcc),
        radius(radius),
        t_low(t_low),
        t_high(t_high)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 LevelsetGrid &phi,
                 const Grid<Vec3> &pAcc,
                 const Grid<Real> &rAcc,
                 const Real radius,
                 const Real t_low,
                 const Real t_high) const
  {
    if (rAcc(i, j, k) <= VECTOR_EPSILON)
      return;  // outside nothing happens
    Real x = pAcc(i, j, k).x;

    // create jacobian of pAcc via central differences
    Matrix3x3f jacobian = Matrix3x3f(0.5 * (pAcc(i + 1, j, k).x - pAcc(i - 1, j, k).x),
                                     0.5 * (pAcc(i, j + 1, k).x - pAcc(i, j - 1, k).x),
                                     0.5 * (pAcc(i, j, k + 1).x - pAcc(i, j, k - 1).x),
                                     0.5 * (pAcc(i + 1, j, k).y - pAcc(i - 1, j, k).y),
                                     0.5 * (pAcc(i, j + 1, k).y - pAcc(i, j - 1, k).y),
                                     0.5 * (pAcc(i, j, k + 1).y - pAcc(i, j, k - 1).y),
                                     0.5 * (pAcc(i + 1, j, k).z - pAcc(i - 1, j, k).z),
                                     0.5 * (pAcc(i, j + 1, k).z - pAcc(i, j - 1, k).z),
                                     0.5 * (pAcc(i, j, k + 1).z - pAcc(i, j, k - 1).z));

    // compute largest eigenvalue of jacobian
    Vec3 EV = jacobian.eigenvalues();
    Real maxEV = std::max(std::max(EV.x, EV.y), EV.z);

    // calculate correction factor
    Real correction = 1;
    if (maxEV >= t_low) {
      Real t = (t_high - maxEV) / (t_high - t_low);
      correction = t * t * t - 3 * t * t + 3 * t;
    }
    correction = (correction < 0) ?
                     0 :
                     correction;  // enforce correction factor to [0,1] (not explicitly in paper)

    const Vec3 gridPos = Vec3(i, j, k) + Vec3(0.5);  // shifted by half cell
    const Real correctedPhi = fabs(norm(gridPos - pAcc(i, j, k))) - rAcc(i, j, k) * correction;
    phi(i, j, k) = (correctedPhi > radius) ?
                       radius :
                       correctedPhi;  // adjust too high outside values when too few particles are
                                      // nearby to make smoothing possible (not in paper)
  }
  inline LevelsetGrid &getArg0()
  {
    return phi;
  }
  typedef LevelsetGrid type0;
  inline const Grid<Vec3> &getArg1()
  {
    return pAcc;
  }
  typedef Grid<Vec3> type1;
  inline const Grid<Real> &getArg2()
  {
    return rAcc;
  }
  typedef Grid<Real> type2;
  inline const Real &getArg3()
  {
    return radius;
  }
  typedef Real type3;
  inline const Real &getArg4()
  {
    return t_low;
  }
  typedef Real type4;
  inline const Real &getArg5()
  {
    return t_high;
  }
  typedef Real type5;
  void runMessage()
  {
    debMsg("Executing kernel correctLevelset ", 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, phi, pAcc, rAcc, radius, t_low, t_high);
    }
    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, phi, pAcc, rAcc, radius, t_low, t_high);
    }
  }
  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);
  }
  LevelsetGrid &phi;
  const Grid<Vec3> &pAcc;
  const Grid<Real> &rAcc;
  const Real radius;
  const Real t_low;
  const Real t_high;
};

//! Approach from "A unified particle model for fluid-solid interactions" by Solenthaler et al. in
//! 2007

void improvedParticleLevelset(const BasicParticleSystem &parts,
                              const ParticleIndexSystem &indexSys,
                              const FlagGrid &flags,
                              const Grid<int> &index,
                              LevelsetGrid &phi,
                              const Real radiusFactor = 1.,
                              const int smoothen = 1,
                              const int smoothenNeg = 1,
                              const Real t_low = 0.4,
                              const Real t_high = 3.5,
                              const ParticleDataImpl<int> *ptype = NULL,
                              const int exclude = 0)
{
  // create temporary grids to store values from levelset weight computation
  Grid<Vec3> save_pAcc(flags.getParent());
  Grid<Real> save_rAcc(flags.getParent());

  const Real radius = 0.5 * calculateRadiusFactor(
                                phi, radiusFactor);  // use half a cell diagonal as base radius
  ComputeAveragedLevelsetWeight(
      parts, index, indexSys, phi, radius, ptype, exclude, &save_pAcc, &save_rAcc);
  correctLevelset(phi, save_pAcc, save_rAcc, radius, t_low, t_high);

  // post-process level-set
  for (int i = 0; i < std::max(smoothen, smoothenNeg); ++i) {
    LevelsetGrid tmp(flags.getParent());
    if (i < smoothen) {
      knSmoothGrid<Real>(phi, tmp, 1. / (phi.is3D() ? 7. : 5.));
      phi.swap(tmp);
    }
    if (i < smoothenNeg) {
      knSmoothGridNeg<Real>(phi, tmp, 1. / (phi.is3D() ? 7. : 5.));
      phi.swap(tmp);
    }
  }
  phi.setBound(0.5, 0);
}
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, "improvedParticleLevelset", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
      const ParticleIndexSystem &indexSys = *_args.getPtr<ParticleIndexSystem>(
          "indexSys", 1, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock);
      const Grid<int> &index = *_args.getPtr<Grid<int>>("index", 3, &_lock);
      LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 4, &_lock);
      const Real radiusFactor = _args.getOpt<Real>("radiusFactor", 5, 1., &_lock);
      const int smoothen = _args.getOpt<int>("smoothen", 6, 1, &_lock);
      const int smoothenNeg = _args.getOpt<int>("smoothenNeg", 7, 1, &_lock);
      const Real t_low = _args.getOpt<Real>("t_low", 8, 0.4, &_lock);
      const Real t_high = _args.getOpt<Real>("t_high", 9, 3.5, &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 10, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 11, 0, &_lock);
      _retval = getPyNone();
      improvedParticleLevelset(parts,
                               indexSys,
                               flags,
                               index,
                               phi,
                               radiusFactor,
                               smoothen,
                               smoothenNeg,
                               t_low,
                               t_high,
                               ptype,
                               exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "improvedParticleLevelset", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("improvedParticleLevelset", e.what());
    return 0;
  }
}
static const Pb::Register _RP_improvedParticleLevelset("", "improvedParticleLevelset", _W_10);
extern "C" {
void PbRegister_improvedParticleLevelset()
{
  KEEP_UNUSED(_RP_improvedParticleLevelset);
}
}

struct knPushOutofObs : public KernelBase {
  knPushOutofObs(BasicParticleSystem &parts,
                 const FlagGrid &flags,
                 const Grid<Real> &phiObs,
                 const Real shift,
                 const Real thresh,
                 const ParticleDataImpl<int> *ptype,
                 const int exclude)
      : KernelBase(parts.size()),
        parts(parts),
        flags(flags),
        phiObs(phiObs),
        shift(shift),
        thresh(thresh),
        ptype(ptype),
        exclude(exclude)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 BasicParticleSystem &parts,
                 const FlagGrid &flags,
                 const Grid<Real> &phiObs,
                 const Real shift,
                 const Real thresh,
                 const ParticleDataImpl<int> *ptype,
                 const int exclude) const
  {
    if (!parts.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
      return;
    Vec3i p = toVec3i(parts.getPos(idx));

    if (!flags.isInBounds(p))
      return;
    Real v = phiObs.getInterpolated(parts.getPos(idx));
    if (v < thresh) {
      Vec3 grad = getGradient(phiObs, p.x, p.y, p.z);
      if (normalize(grad) < VECTOR_EPSILON)
        return;
      parts.setPos(idx, parts.getPos(idx) + grad * (thresh - v + shift));
    }
  }
  inline BasicParticleSystem &getArg0()
  {
    return parts;
  }
  typedef BasicParticleSystem type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const Grid<Real> &getArg2()
  {
    return phiObs;
  }
  typedef Grid<Real> type2;
  inline const Real &getArg3()
  {
    return shift;
  }
  typedef Real type3;
  inline const Real &getArg4()
  {
    return thresh;
  }
  typedef Real type4;
  inline const ParticleDataImpl<int> *getArg5()
  {
    return ptype;
  }
  typedef ParticleDataImpl<int> type5;
  inline const int &getArg6()
  {
    return exclude;
  }
  typedef int type6;
  void runMessage()
  {
    debMsg("Executing kernel knPushOutofObs ", 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, parts, flags, phiObs, shift, thresh, ptype, exclude);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  BasicParticleSystem &parts;
  const FlagGrid &flags;
  const Grid<Real> &phiObs;
  const Real shift;
  const Real thresh;
  const ParticleDataImpl<int> *ptype;
  const int exclude;
};
//! push particles out of obstacle levelset

void pushOutofObs(BasicParticleSystem &parts,
                  const FlagGrid &flags,
                  const Grid<Real> &phiObs,
                  const Real shift = 0,
                  const Real thresh = 0,
                  const ParticleDataImpl<int> *ptype = NULL,
                  const int exclude = 0)
{
  knPushOutofObs(parts, flags, phiObs, shift, thresh, ptype, exclude);
}
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, "pushOutofObs", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 0, &_lock);
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
      const Grid<Real> &phiObs = *_args.getPtr<Grid<Real>>("phiObs", 2, &_lock);
      const Real shift = _args.getOpt<Real>("shift", 3, 0, &_lock);
      const Real thresh = _args.getOpt<Real>("thresh", 4, 0, &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 5, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 6, 0, &_lock);
      _retval = getPyNone();
      pushOutofObs(parts, flags, phiObs, shift, thresh, ptype, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "pushOutofObs", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("pushOutofObs", e.what());
    return 0;
  }
}
static const Pb::Register _RP_pushOutofObs("", "pushOutofObs", _W_11);
extern "C" {
void PbRegister_pushOutofObs()
{
  KEEP_UNUSED(_RP_pushOutofObs);
}
}

//******************************************************************************
// grid interpolation functions

template<class T> struct knSafeDivReal : public KernelBase {
  knSafeDivReal(Grid<T> &me, const Grid<Real> &other, Real cutoff = VECTOR_EPSILON)
      : KernelBase(&me, 0), me(me), other(other), cutoff(cutoff)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 Grid<T> &me,
                 const Grid<Real> &other,
                 Real cutoff = VECTOR_EPSILON) const
  {
    if (other[idx] < cutoff) {
      me[idx] = 0.;
    }
    else {
      T div(other[idx]);
      me[idx] = safeDivide(me[idx], div);
    }
  }
  inline Grid<T> &getArg0()
  {
    return me;
  }
  typedef Grid<T> type0;
  inline const Grid<Real> &getArg1()
  {
    return other;
  }
  typedef Grid<Real> type1;
  inline Real &getArg2()
  {
    return cutoff;
  }
  typedef Real type2;
  void runMessage()
  {
    debMsg("Executing kernel knSafeDivReal ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other, cutoff);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  Grid<T> &me;
  const Grid<Real> &other;
  Real cutoff;
};

// Set velocities on the grid from the particle system

struct knMapLinearVec3ToMACGrid : public KernelBase {
  knMapLinearVec3ToMACGrid(const BasicParticleSystem &p,
                           const FlagGrid &flags,
                           const MACGrid &vel,
                           Grid<Vec3> &tmp,
                           const ParticleDataImpl<Vec3> &pvel,
                           const ParticleDataImpl<int> *ptype,
                           const int exclude)
      : KernelBase(p.size()),
        p(p),
        flags(flags),
        vel(vel),
        tmp(tmp),
        pvel(pvel),
        ptype(ptype),
        exclude(exclude)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const BasicParticleSystem &p,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 Grid<Vec3> &tmp,
                 const ParticleDataImpl<Vec3> &pvel,
                 const ParticleDataImpl<int> *ptype,
                 const int exclude)
  {
    unusedParameter(flags);
    if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
      return;
    vel.setInterpolated(p[idx].pos, pvel[idx], &tmp[0]);
  }
  inline const BasicParticleSystem &getArg0()
  {
    return p;
  }
  typedef BasicParticleSystem type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const MACGrid &getArg2()
  {
    return vel;
  }
  typedef MACGrid type2;
  inline Grid<Vec3> &getArg3()
  {
    return tmp;
  }
  typedef Grid<Vec3> type3;
  inline const ParticleDataImpl<Vec3> &getArg4()
  {
    return pvel;
  }
  typedef ParticleDataImpl<Vec3> type4;
  inline const ParticleDataImpl<int> *getArg5()
  {
    return ptype;
  }
  typedef ParticleDataImpl<int> type5;
  inline const int &getArg6()
  {
    return exclude;
  }
  typedef int type6;
  void runMessage()
  {
    debMsg("Executing kernel knMapLinearVec3ToMACGrid ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void run()
  {
    const IndexInt _sz = size;
    for (IndexInt i = 0; i < _sz; i++)
      op(i, p, flags, vel, tmp, pvel, ptype, exclude);
  }
  const BasicParticleSystem &p;
  const FlagGrid &flags;
  const MACGrid &vel;
  Grid<Vec3> &tmp;
  const ParticleDataImpl<Vec3> &pvel;
  const ParticleDataImpl<int> *ptype;
  const int exclude;
};

// optionally , this function can use an existing vec3 grid to store the weights
// this is useful in combination with the simple extrapolation function

void mapPartsToMAC(const FlagGrid &flags,
                   MACGrid &vel,
                   MACGrid &velOld,
                   const BasicParticleSystem &parts,
                   const ParticleDataImpl<Vec3> &partVel,
                   Grid<Vec3> *weight = NULL,
                   const ParticleDataImpl<int> *ptype = NULL,
                   const int exclude = 0)
{
  // interpol -> grid. tmpgrid for particle contribution weights
  bool freeTmp = false;
  if (!weight) {
    weight = new Grid<Vec3>(flags.getParent());
    freeTmp = true;
  }
  else {
    weight->clear();  // make sure we start with a zero grid!
  }
  vel.clear();
  knMapLinearVec3ToMACGrid(parts, flags, vel, *weight, partVel, ptype, exclude);

  // stomp small values in weight to zero to prevent roundoff errors
  weight->stomp(Vec3(VECTOR_EPSILON));
  vel.safeDivide(*weight);

  // store original state
  velOld.copyFrom(vel);
  if (freeTmp)
    delete weight;
}
static PyObject *_W_12(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, "mapPartsToMAC", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
      MACGrid &velOld = *_args.getPtr<MACGrid>("velOld", 2, &_lock);
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 3, &_lock);
      const ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
          "partVel", 4, &_lock);
      Grid<Vec3> *weight = _args.getPtrOpt<Grid<Vec3>>("weight", 5, NULL, &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 6, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 7, 0, &_lock);
      _retval = getPyNone();
      mapPartsToMAC(flags, vel, velOld, parts, partVel, weight, ptype, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "mapPartsToMAC", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("mapPartsToMAC", e.what());
    return 0;
  }
}
static const Pb::Register _RP_mapPartsToMAC("", "mapPartsToMAC", _W_12);
extern "C" {
void PbRegister_mapPartsToMAC()
{
  KEEP_UNUSED(_RP_mapPartsToMAC);
}
}

template<class T> struct knMapLinear : public KernelBase {
  knMapLinear(const BasicParticleSystem &p,
              const FlagGrid &flags,
              const Grid<T> &target,
              Grid<Real> &gtmp,
              const ParticleDataImpl<T> &psource)
      : KernelBase(p.size()), p(p), flags(flags), target(target), gtmp(gtmp), psource(psource)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const BasicParticleSystem &p,
                 const FlagGrid &flags,
                 const Grid<T> &target,
                 Grid<Real> &gtmp,
                 const ParticleDataImpl<T> &psource)
  {
    unusedParameter(flags);
    if (!p.isActive(idx))
      return;
    target.setInterpolated(p[idx].pos, psource[idx], gtmp);
  }
  inline const BasicParticleSystem &getArg0()
  {
    return p;
  }
  typedef BasicParticleSystem type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const Grid<T> &getArg2()
  {
    return target;
  }
  typedef Grid<T> type2;
  inline Grid<Real> &getArg3()
  {
    return gtmp;
  }
  typedef Grid<Real> type3;
  inline const ParticleDataImpl<T> &getArg4()
  {
    return psource;
  }
  typedef ParticleDataImpl<T> type4;
  void runMessage()
  {
    debMsg("Executing kernel knMapLinear ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void run()
  {
    const IndexInt _sz = size;
    for (IndexInt i = 0; i < _sz; i++)
      op(i, p, flags, target, gtmp, psource);
  }
  const BasicParticleSystem &p;
  const FlagGrid &flags;
  const Grid<T> &target;
  Grid<Real> &gtmp;
  const ParticleDataImpl<T> &psource;
};

template<class T>
void mapLinearRealHelper(const FlagGrid &flags,
                         Grid<T> &target,
                         const BasicParticleSystem &parts,
                         const ParticleDataImpl<T> &source)
{
  Grid<Real> tmp(flags.getParent());
  target.clear();
  knMapLinear<T>(parts, flags, target, tmp, source);
  knSafeDivReal<T>(target, tmp);
}

void mapPartsToGrid(const FlagGrid &flags,
                    Grid<Real> &target,
                    const BasicParticleSystem &parts,
                    const ParticleDataImpl<Real> &source)
{
  mapLinearRealHelper<Real>(flags, target, parts, source);
}
static PyObject *_W_13(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, "mapPartsToGrid", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      Grid<Real> &target = *_args.getPtr<Grid<Real>>("target", 1, &_lock);
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 2, &_lock);
      const ParticleDataImpl<Real> &source = *_args.getPtr<ParticleDataImpl<Real>>(
          "source", 3, &_lock);
      _retval = getPyNone();
      mapPartsToGrid(flags, target, parts, source);
      _args.check();
    }
    pbFinalizePlugin(parent, "mapPartsToGrid", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("mapPartsToGrid", e.what());
    return 0;
  }
}
static const Pb::Register _RP_mapPartsToGrid("", "mapPartsToGrid", _W_13);
extern "C" {
void PbRegister_mapPartsToGrid()
{
  KEEP_UNUSED(_RP_mapPartsToGrid);
}
}

void mapPartsToGridVec3(const FlagGrid &flags,
                        Grid<Vec3> &target,
                        const BasicParticleSystem &parts,
                        const ParticleDataImpl<Vec3> &source)
{
  mapLinearRealHelper<Vec3>(flags, target, parts, source);
}
static PyObject *_W_14(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, "mapPartsToGridVec3", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      Grid<Vec3> &target = *_args.getPtr<Grid<Vec3>>("target", 1, &_lock);
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 2, &_lock);
      const ParticleDataImpl<Vec3> &source = *_args.getPtr<ParticleDataImpl<Vec3>>(
          "source", 3, &_lock);
      _retval = getPyNone();
      mapPartsToGridVec3(flags, target, parts, source);
      _args.check();
    }
    pbFinalizePlugin(parent, "mapPartsToGridVec3", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("mapPartsToGridVec3", e.what());
    return 0;
  }
}
static const Pb::Register _RP_mapPartsToGridVec3("", "mapPartsToGridVec3", _W_14);
extern "C" {
void PbRegister_mapPartsToGridVec3()
{
  KEEP_UNUSED(_RP_mapPartsToGridVec3);
}
}

// integers need "max" mode, not yet implemented
// PYTHON() void mapPartsToGridInt ( FlagGrid& flags, Grid<int >& target , BasicParticleSystem&
// parts , ParticleDataImpl<int >& source ) { 	mapLinearRealHelper<int >(flags,target,parts,source);
//}

template<class T> struct knMapFromGrid : public KernelBase {
  knMapFromGrid(const BasicParticleSystem &p, const Grid<T> &gsrc, ParticleDataImpl<T> &target)
      : KernelBase(p.size()), p(p), gsrc(gsrc), target(target)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const BasicParticleSystem &p,
                 const Grid<T> &gsrc,
                 ParticleDataImpl<T> &target) const
  {
    if (!p.isActive(idx))
      return;
    target[idx] = gsrc.getInterpolated(p[idx].pos);
  }
  inline const BasicParticleSystem &getArg0()
  {
    return p;
  }
  typedef BasicParticleSystem type0;
  inline const Grid<T> &getArg1()
  {
    return gsrc;
  }
  typedef Grid<T> type1;
  inline ParticleDataImpl<T> &getArg2()
  {
    return target;
  }
  typedef ParticleDataImpl<T> type2;
  void runMessage()
  {
    debMsg("Executing kernel knMapFromGrid ", 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, p, gsrc, target);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  const BasicParticleSystem &p;
  const Grid<T> &gsrc;
  ParticleDataImpl<T> &target;
};
void mapGridToParts(const Grid<Real> &source,
                    const BasicParticleSystem &parts,
                    ParticleDataImpl<Real> &target)
{
  knMapFromGrid<Real>(parts, source, target);
}
static PyObject *_W_15(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, "mapGridToParts", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const Grid<Real> &source = *_args.getPtr<Grid<Real>>("source", 0, &_lock);
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 1, &_lock);
      ParticleDataImpl<Real> &target = *_args.getPtr<ParticleDataImpl<Real>>("target", 2, &_lock);
      _retval = getPyNone();
      mapGridToParts(source, parts, target);
      _args.check();
    }
    pbFinalizePlugin(parent, "mapGridToParts", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("mapGridToParts", e.what());
    return 0;
  }
}
static const Pb::Register _RP_mapGridToParts("", "mapGridToParts", _W_15);
extern "C" {
void PbRegister_mapGridToParts()
{
  KEEP_UNUSED(_RP_mapGridToParts);
}
}

void mapGridToPartsVec3(const Grid<Vec3> &source,
                        const BasicParticleSystem &parts,
                        ParticleDataImpl<Vec3> &target)
{
  knMapFromGrid<Vec3>(parts, source, target);
}
static PyObject *_W_16(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, "mapGridToPartsVec3", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const Grid<Vec3> &source = *_args.getPtr<Grid<Vec3>>("source", 0, &_lock);
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 1, &_lock);
      ParticleDataImpl<Vec3> &target = *_args.getPtr<ParticleDataImpl<Vec3>>("target", 2, &_lock);
      _retval = getPyNone();
      mapGridToPartsVec3(source, parts, target);
      _args.check();
    }
    pbFinalizePlugin(parent, "mapGridToPartsVec3", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("mapGridToPartsVec3", e.what());
    return 0;
  }
}
static const Pb::Register _RP_mapGridToPartsVec3("", "mapGridToPartsVec3", _W_16);
extern "C" {
void PbRegister_mapGridToPartsVec3()
{
  KEEP_UNUSED(_RP_mapGridToPartsVec3);
}
}

// Get velocities from grid

struct knMapLinearMACGridToVec3_PIC : public KernelBase {
  knMapLinearMACGridToVec3_PIC(const BasicParticleSystem &p,
                               const FlagGrid &flags,
                               const MACGrid &vel,
                               ParticleDataImpl<Vec3> &pvel,
                               const ParticleDataImpl<int> *ptype,
                               const int exclude)
      : KernelBase(p.size()),
        p(p),
        flags(flags),
        vel(vel),
        pvel(pvel),
        ptype(ptype),
        exclude(exclude)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const BasicParticleSystem &p,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 ParticleDataImpl<Vec3> &pvel,
                 const ParticleDataImpl<int> *ptype,
                 const int exclude) const
  {
    if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
      return;
    // pure PIC
    pvel[idx] = vel.getInterpolated(p[idx].pos);
  }
  inline const BasicParticleSystem &getArg0()
  {
    return p;
  }
  typedef BasicParticleSystem type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const MACGrid &getArg2()
  {
    return vel;
  }
  typedef MACGrid type2;
  inline ParticleDataImpl<Vec3> &getArg3()
  {
    return pvel;
  }
  typedef ParticleDataImpl<Vec3> type3;
  inline const ParticleDataImpl<int> *getArg4()
  {
    return ptype;
  }
  typedef ParticleDataImpl<int> type4;
  inline const int &getArg5()
  {
    return exclude;
  }
  typedef int type5;
  void runMessage()
  {
    debMsg("Executing kernel knMapLinearMACGridToVec3_PIC ", 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, p, flags, vel, pvel, ptype, exclude);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  const BasicParticleSystem &p;
  const FlagGrid &flags;
  const MACGrid &vel;
  ParticleDataImpl<Vec3> &pvel;
  const ParticleDataImpl<int> *ptype;
  const int exclude;
};

void mapMACToParts(const FlagGrid &flags,
                   const MACGrid &vel,
                   const BasicParticleSystem &parts,
                   ParticleDataImpl<Vec3> &partVel,
                   const ParticleDataImpl<int> *ptype = NULL,
                   const int exclude = 0)
{
  knMapLinearMACGridToVec3_PIC(parts, flags, vel, partVel, ptype, exclude);
}
static PyObject *_W_17(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, "mapMACToParts", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 2, &_lock);
      ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
          "partVel", 3, &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 4, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 5, 0, &_lock);
      _retval = getPyNone();
      mapMACToParts(flags, vel, parts, partVel, ptype, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "mapMACToParts", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("mapMACToParts", e.what());
    return 0;
  }
}
static const Pb::Register _RP_mapMACToParts("", "mapMACToParts", _W_17);
extern "C" {
void PbRegister_mapMACToParts()
{
  KEEP_UNUSED(_RP_mapMACToParts);
}
}

// with flip delta interpolation

struct knMapLinearMACGridToVec3_FLIP : public KernelBase {
  knMapLinearMACGridToVec3_FLIP(const BasicParticleSystem &p,
                                const FlagGrid &flags,
                                const MACGrid &vel,
                                const MACGrid &oldVel,
                                ParticleDataImpl<Vec3> &pvel,
                                const Real flipRatio,
                                const ParticleDataImpl<int> *ptype,
                                const int exclude)
      : KernelBase(p.size()),
        p(p),
        flags(flags),
        vel(vel),
        oldVel(oldVel),
        pvel(pvel),
        flipRatio(flipRatio),
        ptype(ptype),
        exclude(exclude)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const BasicParticleSystem &p,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 const MACGrid &oldVel,
                 ParticleDataImpl<Vec3> &pvel,
                 const Real flipRatio,
                 const ParticleDataImpl<int> *ptype,
                 const int exclude) const
  {
    if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
      return;
    Vec3 v = vel.getInterpolated(p[idx].pos);
    Vec3 delta = v - oldVel.getInterpolated(p[idx].pos);
    pvel[idx] = flipRatio * (pvel[idx] + delta) + (1.0 - flipRatio) * v;
  }
  inline const BasicParticleSystem &getArg0()
  {
    return p;
  }
  typedef BasicParticleSystem type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const MACGrid &getArg2()
  {
    return vel;
  }
  typedef MACGrid type2;
  inline const MACGrid &getArg3()
  {
    return oldVel;
  }
  typedef MACGrid type3;
  inline ParticleDataImpl<Vec3> &getArg4()
  {
    return pvel;
  }
  typedef ParticleDataImpl<Vec3> type4;
  inline const Real &getArg5()
  {
    return flipRatio;
  }
  typedef Real type5;
  inline const ParticleDataImpl<int> *getArg6()
  {
    return ptype;
  }
  typedef ParticleDataImpl<int> type6;
  inline const int &getArg7()
  {
    return exclude;
  }
  typedef int type7;
  void runMessage()
  {
    debMsg("Executing kernel knMapLinearMACGridToVec3_FLIP ", 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, p, flags, vel, oldVel, pvel, flipRatio, ptype, exclude);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  const BasicParticleSystem &p;
  const FlagGrid &flags;
  const MACGrid &vel;
  const MACGrid &oldVel;
  ParticleDataImpl<Vec3> &pvel;
  const Real flipRatio;
  const ParticleDataImpl<int> *ptype;
  const int exclude;
};

void flipVelocityUpdate(const FlagGrid &flags,
                        const MACGrid &vel,
                        const MACGrid &velOld,
                        const BasicParticleSystem &parts,
                        ParticleDataImpl<Vec3> &partVel,
                        const Real flipRatio,
                        const ParticleDataImpl<int> *ptype = NULL,
                        const int exclude = 0)
{
  knMapLinearMACGridToVec3_FLIP(parts, flags, vel, velOld, partVel, flipRatio, ptype, exclude);
}
static PyObject *_W_18(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, "flipVelocityUpdate", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
      const MACGrid &velOld = *_args.getPtr<MACGrid>("velOld", 2, &_lock);
      const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 3, &_lock);
      ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
          "partVel", 4, &_lock);
      const Real flipRatio = _args.get<Real>("flipRatio", 5, &_lock);
      const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
          "ptype", 6, NULL, &_lock);
      const int exclude = _args.getOpt<int>("exclude", 7, 0, &_lock);
      _retval = getPyNone();
      flipVelocityUpdate(flags, vel, velOld, parts, partVel, flipRatio, ptype, exclude);
      _args.check();
    }
    pbFinalizePlugin(parent, "flipVelocityUpdate", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("flipVelocityUpdate", e.what());
    return 0;
  }
}
static const Pb::Register _RP_flipVelocityUpdate("", "flipVelocityUpdate", _W_18);
extern "C" {
void PbRegister_flipVelocityUpdate()
{
  KEEP_UNUSED(_RP_flipVelocityUpdate);
}
}

//******************************************************************************
// narrow band

struct knCombineVels : public KernelBase {
  knCombineVels(MACGrid &vel,
                const Grid<Vec3> &w,
                MACGrid &combineVel,
                const LevelsetGrid *phi,
                Real narrowBand,
                Real thresh)
      : KernelBase(&vel, 0),
        vel(vel),
        w(w),
        combineVel(combineVel),
        phi(phi),
        narrowBand(narrowBand),
        thresh(thresh)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 MACGrid &vel,
                 const Grid<Vec3> &w,
                 MACGrid &combineVel,
                 const LevelsetGrid *phi,
                 Real narrowBand,
                 Real thresh) const
  {
    int idx = vel.index(i, j, k);

    for (int c = 0; c < 3; ++c) {
      // Correct narrow-band FLIP
      if (phi) {
        Vec3 pos(i, j, k);
        pos[(c + 1) % 3] += Real(0.5);
        pos[(c + 2) % 3] += Real(0.5);
        Real p = phi->getInterpolated(pos);
        if (p < -narrowBand) {
          vel[idx][c] = 0;
          continue;
        }
      }

      if (w[idx][c] > thresh) {
        combineVel[idx][c] = vel[idx][c];
        vel[idx][c] = -1;
      }
      else {
        vel[idx][c] = 0;
      }
    }
  }
  inline MACGrid &getArg0()
  {
    return vel;
  }
  typedef MACGrid type0;
  inline const Grid<Vec3> &getArg1()
  {
    return w;
  }
  typedef Grid<Vec3> type1;
  inline MACGrid &getArg2()
  {
    return combineVel;
  }
  typedef MACGrid type2;
  inline const LevelsetGrid *getArg3()
  {
    return phi;
  }
  typedef LevelsetGrid type3;
  inline Real &getArg4()
  {
    return narrowBand;
  }
  typedef Real type4;
  inline Real &getArg5()
  {
    return thresh;
  }
  typedef Real type5;
  void runMessage()
  {
    debMsg("Executing kernel knCombineVels ", 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, vel, w, combineVel, phi, narrowBand, thresh);
    }
    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, vel, w, combineVel, phi, narrowBand, thresh);
    }
  }
  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 &vel;
  const Grid<Vec3> &w;
  MACGrid &combineVel;
  const LevelsetGrid *phi;
  Real narrowBand;
  Real thresh;
};

//! narrow band velocity combination

void combineGridVel(MACGrid &vel,
                    const Grid<Vec3> &weight,
                    MACGrid &combineVel,
                    const LevelsetGrid *phi = NULL,
                    Real narrowBand = 0.0,
                    Real thresh = 0.0)
{
  knCombineVels(vel, weight, combineVel, phi, narrowBand, thresh);
}
static PyObject *_W_19(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, "combineGridVel", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
      const Grid<Vec3> &weight = *_args.getPtr<Grid<Vec3>>("weight", 1, &_lock);
      MACGrid &combineVel = *_args.getPtr<MACGrid>("combineVel", 2, &_lock);
      const LevelsetGrid *phi = _args.getPtrOpt<LevelsetGrid>("phi", 3, NULL, &_lock);
      Real narrowBand = _args.getOpt<Real>("narrowBand", 4, 0.0, &_lock);
      Real thresh = _args.getOpt<Real>("thresh", 5, 0.0, &_lock);
      _retval = getPyNone();
      combineGridVel(vel, weight, combineVel, phi, narrowBand, thresh);
      _args.check();
    }
    pbFinalizePlugin(parent, "combineGridVel", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("combineGridVel", e.what());
    return 0;
  }
}
static const Pb::Register _RP_combineGridVel("", "combineGridVel", _W_19);
extern "C" {
void PbRegister_combineGridVel()
{
  KEEP_UNUSED(_RP_combineGridVel);
}
}

//! surface tension helper
void getLaplacian(Grid<Real> &laplacian, const Grid<Real> &grid)
{
  LaplaceOp(laplacian, grid);
}
static PyObject *_W_20(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, "getLaplacian", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &laplacian = *_args.getPtr<Grid<Real>>("laplacian", 0, &_lock);
      const Grid<Real> &grid = *_args.getPtr<Grid<Real>>("grid", 1, &_lock);
      _retval = getPyNone();
      getLaplacian(laplacian, grid);
      _args.check();
    }
    pbFinalizePlugin(parent, "getLaplacian", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("getLaplacian", e.what());
    return 0;
  }
}
static const Pb::Register _RP_getLaplacian("", "getLaplacian", _W_20);
extern "C" {
void PbRegister_getLaplacian()
{
  KEEP_UNUSED(_RP_getLaplacian);
}
}

void getCurvature(Grid<Real> &curv, const Grid<Real> &grid, const Real h = 1.0)
{
  CurvatureOp(curv, grid, h);
}
static PyObject *_W_21(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, "getCurvature", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &curv = *_args.getPtr<Grid<Real>>("curv", 0, &_lock);
      const Grid<Real> &grid = *_args.getPtr<Grid<Real>>("grid", 1, &_lock);
      const Real h = _args.getOpt<Real>("h", 2, 1.0, &_lock);
      _retval = getPyNone();
      getCurvature(curv, grid, h);
      _args.check();
    }
    pbFinalizePlugin(parent, "getCurvature", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("getCurvature", e.what());
    return 0;
  }
}
static const Pb::Register _RP_getCurvature("", "getCurvature", _W_21);
extern "C" {
void PbRegister_getCurvature()
{
  KEEP_UNUSED(_RP_getCurvature);
}
}

}  // namespace Manta