// 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
 *
 * Functions for calculating wavelet turbulence,
 * plus helpers to compute vorticity, and strain rate magnitude
 *
 ******************************************************************************/

#include "vectorbase.h"
#include "shapes.h"
#include "commonkernels.h"
#include "noisefield.h"

using namespace std;

namespace Manta {

//*****************************************************************************

// first some fairly generic interpolation functions for grids with multiple sizes

//! same as in grid.h , but takes an additional optional "desired" size
inline void calcGridSizeFactorMod(
    Vec3i s1, Vec3i s2, Vec3i optSize, Vec3 scale, Vec3 &sourceFactor, Vec3 &retOff)
{
  for (int c = 0; c < 3; c++) {
    if (optSize[c] > 0) {
      s2[c] = optSize[c];
    }
  }
  sourceFactor = calcGridSizeFactor(s1, s2) / scale;
  retOff = -retOff * sourceFactor + sourceFactor * 0.5;
}

void interpolateGrid(Grid<Real> &target,
                     const Grid<Real> &source,
                     Vec3 scale = Vec3(1.),
                     Vec3 offset = Vec3(0.),
                     Vec3i size = Vec3i(-1, -1, -1),
                     int orderSpace = 1)
{
  Vec3 sourceFactor(1.), off2 = offset;
  calcGridSizeFactorMod(source.getSize(), target.getSize(), size, scale, sourceFactor, off2);

  // a brief note on a mantaflow specialty: the target grid has to be the first argument here!
  // the parent fluidsolver object is taken from the first grid, and it determines the size of the
  // loop for the kernel call. as we're writing into target, it's important to loop exactly over
  // all cells of the target grid... (note, when calling the plugin in python, it doesnt matter
  // anymore).

  // sourceFactor offset necessary to shift eval points by half a small cell width
  knInterpolateGridTempl<Real>(target, source, sourceFactor, off2, orderSpace);
}
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, "interpolateGrid", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &target = *_args.getPtr<Grid<Real>>("target", 0, &_lock);
      const Grid<Real> &source = *_args.getPtr<Grid<Real>>("source", 1, &_lock);
      Vec3 scale = _args.getOpt<Vec3>("scale", 2, Vec3(1.), &_lock);
      Vec3 offset = _args.getOpt<Vec3>("offset", 3, Vec3(0.), &_lock);
      Vec3i size = _args.getOpt<Vec3i>("size", 4, Vec3i(-1, -1, -1), &_lock);
      int orderSpace = _args.getOpt<int>("orderSpace", 5, 1, &_lock);
      _retval = getPyNone();
      interpolateGrid(target, source, scale, offset, size, orderSpace);
      _args.check();
    }
    pbFinalizePlugin(parent, "interpolateGrid", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("interpolateGrid", e.what());
    return 0;
  }
}
static const Pb::Register _RP_interpolateGrid("", "interpolateGrid", _W_0);
extern "C" {
void PbRegister_interpolateGrid()
{
  KEEP_UNUSED(_RP_interpolateGrid);
}
}

void interpolateGridVec3(Grid<Vec3> &target,
                         const Grid<Vec3> &source,
                         Vec3 scale = Vec3(1.),
                         Vec3 offset = Vec3(0.),
                         Vec3i size = Vec3i(-1, -1, -1),
                         int orderSpace = 1)
{
  Vec3 sourceFactor(1.), off2 = offset;
  calcGridSizeFactorMod(source.getSize(), target.getSize(), size, scale, sourceFactor, off2);
  knInterpolateGridTempl<Vec3>(target, source, sourceFactor, off2, orderSpace);
}
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, "interpolateGridVec3", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Vec3> &target = *_args.getPtr<Grid<Vec3>>("target", 0, &_lock);
      const Grid<Vec3> &source = *_args.getPtr<Grid<Vec3>>("source", 1, &_lock);
      Vec3 scale = _args.getOpt<Vec3>("scale", 2, Vec3(1.), &_lock);
      Vec3 offset = _args.getOpt<Vec3>("offset", 3, Vec3(0.), &_lock);
      Vec3i size = _args.getOpt<Vec3i>("size", 4, Vec3i(-1, -1, -1), &_lock);
      int orderSpace = _args.getOpt<int>("orderSpace", 5, 1, &_lock);
      _retval = getPyNone();
      interpolateGridVec3(target, source, scale, offset, size, orderSpace);
      _args.check();
    }
    pbFinalizePlugin(parent, "interpolateGridVec3", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("interpolateGridVec3", e.what());
    return 0;
  }
}
static const Pb::Register _RP_interpolateGridVec3("", "interpolateGridVec3", _W_1);
extern "C" {
void PbRegister_interpolateGridVec3()
{
  KEEP_UNUSED(_RP_interpolateGridVec3);
}
}

//! interpolate a mac velocity grid from one size to another size

struct KnInterpolateMACGrid : public KernelBase {
  KnInterpolateMACGrid(MACGrid &target,
                       const MACGrid &source,
                       const Vec3 &sourceFactor,
                       const Vec3 &off,
                       int orderSpace)
      : KernelBase(&target, 0),
        target(target),
        source(source),
        sourceFactor(sourceFactor),
        off(off),
        orderSpace(orderSpace)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 MACGrid &target,
                 const MACGrid &source,
                 const Vec3 &sourceFactor,
                 const Vec3 &off,
                 int orderSpace) const
  {
    Vec3 pos = Vec3(i, j, k) * sourceFactor + off;

    Real vx = source.getInterpolatedHi(pos - Vec3(0.5, 0, 0), orderSpace)[0];
    Real vy = source.getInterpolatedHi(pos - Vec3(0, 0.5, 0), orderSpace)[1];
    Real vz = 0.f;
    if (source.is3D())
      vz = source.getInterpolatedHi(pos - Vec3(0, 0, 0.5), orderSpace)[2];

    target(i, j, k) = Vec3(vx, vy, vz);
  }
  inline MACGrid &getArg0()
  {
    return target;
  }
  typedef MACGrid type0;
  inline const MACGrid &getArg1()
  {
    return source;
  }
  typedef MACGrid type1;
  inline const Vec3 &getArg2()
  {
    return sourceFactor;
  }
  typedef Vec3 type2;
  inline const Vec3 &getArg3()
  {
    return off;
  }
  typedef Vec3 type3;
  inline int &getArg4()
  {
    return orderSpace;
  }
  typedef int type4;
  void runMessage()
  {
    debMsg("Executing kernel KnInterpolateMACGrid ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 0; j < _maxY; j++)
          for (int i = 0; i < _maxX; i++)
            op(i, j, k, target, source, sourceFactor, off, orderSpace);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 0; i < _maxX; i++)
          op(i, j, k, target, source, sourceFactor, off, orderSpace);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  MACGrid &target;
  const MACGrid &source;
  const Vec3 &sourceFactor;
  const Vec3 &off;
  int orderSpace;
};

void interpolateMACGrid(MACGrid &target,
                        const MACGrid &source,
                        Vec3 scale = Vec3(1.),
                        Vec3 offset = Vec3(0.),
                        Vec3i size = Vec3i(-1, -1, -1),
                        int orderSpace = 1)
{
  Vec3 sourceFactor(1.), off2 = offset;
  calcGridSizeFactorMod(source.getSize(), target.getSize(), size, scale, sourceFactor, off2);
  KnInterpolateMACGrid(target, source, sourceFactor, off2, orderSpace);
}
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, "interpolateMACGrid", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      MACGrid &target = *_args.getPtr<MACGrid>("target", 0, &_lock);
      const MACGrid &source = *_args.getPtr<MACGrid>("source", 1, &_lock);
      Vec3 scale = _args.getOpt<Vec3>("scale", 2, Vec3(1.), &_lock);
      Vec3 offset = _args.getOpt<Vec3>("offset", 3, Vec3(0.), &_lock);
      Vec3i size = _args.getOpt<Vec3i>("size", 4, Vec3i(-1, -1, -1), &_lock);
      int orderSpace = _args.getOpt<int>("orderSpace", 5, 1, &_lock);
      _retval = getPyNone();
      interpolateMACGrid(target, source, scale, offset, size, orderSpace);
      _args.check();
    }
    pbFinalizePlugin(parent, "interpolateMACGrid", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("interpolateMACGrid", e.what());
    return 0;
  }
}
static const Pb::Register _RP_interpolateMACGrid("", "interpolateMACGrid", _W_2);
extern "C" {
void PbRegister_interpolateMACGrid()
{
  KEEP_UNUSED(_RP_interpolateMACGrid);
}
}

//*****************************************************************************

//! Apply vector noise to grid, this is a simplified version - no position scaling or UVs

struct knApplySimpleNoiseVec3 : public KernelBase {
  knApplySimpleNoiseVec3(const FlagGrid &flags,
                         Grid<Vec3> &target,
                         const WaveletNoiseField &noise,
                         Real scale,
                         const Grid<Real> *weight)
      : KernelBase(&flags, 0),
        flags(flags),
        target(target),
        noise(noise),
        scale(scale),
        weight(weight)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 Grid<Vec3> &target,
                 const WaveletNoiseField &noise,
                 Real scale,
                 const Grid<Real> *weight) const
  {
    if (!flags.isFluid(i, j, k))
      return;
    Real factor = 1;
    if (weight)
      factor = (*weight)(i, j, k);
    target(i, j, k) += noise.evaluateCurl(Vec3(i, j, k) + Vec3(0.5)) * scale * factor;
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<Vec3> &getArg1()
  {
    return target;
  }
  typedef Grid<Vec3> type1;
  inline const WaveletNoiseField &getArg2()
  {
    return noise;
  }
  typedef WaveletNoiseField type2;
  inline Real &getArg3()
  {
    return scale;
  }
  typedef Real type3;
  inline const Grid<Real> *getArg4()
  {
    return weight;
  }
  typedef Grid<Real> type4;
  void runMessage()
  {
    debMsg("Executing kernel knApplySimpleNoiseVec3 ", 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, target, noise, scale, weight);
    }
    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, target, noise, scale, weight);
    }
  }
  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 FlagGrid &flags;
  Grid<Vec3> &target;
  const WaveletNoiseField &noise;
  Real scale;
  const Grid<Real> *weight;
};

void applySimpleNoiseVec3(const FlagGrid &flags,
                          Grid<Vec3> &target,
                          const WaveletNoiseField &noise,
                          Real scale = 1.0,
                          const Grid<Real> *weight = NULL)
{
  // note - passing a MAC grid here is slightly inaccurate, we should evaluate each component
  // separately
  knApplySimpleNoiseVec3(flags, target, noise, scale, weight);
}
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, "applySimpleNoiseVec3", !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 WaveletNoiseField &noise = *_args.getPtr<WaveletNoiseField>("noise", 2, &_lock);
      Real scale = _args.getOpt<Real>("scale", 3, 1.0, &_lock);
      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 4, NULL, &_lock);
      _retval = getPyNone();
      applySimpleNoiseVec3(flags, target, noise, scale, weight);
      _args.check();
    }
    pbFinalizePlugin(parent, "applySimpleNoiseVec3", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("applySimpleNoiseVec3", e.what());
    return 0;
  }
}
static const Pb::Register _RP_applySimpleNoiseVec3("", "applySimpleNoiseVec3", _W_3);
extern "C" {
void PbRegister_applySimpleNoiseVec3()
{
  KEEP_UNUSED(_RP_applySimpleNoiseVec3);
}
}

//! Simple noise for a real grid , follows applySimpleNoiseVec3

struct knApplySimpleNoiseReal : public KernelBase {
  knApplySimpleNoiseReal(const FlagGrid &flags,
                         Grid<Real> &target,
                         const WaveletNoiseField &noise,
                         Real scale,
                         const Grid<Real> *weight)
      : KernelBase(&flags, 0),
        flags(flags),
        target(target),
        noise(noise),
        scale(scale),
        weight(weight)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 Grid<Real> &target,
                 const WaveletNoiseField &noise,
                 Real scale,
                 const Grid<Real> *weight) const
  {
    if (!flags.isFluid(i, j, k))
      return;
    Real factor = 1;
    if (weight)
      factor = (*weight)(i, j, k);
    target(i, j, k) += noise.evaluate(Vec3(i, j, k) + Vec3(0.5)) * scale * factor;
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<Real> &getArg1()
  {
    return target;
  }
  typedef Grid<Real> type1;
  inline const WaveletNoiseField &getArg2()
  {
    return noise;
  }
  typedef WaveletNoiseField type2;
  inline Real &getArg3()
  {
    return scale;
  }
  typedef Real type3;
  inline const Grid<Real> *getArg4()
  {
    return weight;
  }
  typedef Grid<Real> type4;
  void runMessage()
  {
    debMsg("Executing kernel knApplySimpleNoiseReal ", 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, target, noise, scale, weight);
    }
    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, target, noise, scale, weight);
    }
  }
  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 FlagGrid &flags;
  Grid<Real> &target;
  const WaveletNoiseField &noise;
  Real scale;
  const Grid<Real> *weight;
};

void applySimpleNoiseReal(const FlagGrid &flags,
                          Grid<Real> &target,
                          const WaveletNoiseField &noise,
                          Real scale = 1.0,
                          const Grid<Real> *weight = NULL)
{
  knApplySimpleNoiseReal(flags, target, noise, scale, weight);
}
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, "applySimpleNoiseReal", !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 WaveletNoiseField &noise = *_args.getPtr<WaveletNoiseField>("noise", 2, &_lock);
      Real scale = _args.getOpt<Real>("scale", 3, 1.0, &_lock);
      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 4, NULL, &_lock);
      _retval = getPyNone();
      applySimpleNoiseReal(flags, target, noise, scale, weight);
      _args.check();
    }
    pbFinalizePlugin(parent, "applySimpleNoiseReal", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("applySimpleNoiseReal", e.what());
    return 0;
  }
}
static const Pb::Register _RP_applySimpleNoiseReal("", "applySimpleNoiseReal", _W_4);
extern "C" {
void PbRegister_applySimpleNoiseReal()
{
  KEEP_UNUSED(_RP_applySimpleNoiseReal);
}
}

//! Apply vector-based wavelet noise to target grid
//! This is the version with more functionality - supports uv grids, and on-the-fly interpolation
//! of input grids.

struct knApplyNoiseVec3 : public KernelBase {
  knApplyNoiseVec3(const FlagGrid &flags,
                   Grid<Vec3> &target,
                   const WaveletNoiseField &noise,
                   Real scale,
                   Real scaleSpatial,
                   const Grid<Real> *weight,
                   const Grid<Vec3> *uv,
                   bool uvInterpol,
                   const Vec3 &sourceFactor)
      : KernelBase(&flags, 0),
        flags(flags),
        target(target),
        noise(noise),
        scale(scale),
        scaleSpatial(scaleSpatial),
        weight(weight),
        uv(uv),
        uvInterpol(uvInterpol),
        sourceFactor(sourceFactor)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 Grid<Vec3> &target,
                 const WaveletNoiseField &noise,
                 Real scale,
                 Real scaleSpatial,
                 const Grid<Real> *weight,
                 const Grid<Vec3> *uv,
                 bool uvInterpol,
                 const Vec3 &sourceFactor) const
  {
    if (!flags.isFluid(i, j, k))
      return;

    // get weighting, interpolate if necessary
    Real w = 1;
    if (weight) {
      if (!uvInterpol) {
        w = (*weight)(i, j, k);
      }
      else {
        w = weight->getInterpolated(Vec3(i, j, k) * sourceFactor);
      }
    }

    // compute position where to evaluate the noise
    Vec3 pos = Vec3(i, j, k) + Vec3(0.5);
    if (uv) {
      if (!uvInterpol) {
        pos = (*uv)(i, j, k);
      }
      else {
        pos = uv->getInterpolated(Vec3(i, j, k) * sourceFactor);
        // uv coordinates are in local space - so we need to adjust the values of the positions
        pos /= sourceFactor;
      }
    }
    pos *= scaleSpatial;

    Vec3 noiseVec3 = noise.evaluateCurl(pos) * scale * w;
    // noiseVec3=pos; // debug , show interpolated positions
    target(i, j, k) += noiseVec3;
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<Vec3> &getArg1()
  {
    return target;
  }
  typedef Grid<Vec3> type1;
  inline const WaveletNoiseField &getArg2()
  {
    return noise;
  }
  typedef WaveletNoiseField type2;
  inline Real &getArg3()
  {
    return scale;
  }
  typedef Real type3;
  inline Real &getArg4()
  {
    return scaleSpatial;
  }
  typedef Real type4;
  inline const Grid<Real> *getArg5()
  {
    return weight;
  }
  typedef Grid<Real> type5;
  inline const Grid<Vec3> *getArg6()
  {
    return uv;
  }
  typedef Grid<Vec3> type6;
  inline bool &getArg7()
  {
    return uvInterpol;
  }
  typedef bool type7;
  inline const Vec3 &getArg8()
  {
    return sourceFactor;
  }
  typedef Vec3 type8;
  void runMessage()
  {
    debMsg("Executing kernel knApplyNoiseVec3 ", 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,
               target,
               noise,
               scale,
               scaleSpatial,
               weight,
               uv,
               uvInterpol,
               sourceFactor);
    }
    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,
             target,
             noise,
             scale,
             scaleSpatial,
             weight,
             uv,
             uvInterpol,
             sourceFactor);
    }
  }
  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 FlagGrid &flags;
  Grid<Vec3> &target;
  const WaveletNoiseField &noise;
  Real scale;
  Real scaleSpatial;
  const Grid<Real> *weight;
  const Grid<Vec3> *uv;
  bool uvInterpol;
  const Vec3 &sourceFactor;
};

void applyNoiseVec3(const FlagGrid &flags,
                    Grid<Vec3> &target,
                    const WaveletNoiseField &noise,
                    Real scale = 1.0,
                    Real scaleSpatial = 1.0,
                    const Grid<Real> *weight = NULL,
                    const Grid<Vec3> *uv = NULL)
{
  // check whether the uv grid has a different resolution
  bool uvInterpol = false;
  // and pre-compute conversion (only used if uvInterpol==true)
  // used for both uv and weight grid...
  Vec3 sourceFactor = Vec3(1.);
  if (uv) {
    uvInterpol = (target.getSize() != uv->getSize());
    sourceFactor = calcGridSizeFactor(uv->getSize(), target.getSize());
  }
  else if (weight) {
    uvInterpol = (target.getSize() != weight->getSize());
    sourceFactor = calcGridSizeFactor(weight->getSize(), target.getSize());
  }
  if (uv && weight)
    assertMsg(uv->getSize() == weight->getSize(), "UV and weight grid have to match!");

  // note - passing a MAC grid here is slightly inaccurate, we should evaluate each component
  // separately
  knApplyNoiseVec3(
      flags, target, noise, scale, scaleSpatial, weight, uv, uvInterpol, sourceFactor);
}
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, "applyNoiseVec3", !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 WaveletNoiseField &noise = *_args.getPtr<WaveletNoiseField>("noise", 2, &_lock);
      Real scale = _args.getOpt<Real>("scale", 3, 1.0, &_lock);
      Real scaleSpatial = _args.getOpt<Real>("scaleSpatial", 4, 1.0, &_lock);
      const Grid<Real> *weight = _args.getPtrOpt<Grid<Real>>("weight", 5, NULL, &_lock);
      const Grid<Vec3> *uv = _args.getPtrOpt<Grid<Vec3>>("uv", 6, NULL, &_lock);
      _retval = getPyNone();
      applyNoiseVec3(flags, target, noise, scale, scaleSpatial, weight, uv);
      _args.check();
    }
    pbFinalizePlugin(parent, "applyNoiseVec3", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("applyNoiseVec3", e.what());
    return 0;
  }
}
static const Pb::Register _RP_applyNoiseVec3("", "applyNoiseVec3", _W_5);
extern "C" {
void PbRegister_applyNoiseVec3()
{
  KEEP_UNUSED(_RP_applyNoiseVec3);
}
}

//! Compute energy of a staggered velocity field (at cell center)

struct KnApplyComputeEnergy : public KernelBase {
  KnApplyComputeEnergy(const FlagGrid &flags, const MACGrid &vel, Grid<Real> &energy)
      : KernelBase(&flags, 0), flags(flags), vel(vel), energy(energy)
  {
    runMessage();
    run();
  }
  inline void op(
      int i, int j, int k, const FlagGrid &flags, const MACGrid &vel, Grid<Real> &energy) const
  {
    Real e = 0.f;
    if (flags.isFluid(i, j, k)) {
      Vec3 v = vel.getCentered(i, j, k);
      e = 0.5 * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
    }
    energy(i, j, k) = e;
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline const MACGrid &getArg1()
  {
    return vel;
  }
  typedef MACGrid type1;
  inline Grid<Real> &getArg2()
  {
    return energy;
  }
  typedef Grid<Real> type2;
  void runMessage()
  {
    debMsg("Executing kernel KnApplyComputeEnergy ", 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, vel, energy);
    }
    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, vel, energy);
    }
  }
  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 FlagGrid &flags;
  const MACGrid &vel;
  Grid<Real> &energy;
};

void computeEnergy(const FlagGrid &flags, const MACGrid &vel, Grid<Real> &energy)
{
  KnApplyComputeEnergy(flags, vel, energy);
}
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, "computeEnergy", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
      Grid<Real> &energy = *_args.getPtr<Grid<Real>>("energy", 2, &_lock);
      _retval = getPyNone();
      computeEnergy(flags, vel, energy);
      _args.check();
    }
    pbFinalizePlugin(parent, "computeEnergy", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("computeEnergy", e.what());
    return 0;
  }
}
static const Pb::Register _RP_computeEnergy("", "computeEnergy", _W_6);
extern "C" {
void PbRegister_computeEnergy()
{
  KEEP_UNUSED(_RP_computeEnergy);
}
}

void computeWaveletCoeffs(Grid<Real> &input)
{
  Grid<Real> temp1(input.getParent()), temp2(input.getParent());
  WaveletNoiseField::computeCoefficients(input, temp1, temp2);
}
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, "computeWaveletCoeffs", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      Grid<Real> &input = *_args.getPtr<Grid<Real>>("input", 0, &_lock);
      _retval = getPyNone();
      computeWaveletCoeffs(input);
      _args.check();
    }
    pbFinalizePlugin(parent, "computeWaveletCoeffs", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("computeWaveletCoeffs", e.what());
    return 0;
  }
}
static const Pb::Register _RP_computeWaveletCoeffs("", "computeWaveletCoeffs", _W_7);
extern "C" {
void PbRegister_computeWaveletCoeffs()
{
  KEEP_UNUSED(_RP_computeWaveletCoeffs);
}
}

// note - alomst the same as for vorticity confinement
void computeVorticity(const MACGrid &vel, Grid<Vec3> &vorticity, Grid<Real> *norm = NULL)
{
  Grid<Vec3> velCenter(vel.getParent());
  GetCentered(velCenter, vel);
  CurlOp(velCenter, vorticity);
  if (norm)
    GridNorm(*norm, vorticity);
}
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, "computeVorticity", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
      Grid<Vec3> &vorticity = *_args.getPtr<Grid<Vec3>>("vorticity", 1, &_lock);
      Grid<Real> *norm = _args.getPtrOpt<Grid<Real>>("norm", 2, NULL, &_lock);
      _retval = getPyNone();
      computeVorticity(vel, vorticity, norm);
      _args.check();
    }
    pbFinalizePlugin(parent, "computeVorticity", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("computeVorticity", e.what());
    return 0;
  }
}
static const Pb::Register _RP_computeVorticity("", "computeVorticity", _W_8);
extern "C" {
void PbRegister_computeVorticity()
{
  KEEP_UNUSED(_RP_computeVorticity);
}
}

// note - very similar to KnComputeProductionStrain, but for use as wavelet turb weighting

struct KnComputeStrainRateMag : public KernelBase {
  KnComputeStrainRateMag(const MACGrid &vel, const Grid<Vec3> &velCenter, Grid<Real> &prod)
      : KernelBase(&vel, 1), vel(vel), velCenter(velCenter), prod(prod)
  {
    runMessage();
    run();
  }
  inline void op(
      int i, int j, int k, const MACGrid &vel, const Grid<Vec3> &velCenter, Grid<Real> &prod) const
  {
    // compute Sij = 1/2 * (dU_i/dx_j + dU_j/dx_i)
    Vec3 diag = Vec3(vel(i + 1, j, k).x, vel(i, j + 1, k).y, 0.) - vel(i, j, k);
    if (vel.is3D())
      diag[2] += vel(i, j, k + 1).z;
    else
      diag[2] = 0.;

    Vec3 ux = 0.5 * (velCenter(i + 1, j, k) - velCenter(i - 1, j, k));
    Vec3 uy = 0.5 * (velCenter(i, j + 1, k) - velCenter(i, j - 1, k));
    Vec3 uz;
    if (vel.is3D())
      uz = 0.5 * (velCenter(i, j, k + 1) - velCenter(i, j, k - 1));

    Real S12 = 0.5 * (ux.y + uy.x);
    Real S13 = 0.5 * (ux.z + uz.x);
    Real S23 = 0.5 * (uy.z + uz.y);
    Real S2 = square(diag.x) + square(diag.y) + square(diag.z) + 2.0 * square(S12) +
              2.0 * square(S13) + 2.0 * square(S23);
    prod(i, j, k) = S2;
  }
  inline const MACGrid &getArg0()
  {
    return vel;
  }
  typedef MACGrid type0;
  inline const Grid<Vec3> &getArg1()
  {
    return velCenter;
  }
  typedef Grid<Vec3> type1;
  inline Grid<Real> &getArg2()
  {
    return prod;
  }
  typedef Grid<Real> type2;
  void runMessage()
  {
    debMsg("Executing kernel KnComputeStrainRateMag ", 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, vel, velCenter, prod);
    }
    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, vel, velCenter, prod);
    }
  }
  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 MACGrid &vel;
  const Grid<Vec3> &velCenter;
  Grid<Real> &prod;
};
void computeStrainRateMag(const MACGrid &vel, Grid<Real> &mag)
{
  Grid<Vec3> velCenter(vel.getParent());
  GetCentered(velCenter, vel);
  KnComputeStrainRateMag(vel, velCenter, mag);
}
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, "computeStrainRateMag", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
      Grid<Real> &mag = *_args.getPtr<Grid<Real>>("mag", 1, &_lock);
      _retval = getPyNone();
      computeStrainRateMag(vel, mag);
      _args.check();
    }
    pbFinalizePlugin(parent, "computeStrainRateMag", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("computeStrainRateMag", e.what());
    return 0;
  }
}
static const Pb::Register _RP_computeStrainRateMag("", "computeStrainRateMag", _W_9);
extern "C" {
void PbRegister_computeStrainRateMag()
{
  KEEP_UNUSED(_RP_computeStrainRateMag);
}
}

// extrapolate a real grid into a flagged region (based on initial flags)
// by default extrapolates from fluid to obstacle cells
template<class T>
void extrapolSimpleFlagsHelper(const FlagGrid &flags,
                               Grid<T> &val,
                               int distance = 4,
                               int flagFrom = FlagGrid::TypeFluid,
                               int flagTo = FlagGrid::TypeObstacle)
{
  Grid<int> tmp(flags.getParent());
  int dim = (flags.is3D() ? 3 : 2);
  const Vec3i nb[6] = {Vec3i(1, 0, 0),
                       Vec3i(-1, 0, 0),
                       Vec3i(0, 1, 0),
                       Vec3i(0, -1, 0),
                       Vec3i(0, 0, 1),
                       Vec3i(0, 0, -1)};

  // remove all fluid cells (set to 1)
  tmp.clear();
  bool foundTarget = false;
  FOR_IJK_BND(flags, 0)
  {
    if (flags(i, j, k) & flagFrom)
      tmp(Vec3i(i, j, k)) = 1;
    if (!foundTarget && (flags(i, j, k) & flagTo))
      foundTarget = true;
  }
  // optimization, skip extrapolation if we dont have any cells to extrapolate to
  if (!foundTarget) {
    debMsg("No target cells found, skipping extrapolation", 1);
    return;
  }

  // extrapolate for given distance
  for (int d = 1; d < 1 + distance; ++d) {

    // TODO, parallelize
    FOR_IJK_BND(flags, 1)
    {
      if (tmp(i, j, k) != 0)
        continue;
      if (!(flags(i, j, k) & flagTo))
        continue;

      // copy from initialized neighbors
      Vec3i p(i, j, k);
      int nbs = 0;
      T avgVal = 0.;
      for (int n = 0; n < 2 * dim; ++n) {
        if (tmp(p + nb[n]) == d) {
          avgVal += val(p + nb[n]);
          nbs++;
        }
      }

      if (nbs > 0) {
        tmp(p) = d + 1;
        val(p) = avgVal / nbs;
      }
    }

  }  // distance
}

void extrapolateSimpleFlags(const FlagGrid &flags,
                            GridBase *val,
                            int distance = 4,
                            int flagFrom = FlagGrid::TypeFluid,
                            int flagTo = FlagGrid::TypeObstacle)
{
  if (val->getType() & GridBase::TypeReal) {
    extrapolSimpleFlagsHelper<Real>(flags, *((Grid<Real> *)val), distance, flagFrom, flagTo);
  }
  else if (val->getType() & GridBase::TypeInt) {
    extrapolSimpleFlagsHelper<int>(flags, *((Grid<int> *)val), distance, flagFrom, flagTo);
  }
  else if (val->getType() & GridBase::TypeVec3) {
    extrapolSimpleFlagsHelper<Vec3>(flags, *((Grid<Vec3> *)val), distance, flagFrom, flagTo);
  }
  else
    errMsg("extrapolateSimpleFlags: Grid Type is not supported (only int, Real, Vec3)");
}
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, "extrapolateSimpleFlags", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
      GridBase *val = _args.getPtr<GridBase>("val", 1, &_lock);
      int distance = _args.getOpt<int>("distance", 2, 4, &_lock);
      int flagFrom = _args.getOpt<int>("flagFrom", 3, FlagGrid::TypeFluid, &_lock);
      int flagTo = _args.getOpt<int>("flagTo", 4, FlagGrid::TypeObstacle, &_lock);
      _retval = getPyNone();
      extrapolateSimpleFlags(flags, val, distance, flagFrom, flagTo);
      _args.check();
    }
    pbFinalizePlugin(parent, "extrapolateSimpleFlags", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("extrapolateSimpleFlags", e.what());
    return 0;
  }
}
static const Pb::Register _RP_extrapolateSimpleFlags("", "extrapolateSimpleFlags", _W_10);
extern "C" {
void PbRegister_extrapolateSimpleFlags()
{
  KEEP_UNUSED(_RP_extrapolateSimpleFlags);
}
}

//! convert vel to a centered grid, then compute its curl
void getCurl(const MACGrid &vel, Grid<Real> &vort, int comp)
{
  Grid<Vec3> velCenter(vel.getParent()), curl(vel.getParent());

  GetCentered(velCenter, vel);
  CurlOp(velCenter, curl);
  GetComponent(curl, vort, comp);
}
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, "getCurl", !noTiming);
    PyObject *_retval = 0;
    {
      ArgLocker _lock;
      const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
      Grid<Real> &vort = *_args.getPtr<Grid<Real>>("vort", 1, &_lock);
      int comp = _args.get<int>("comp", 2, &_lock);
      _retval = getPyNone();
      getCurl(vel, vort, comp);
      _args.check();
    }
    pbFinalizePlugin(parent, "getCurl", !noTiming);
    return _retval;
  }
  catch (std::exception &e) {
    pbSetError("getCurl", e.what());
    return 0;
  }
}
static const Pb::Register _RP_getCurl("", "getCurl", _W_11);
extern "C" {
void PbRegister_getCurl()
{
  KEEP_UNUSED(_RP_getCurl);
}
}

}  // namespace Manta