path: root/extern/mantaflow/preprocessed/levelset.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
 *
 * Levelset
 *
 ******************************************************************************/

#include "levelset.h"
#include "fastmarch.h"
#include "kernel.h"
#include "mcubes.h"
#include "mesh.h"
#include <stack>

using namespace std;
namespace Manta {

//************************************************************************
// Helper functions and kernels for marching

static const int FlagInited = FastMarch<FmHeapEntryOut, +1>::FlagInited;

// neighbor lookup vectors
static const Vec3i neighbors[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)};

struct InitFmIn : public KernelBase {
  InitFmIn(const FlagGrid &flags,
           Grid<int> &fmFlags,
           Grid<Real> &phi,
           bool ignoreWalls,
           int obstacleType)
      : KernelBase(&flags, 1),
        flags(flags),
        fmFlags(fmFlags),
        phi(phi),
        ignoreWalls(ignoreWalls),
        obstacleType(obstacleType)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 Grid<int> &fmFlags,
                 Grid<Real> &phi,
                 bool ignoreWalls,
                 int obstacleType) const
  {
    const IndexInt idx = flags.index(i, j, k);
    const Real v = phi[idx];
    if (ignoreWalls) {
      if (v >= 0. && ((flags[idx] & obstacleType) == 0))
        fmFlags[idx] = FlagInited;
      else
        fmFlags[idx] = 0;
    }
    else {
      if (v >= 0)
        fmFlags[idx] = FlagInited;
      else
        fmFlags[idx] = 0;
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<int> &getArg1()
  {
    return fmFlags;
  }
  typedef Grid<int> type1;
  inline Grid<Real> &getArg2()
  {
    return phi;
  }
  typedef Grid<Real> type2;
  inline bool &getArg3()
  {
    return ignoreWalls;
  }
  typedef bool type3;
  inline int &getArg4()
  {
    return obstacleType;
  }
  typedef int type4;
  void runMessage()
  {
    debMsg("Executing kernel InitFmIn ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 1; i < _maxX; i++)
          op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
  }
  const FlagGrid &flags;
  Grid<int> &fmFlags;
  Grid<Real> &phi;
  bool ignoreWalls;
  int obstacleType;
};

struct InitFmOut : public KernelBase {
  InitFmOut(const FlagGrid &flags,
            Grid<int> &fmFlags,
            Grid<Real> &phi,
            bool ignoreWalls,
            int obstacleType)
      : KernelBase(&flags, 1),
        flags(flags),
        fmFlags(fmFlags),
        phi(phi),
        ignoreWalls(ignoreWalls),
        obstacleType(obstacleType)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const FlagGrid &flags,
                 Grid<int> &fmFlags,
                 Grid<Real> &phi,
                 bool ignoreWalls,
                 int obstacleType) const
  {
    const IndexInt idx = flags.index(i, j, k);
    const Real v = phi[idx];
    if (ignoreWalls) {
      fmFlags[idx] = (v < 0) ? FlagInited : 0;
      if ((flags[idx] & obstacleType) != 0) {
        fmFlags[idx] = 0;
        phi[idx] = 0;
      }
    }
    else {
      fmFlags[idx] = (v < 0) ? FlagInited : 0;
    }
  }
  inline const FlagGrid &getArg0()
  {
    return flags;
  }
  typedef FlagGrid type0;
  inline Grid<int> &getArg1()
  {
    return fmFlags;
  }
  typedef Grid<int> type1;
  inline Grid<Real> &getArg2()
  {
    return phi;
  }
  typedef Grid<Real> type2;
  inline bool &getArg3()
  {
    return ignoreWalls;
  }
  typedef bool type3;
  inline int &getArg4()
  {
    return obstacleType;
  }
  typedef int type4;
  void runMessage()
  {
    debMsg("Executing kernel InitFmOut ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 1; i < _maxX; i++)
          op(i, j, k, flags, fmFlags, phi, ignoreWalls, obstacleType);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
  }
  const FlagGrid &flags;
  Grid<int> &fmFlags;
  Grid<Real> &phi;
  bool ignoreWalls;
  int obstacleType;
};

struct SetUninitialized : public KernelBase {
  SetUninitialized(const Grid<int> &flags,
                   Grid<int> &fmFlags,
                   Grid<Real> &phi,
                   const Real val,
                   int ignoreWalls,
                   int obstacleType)
      : KernelBase(&flags, 1),
        flags(flags),
        fmFlags(fmFlags),
        phi(phi),
        val(val),
        ignoreWalls(ignoreWalls),
        obstacleType(obstacleType)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const Grid<int> &flags,
                 Grid<int> &fmFlags,
                 Grid<Real> &phi,
                 const Real val,
                 int ignoreWalls,
                 int obstacleType) const
  {
    if (ignoreWalls) {
      if ((fmFlags(i, j, k) != FlagInited) && ((flags(i, j, k) & obstacleType) == 0)) {
        phi(i, j, k) = val;
      }
    }
    else {
      if ((fmFlags(i, j, k) != FlagInited))
        phi(i, j, k) = val;
    }
  }
  inline const Grid<int> &getArg0()
  {
    return flags;
  }
  typedef Grid<int> type0;
  inline Grid<int> &getArg1()
  {
    return fmFlags;
  }
  typedef Grid<int> type1;
  inline Grid<Real> &getArg2()
  {
    return phi;
  }
  typedef Grid<Real> type2;
  inline const Real &getArg3()
  {
    return val;
  }
  typedef Real type3;
  inline int &getArg4()
  {
    return ignoreWalls;
  }
  typedef int type4;
  inline int &getArg5()
  {
    return obstacleType;
  }
  typedef int type5;
  void runMessage()
  {
    debMsg("Executing kernel SetUninitialized ", 3);
    debMsg("Kernel range"
               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    const int _maxX = maxX;
    const int _maxY = maxY;
    if (maxZ > 1) {
      for (int k = __r.begin(); k != (int)__r.end(); k++)
        for (int j = 1; j < _maxY; j++)
          for (int i = 1; i < _maxX; i++)
            op(i, j, k, flags, fmFlags, phi, val, ignoreWalls, obstacleType);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = 1; i < _maxX; i++)
          op(i, j, k, flags, fmFlags, phi, val, ignoreWalls, obstacleType);
    }
  }
  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<int> &flags;
  Grid<int> &fmFlags;
  Grid<Real> &phi;
  const Real val;
  int ignoreWalls;
  int obstacleType;
};

template<bool inward>
inline bool isAtInterface(const Grid<int> &fmFlags, Grid<Real> &phi, const Vec3i &p)
{
  // check for interface
  int max = phi.is3D() ? 6 : 4;
  for (int nb = 0; nb < max; nb++) {
    const Vec3i pn(p + neighbors[nb]);
    if (!fmFlags.isInBounds(pn))
      continue;

    if (fmFlags(pn) != FlagInited)
      continue;
    if ((inward && phi(pn) >= 0.) || (!inward && phi(pn) < 0.))
      return true;
  }
  return false;
}

//************************************************************************
// Levelset class def

LevelsetGrid::LevelsetGrid(FluidSolver *parent, bool show) : Grid<Real>(parent, show)
{
  mType = (GridType)(TypeLevelset | TypeReal);
}

LevelsetGrid::LevelsetGrid(FluidSolver *parent, Real *data, bool show)
    : Grid<Real>(parent, data, show)
{
  mType = (GridType)(TypeLevelset | TypeReal);
}

Real LevelsetGrid::invalidTimeValue()
{
  return FastMarch<FmHeapEntryOut, 1>::InvalidTime();
}

//! Kernel: perform levelset union
struct KnJoin : public KernelBase {
  KnJoin(Grid<Real> &a, const Grid<Real> &b) : KernelBase(&a, 0), a(a), b(b)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, Grid<Real> &a, const Grid<Real> &b) const
  {
    a[idx] = min(a[idx], b[idx]);
  }
  inline Grid<Real> &getArg0()
  {
    return a;
  }
  typedef Grid<Real> type0;
  inline const Grid<Real> &getArg1()
  {
    return b;
  }
  typedef Grid<Real> type1;
  void runMessage()
  {
    debMsg("Executing kernel KnJoin ", 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, a, b);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  Grid<Real> &a;
  const Grid<Real> &b;
};
void LevelsetGrid::join(const LevelsetGrid &o)
{
  KnJoin(*this, o);
}

//! subtract b, note does not preserve SDF!
struct KnSubtract : public KernelBase {
  KnSubtract(Grid<Real> &a, const Grid<Real> &b, const FlagGrid *flags, int subtractType)
      : KernelBase(&a, 0), a(a), b(b), flags(flags), subtractType(subtractType)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 Grid<Real> &a,
                 const Grid<Real> &b,
                 const FlagGrid *flags,
                 int subtractType) const
  {
    if (flags && ((*flags)(idx)&subtractType) == 0)
      return;
    if (b[idx] < 0.)
      a[idx] = b[idx] * -1.;
  }
  inline Grid<Real> &getArg0()
  {
    return a;
  }
  typedef Grid<Real> type0;
  inline const Grid<Real> &getArg1()
  {
    return b;
  }
  typedef Grid<Real> type1;
  inline const FlagGrid *getArg2()
  {
    return flags;
  }
  typedef FlagGrid type2;
  inline int &getArg3()
  {
    return subtractType;
  }
  typedef int type3;
  void runMessage()
  {
    debMsg("Executing kernel KnSubtract ", 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, a, b, flags, subtractType);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  Grid<Real> &a;
  const Grid<Real> &b;
  const FlagGrid *flags;
  int subtractType;
};
void LevelsetGrid::subtract(const LevelsetGrid &o, const FlagGrid *flags, const int subtractType)
{
  KnSubtract(*this, o, flags, subtractType);
}

//! re-init levelset and extrapolate velocities (in & out)
//  note - uses flags to identify border (could also be done based on ls values)
static void doReinitMarch(Grid<Real> &phi,
                          const FlagGrid &flags,
                          Real maxTime,
                          MACGrid *velTransport,
                          bool ignoreWalls,
                          bool correctOuterLayer,
                          int obstacleType)
{
  const int dim = (phi.is3D() ? 3 : 2);
  Grid<int> fmFlags(phi.getParent());

  FastMarch<FmHeapEntryIn, -1> marchIn(flags, fmFlags, phi, maxTime, nullptr);

  // march inside
  InitFmIn(flags, fmFlags, phi, ignoreWalls, obstacleType);

  FOR_IJK_BND(flags, 1)
  {
    if (fmFlags(i, j, k) == FlagInited)
      continue;
    if (ignoreWalls && ((flags(i, j, k) & obstacleType) != 0))
      continue;
    const Vec3i p(i, j, k);

    if (isAtInterface<true>(fmFlags, phi, p)) {
      // set value
      fmFlags(p) = FlagInited;

      // add neighbors that are not at the interface
      for (int nb = 0; nb < 2 * dim; nb++) {
        const Vec3i pn(p + neighbors[nb]);  // index always valid due to bnd=1
        if (ignoreWalls && ((flags.get(pn) & obstacleType) != 0))
          continue;

        // check neighbors of neighbor
        if (phi(pn) < 0. && !isAtInterface<true>(fmFlags, phi, pn)) {
          marchIn.addToList(pn, p);
        }
      }
    }
  }
  marchIn.performMarching();
  // done with inwards marching

  // now march out...

  // set un initialized regions
  SetUninitialized(flags, fmFlags, phi, -maxTime - 1., ignoreWalls, obstacleType);

  InitFmOut(flags, fmFlags, phi, ignoreWalls, obstacleType);

  FastMarch<FmHeapEntryOut, +1> marchOut(flags, fmFlags, phi, maxTime, velTransport);

  // by default, correctOuterLayer is on
  if (correctOuterLayer) {
    // normal version, inwards march is done, now add all outside values (0..2] to list
    // note, this might move the interface a bit! but keeps a nice signed distance field...
    FOR_IJK_BND(flags, 1)
    {
      if (ignoreWalls && ((flags(i, j, k) & obstacleType) != 0))
        continue;
      const Vec3i p(i, j, k);

      // check nbs
      for (int nb = 0; nb < 2 * dim; nb++) {
        const Vec3i pn(p + neighbors[nb]);  // index always valid due to bnd=1

        if (fmFlags(pn) != FlagInited)
          continue;
        if (ignoreWalls && ((flags.get(pn) & obstacleType)) != 0)
          continue;

        const Real nbPhi = phi(pn);

        // only add nodes near interface, not e.g. outer boundary vs. invalid region
        if (nbPhi < 0 && nbPhi >= -2)
          marchOut.addToList(p, pn);
      }
    }
  }
  else {
    // alternative version, keep interface, do not distort outer cells
    // add all ouside values, but not those at the IF layer
    FOR_IJK_BND(flags, 1)
    {
      if (ignoreWalls && ((flags(i, j, k) & obstacleType) != 0))
        continue;

      // only look at ouside values
      const Vec3i p(i, j, k);
      if (phi(p) < 0)
        continue;

      if (isAtInterface<false>(fmFlags, phi, p)) {
        // now add all non, interface neighbors
        fmFlags(p) = FlagInited;

        // add neighbors that are not at the interface
        for (int nb = 0; nb < 2 * dim; nb++) {
          const Vec3i pn(p + neighbors[nb]);  // index always valid due to bnd=1
          if (ignoreWalls && ((flags.get(pn) & obstacleType) != 0))
            continue;

          // check neighbors of neighbor
          if (phi(pn) > 0. && !isAtInterface<false>(fmFlags, phi, pn)) {
            marchOut.addToList(pn, p);
          }
        }
      }
    }
  }
  marchOut.performMarching();

  // set un initialized regions
  SetUninitialized(flags, fmFlags, phi, +maxTime + 1., ignoreWalls, obstacleType);
}

//! call for levelset grids & external real grids

void LevelsetGrid::reinitMarching(const FlagGrid &flags,
                                  Real maxTime,
                                  MACGrid *velTransport,
                                  bool ignoreWalls,
                                  bool correctOuterLayer,
                                  int obstacleType)
{
  doReinitMarch(*this, flags, maxTime, velTransport, ignoreWalls, correctOuterLayer, obstacleType);
}

void LevelsetGrid::initFromFlags(const FlagGrid &flags, bool ignoreWalls)
{
  FOR_IDX(*this)
  {
    if (flags.isFluid(idx) || (ignoreWalls && flags.isObstacle(idx)))
      mData[idx] = -0.5;
    else
      mData[idx] = 0.5;
  }
}

/* Helper variables that are used in flood-fill functions. */
static const int ID_UNKNOWN = 0;
static const int ID_VISITED = 1;

/* Fills all cells in the target grid that have not been marked during a flood-fill. */

struct KnFillApply : public KernelBase {
  KnFillApply(Grid<Real> &target,
              Grid<int> &visited,
              const Real value,
              const int boundaryWidth,
              const bool outside)
      : KernelBase(&target, boundaryWidth),
        target(target),
        visited(visited),
        value(value),
        boundaryWidth(boundaryWidth),
        outside(outside)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 Grid<Real> &target,
                 Grid<int> &visited,
                 const Real value,
                 const int boundaryWidth,
                 const bool outside) const
  {

    if (visited(i, j, k) == ID_VISITED)
      return;
    if (outside && target(i, j, k) < 0)
      return;
    if (!outside && target(i, j, k) >= 0)
      return;

    /* Actual flood-fill override. */
    target(i, j, k) = value;
  }
  inline Grid<Real> &getArg0()
  {
    return target;
  }
  typedef Grid<Real> type0;
  inline Grid<int> &getArg1()
  {
    return visited;
  }
  typedef Grid<int> type1;
  inline const Real &getArg2()
  {
    return value;
  }
  typedef Real type2;
  inline const int &getArg3()
  {
    return boundaryWidth;
  }
  typedef int type3;
  inline const bool &getArg4()
  {
    return outside;
  }
  typedef bool type4;
  void runMessage()
  {
    debMsg("Executing kernel KnFillApply ", 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 = boundaryWidth; j < _maxY; j++)
          for (int i = boundaryWidth; i < _maxX; i++)
            op(i, j, k, target, visited, value, boundaryWidth, outside);
    }
    else {
      const int k = 0;
      for (int j = __r.begin(); j != (int)__r.end(); j++)
        for (int i = boundaryWidth; i < _maxX; i++)
          op(i, j, k, target, visited, value, boundaryWidth, outside);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(boundaryWidth, maxY), *this);
  }
  Grid<Real> &target;
  Grid<int> &visited;
  const Real value;
  const int boundaryWidth;
  const bool outside;
};

/* Basic flood fill implementation used to fill inside / outside areas of levelset.
 * Calling this function will ensure that there are no fluid cells inside obstacles.
 * I.e. starting from walls, cells will be tagged in flood-fill fashion, stopping at 0 borders.
 * All remaining cells will be filled with the fill value. Outside mode inverts search behavior. */
void LevelsetGrid::floodFill(const Real value, const bool outside, const int boundaryWidth)
{

  /* Sanity check: Filling mode and filling value need to "match". */
  if (outside) {
    assertMsg(value < 0, "Cannot fill outside with (positive) value " << value);
  }
  else {
    assertMsg(value >= 0, "Cannot fill inside with (negative) value " << value);
  }

  Grid<Real> levelsetCopy(this->getParent());
  Grid<int> visited(this->getParent());
  std::stack<Vec3i> todoPos;

  const int maxNeighbors = this->is3D() ? 6 : 4;
  const Vec3i maxSize(this->getSize() - 1);

  Vec3i bnd(2 * boundaryWidth);
  if (!this->is3D())
    bnd.z = 0;
  const int cellCntNoBnd = (this->getSizeX() - bnd.x) * (this->getSizeY() - bnd.y) *
                           (this->getSizeZ() - bnd.z);

  /* Initialize temporary helper grids. */
  levelsetCopy.copyFrom(*this);
  visited.setConst(ID_UNKNOWN);

  FOR_IJK_BND(visited, boundaryWidth)
  {

    /* Skip inside / outside cells depending on search mode. */
    if (outside && levelsetCopy(i, j, k) < 0)
      continue;
    if (!outside && levelsetCopy(i, j, k) >= 0)
      continue;
    /* Skip cell if it already has been visited. */
    if (visited(i, j, k) == ID_VISITED)
      continue;

    Vec3i c(i, j, k);

    bool isWallCell = (c.x - boundaryWidth == 0 || c.x == maxSize.x - boundaryWidth);
    isWallCell |= (c.y - boundaryWidth == 0 || c.y == maxSize.y - boundaryWidth);
    if (this->is3D())
      isWallCell |= (c.z - boundaryWidth == 0 || c.z == maxSize.z - boundaryWidth);

    /* Only start searching from borders. */
    if (!isWallCell)
      continue;

    /* Start flood-fill loop by initializing todo stack with current cell. */
    todoPos.push(c);
    visited(c) = ID_VISITED;

    while (!todoPos.empty()) {
      c = todoPos.top();
      todoPos.pop();

      /* Add all neighbor cells to search stack. */
      for (int nb = 0; nb < maxNeighbors; nb++) {
        const Vec3i neigh(c + neighbors[nb]);

        if (!visited.isInBounds(neigh, boundaryWidth))
          continue;
        /* Skip inside / outside area depening on what we search for. */
        if (outside && levelsetCopy(neigh) < 0)
          continue;
        if (!outside && levelsetCopy(neigh) >= 0)
          continue;
        /* Skip neighbor if it already has been visited. */
        if (visited(neigh) == ID_VISITED)
          continue;

        assertMsg(visited(neigh) == ID_UNKNOWN,
                  "Cell must be of type 'unknown' at this point in the loop");
        todoPos.push(neigh);
        visited(neigh) = ID_VISITED;
      }
      assertMsg(todoPos.size() <= cellCntNoBnd,
                "Flood-fill todo stack cannot be greater than domain cell count - "
                    << todoPos.size() << " vs " << cellCntNoBnd);
    }
  }
  KnFillApply(*this, visited, value, boundaryWidth, outside);
}

/* Deprecated: Use floodFill() function instead. */
void LevelsetGrid::fillHoles(int maxDepth, int boundaryWidth)
{
  debMsg("Deprecated - do not use fillHoles() ... use floodFill() instead", 1);

  Real curVal, i1, i2, j1, j2, k1, k2;
  Vec3i c, cTmp;
  std::stack<Vec3i> undoPos;
  std::stack<Real> undoVal;
  std::stack<Vec3i> todoPos;

  FOR_IJK_BND(*this, boundaryWidth)
  {

    curVal = mData[index(i, j, k)];
    i1 = mData[index(i - 1, j, k)];
    i2 = mData[index(i + 1, j, k)];
    j1 = mData[index(i, j - 1, k)];
    j2 = mData[index(i, j + 1, k)];
    k1 = mData[index(i, j, k - 1)];
    k2 = mData[index(i, j, k + 1)];

    /* Skip cells inside and cells outside with no inside neighbours early */
    if (curVal < 0.)
      continue;
    if (curVal > 0. && i1 > 0. && i2 > 0. && j1 > 0. && j2 > 0. && k1 > 0. && k2 > 0.)
      continue;

    /* Cell at c is positive (outside) and has at least one negative (inside) neighbour cell */
    c = Vec3i(i, j, k);

    /* Current cell is outside and has inside neighbour(s) */
    undoPos.push(c);
    undoVal.push(curVal);
    todoPos.push(c);

    /* Enforce negative cell - if search depth gets exceeded this will be reverted to the original
     * value */
    mData[index(c.x, c.y, c.z)] = -0.5;

    while (!todoPos.empty()) {
      todoPos.pop();

      /* Add neighbouring positive (inside) cells to stacks and set negavtive cell value */
      if (c.x > 0 && mData[index(c.x - 1, c.y, c.z)] > 0.) {
        cTmp = Vec3i(c.x - 1, c.y, c.z);
        undoPos.push(cTmp);
        undoVal.push(mData[index(cTmp)]);
        todoPos.push(cTmp);
        mData[index(cTmp)] = -0.5;
      }
      if (c.y > 0 && mData[index(c.x, c.y - 1, c.z)] > 0.) {
        cTmp = Vec3i(c.x, c.y - 1, c.z);
        undoPos.push(cTmp);
        undoVal.push(mData[index(cTmp)]);
        todoPos.push(cTmp);
        mData[index(cTmp)] = -0.5;
      }
      if (c.z > 0 && mData[index(c.x, c.y, c.z - 1)] > 0.) {
        cTmp = Vec3i(c.x, c.y, c.z - 1);
        undoPos.push(cTmp);
        undoVal.push(mData[index(cTmp)]);
        todoPos.push(cTmp);
        mData[index(cTmp)] = -0.5;
      }
      if (c.x < (*this).getSizeX() - 1 && mData[index(c.x + 1, c.y, c.z)] > 0.) {
        cTmp = Vec3i(c.x + 1, c.y, c.z);
        undoPos.push(cTmp);
        undoVal.push(mData[index(cTmp)]);
        todoPos.push(cTmp);
        mData[index(cTmp)] = -0.5;
      }
      if (c.y < (*this).getSizeY() - 1 && mData[index(c.x, c.y + 1, c.z)] > 0.) {
        cTmp = Vec3i(c.x, c.y + 1, c.z);
        undoPos.push(cTmp);
        undoVal.push(mData[index(cTmp)]);
        todoPos.push(cTmp);
        mData[index(cTmp)] = -0.5;
      }
      if (c.z < (*this).getSizeZ() - 1 && mData[index(c.x, c.y, c.z + 1)] > 0.) {
        cTmp = Vec3i(c.x, c.y, c.z + 1);
        undoPos.push(cTmp);
        undoVal.push(mData[index(cTmp)]);
        todoPos.push(cTmp);
        mData[index(cTmp)] = -0.5;
      }

      /* Restore original value in cells if undo needed ie once cell undo count exceeds given limit
       */
      if (undoPos.size() > maxDepth) {
        /* Clear todo stack */
        while (!todoPos.empty()) {
          todoPos.pop();
        }
        /* Clear undo stack and revert value */
        while (!undoPos.empty()) {
          c = undoPos.top();
          curVal = undoVal.top();
          undoPos.pop();
          undoVal.pop();
          mData[index(c.x, c.y, c.z)] = curVal;
        }
        break;
      }

      /* Ensure that undo stack is cleared at the end if no more items in todo stack left */
      if (todoPos.empty()) {
        while (!undoPos.empty()) {
          undoPos.pop();
        }
        while (!undoVal.empty()) {
          undoVal.pop();
        }
      }
      /* Pop value for next while iteration */
      else {
        c = todoPos.top();
      }
    }
  }
}

//! run marching cubes to create a mesh for the 0-levelset
void LevelsetGrid::createMesh(Mesh &mesh)
{
  assertMsg(is3D(), "Only 3D grids supported so far");

  mesh.clear();

  const Real invalidTime = invalidTimeValue();
  const Real isoValue = 1e-4;

  // create some temp grids
  Grid<int> edgeVX(mParent);
  Grid<int> edgeVY(mParent);
  Grid<int> edgeVZ(mParent);

  for (int k = 0; k < mSize.z - 1; k++)
    for (int j = 0; j < mSize.y - 1; j++)
      for (int i = 0; i < mSize.x - 1; i++) {
        Real value[8] = {get(i, j, k),
                         get(i + 1, j, k),
                         get(i + 1, j + 1, k),
                         get(i, j + 1, k),
                         get(i, j, k + 1),
                         get(i + 1, j, k + 1),
                         get(i + 1, j + 1, k + 1),
                         get(i, j + 1, k + 1)};

        // build lookup index, check for invalid times
        bool skip = false;
        int cubeIdx = 0;
        for (int l = 0; l < 8; l++) {
          value[l] *= -1;
          if (-value[l] <= invalidTime)
            skip = true;
          if (value[l] < isoValue)
            cubeIdx |= 1 << l;
        }
        if (skip || (mcEdgeTable[cubeIdx] == 0))
          continue;

        // where to look up if this point already exists
        int triIndices[12];
        int *eVert[12] = {&edgeVX(i, j, k),
                          &edgeVY(i + 1, j, k),
                          &edgeVX(i, j + 1, k),
                          &edgeVY(i, j, k),
                          &edgeVX(i, j, k + 1),
                          &edgeVY(i + 1, j, k + 1),
                          &edgeVX(i, j + 1, k + 1),
                          &edgeVY(i, j, k + 1),
                          &edgeVZ(i, j, k),
                          &edgeVZ(i + 1, j, k),
                          &edgeVZ(i + 1, j + 1, k),
                          &edgeVZ(i, j + 1, k)};

        const Vec3 pos[9] = {Vec3(i, j, k),
                             Vec3(i + 1, j, k),
                             Vec3(i + 1, j + 1, k),
                             Vec3(i, j + 1, k),
                             Vec3(i, j, k + 1),
                             Vec3(i + 1, j, k + 1),
                             Vec3(i + 1, j + 1, k + 1),
                             Vec3(i, j + 1, k + 1)};

        for (int e = 0; e < 12; e++) {
          if (mcEdgeTable[cubeIdx] & (1 << e)) {
            // vertex already calculated ?
            if (*eVert[e] == 0) {
              // interpolate edge
              const int e1 = mcEdges[e * 2];
              const int e2 = mcEdges[e * 2 + 1];
              const Vec3 p1 = pos[e1];        // scalar field pos 1
              const Vec3 p2 = pos[e2];        // scalar field pos 2
              const float valp1 = value[e1];  // scalar field val 1
              const float valp2 = value[e2];  // scalar field val 2
              const float mu = (isoValue - valp1) / (valp2 - valp1);

              // init isolevel vertex
              Node vertex;
              vertex.pos = p1 + (p2 - p1) * mu + Vec3(Real(0.5));
              vertex.normal = getNormalized(
                  getGradient(
                      *this, i + cubieOffsetX[e1], j + cubieOffsetY[e1], k + cubieOffsetZ[e1]) *
                      (1.0 - mu) +
                  getGradient(
                      *this, i + cubieOffsetX[e2], j + cubieOffsetY[e2], k + cubieOffsetZ[e2]) *
                      (mu));

              triIndices[e] = mesh.addNode(vertex) + 1;

              // store vertex
              *eVert[e] = triIndices[e];
            }
            else {
              // retrieve  from vert array
              triIndices[e] = *eVert[e];
            }
          }
        }

        // Create the triangles...
        for (int e = 0; mcTriTable[cubeIdx][e] != -1; e += 3) {
          mesh.addTri(Triangle(triIndices[mcTriTable[cubeIdx][e + 0]] - 1,
                               triIndices[mcTriTable[cubeIdx][e + 1]] - 1,
                               triIndices[mcTriTable[cubeIdx][e + 2]] - 1));
        }
      }

  // mesh.rebuildCorners();
  // mesh.rebuildLookup();

  // Update mdata fields
  mesh.updateDataFields();
}

}  // namespace Manta