// 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
 *
 * Shape classes
 *
 ******************************************************************************/

#include "shapes.h"
#include "commonkernels.h"
#include "mesh.h"

using namespace std;
namespace Manta {

//******************************************************************************
// Shape class members

Shape::Shape(FluidSolver *parent) : PbClass(parent), mType(TypeNone)
{
}

LevelsetGrid Shape::computeLevelset()
{
  // note - 3d check deactivated! TODO double check...
  LevelsetGrid phi(getParent());
  generateLevelset(phi);
  return phi;
}

bool Shape::isInside(const Vec3 &pos) const
{
  return false;
}

//! Kernel: Apply a shape to a grid, setting value inside

template<class T> struct ApplyShapeToGrid : public KernelBase {
  ApplyShapeToGrid(Grid<T> *grid, Shape *shape, T value, FlagGrid *respectFlags)
      : KernelBase(grid, 0), grid(grid), shape(shape), value(value), respectFlags(respectFlags)
  {
    runMessage();
    run();
  }
  inline void op(
      int i, int j, int k, Grid<T> *grid, Shape *shape, T value, FlagGrid *respectFlags) const
  {
    if (respectFlags && respectFlags->isObstacle(i, j, k))
      return;
    if (shape->isInsideGrid(i, j, k))
      (*grid)(i, j, k) = value;
  }
  inline Grid<T> *getArg0()
  {
    return grid;
  }
  typedef Grid<T> type0;
  inline Shape *getArg1()
  {
    return shape;
  }
  typedef Shape type1;
  inline T &getArg2()
  {
    return value;
  }
  typedef T type2;
  inline FlagGrid *getArg3()
  {
    return respectFlags;
  }
  typedef FlagGrid type3;
  void runMessage()
  {
    debMsg("Executing kernel ApplyShapeToGrid ", 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, grid, shape, value, respectFlags);
    }
    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, grid, shape, value, respectFlags);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  Grid<T> *grid;
  Shape *shape;
  T value;
  FlagGrid *respectFlags;
};

//! Kernel: Apply a shape to a grid, setting value inside (scaling by SDF value)

template<class T> struct ApplyShapeToGridSmooth : public KernelBase {
  ApplyShapeToGridSmooth(
      Grid<T> *grid, Grid<Real> &phi, Real sigma, Real shift, T value, FlagGrid *respectFlags)
      : KernelBase(grid, 0),
        grid(grid),
        phi(phi),
        sigma(sigma),
        shift(shift),
        value(value),
        respectFlags(respectFlags)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 Grid<T> *grid,
                 Grid<Real> &phi,
                 Real sigma,
                 Real shift,
                 T value,
                 FlagGrid *respectFlags) const
  {
    if (respectFlags && respectFlags->isObstacle(i, j, k))
      return;
    const Real p = phi(i, j, k) - shift;
    if (p < -sigma)
      (*grid)(i, j, k) = value;
    else if (p < sigma)
      (*grid)(i, j, k) = value * (0.5f * (1.0f - p / sigma));
  }
  inline Grid<T> *getArg0()
  {
    return grid;
  }
  typedef Grid<T> type0;
  inline Grid<Real> &getArg1()
  {
    return phi;
  }
  typedef Grid<Real> type1;
  inline Real &getArg2()
  {
    return sigma;
  }
  typedef Real type2;
  inline Real &getArg3()
  {
    return shift;
  }
  typedef Real type3;
  inline T &getArg4()
  {
    return value;
  }
  typedef T type4;
  inline FlagGrid *getArg5()
  {
    return respectFlags;
  }
  typedef FlagGrid type5;
  void runMessage()
  {
    debMsg("Executing kernel ApplyShapeToGridSmooth ", 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, grid, phi, sigma, shift, value, respectFlags);
    }
    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, grid, phi, sigma, shift, value, respectFlags);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  Grid<T> *grid;
  Grid<Real> &phi;
  Real sigma;
  Real shift;
  T value;
  FlagGrid *respectFlags;
};

//! Kernel: Apply a shape to a MAC grid, setting value inside

struct ApplyShapeToMACGrid : public KernelBase {
  ApplyShapeToMACGrid(MACGrid *grid, Shape *shape, Vec3 value, FlagGrid *respectFlags)
      : KernelBase(grid, 0), grid(grid), shape(shape), value(value), respectFlags(respectFlags)
  {
    runMessage();
    run();
  }
  inline void op(
      int i, int j, int k, MACGrid *grid, Shape *shape, Vec3 value, FlagGrid *respectFlags) const
  {
    if (respectFlags && respectFlags->isObstacle(i, j, k))
      return;
    if (shape->isInside(Vec3(i, j + 0.5, k + 0.5)))
      (*grid)(i, j, k).x = value.x;
    if (shape->isInside(Vec3(i + 0.5, j, k + 0.5)))
      (*grid)(i, j, k).y = value.y;
    if (shape->isInside(Vec3(i + 0.5, j + 0.5, k)))
      (*grid)(i, j, k).z = value.z;
  }
  inline MACGrid *getArg0()
  {
    return grid;
  }
  typedef MACGrid type0;
  inline Shape *getArg1()
  {
    return shape;
  }
  typedef Shape type1;
  inline Vec3 &getArg2()
  {
    return value;
  }
  typedef Vec3 type2;
  inline FlagGrid *getArg3()
  {
    return respectFlags;
  }
  typedef FlagGrid type3;
  void runMessage()
  {
    debMsg("Executing kernel ApplyShapeToMACGrid ", 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, grid, shape, value, respectFlags);
    }
    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, grid, shape, value, respectFlags);
    }
  }
  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 *grid;
  Shape *shape;
  Vec3 value;
  FlagGrid *respectFlags;
};

void Shape::applyToGrid(GridBase *grid, FlagGrid *respectFlags)
{
#if NOPYTHON != 1
  if (grid->getType() & GridBase::TypeInt)
    ApplyShapeToGrid<int>((Grid<int> *)grid, this, _args.get<int>("value"), respectFlags);
  else if (grid->getType() & GridBase::TypeReal)
    ApplyShapeToGrid<Real>((Grid<Real> *)grid, this, _args.get<Real>("value"), respectFlags);
  else if (grid->getType() & GridBase::TypeMAC)
    ApplyShapeToMACGrid((MACGrid *)grid, this, _args.get<Vec3>("value"), respectFlags);
  else if (grid->getType() & GridBase::TypeVec3)
    ApplyShapeToGrid<Vec3>((Grid<Vec3> *)grid, this, _args.get<Vec3>("value"), respectFlags);
  else
    errMsg("Shape::applyToGrid(): unknown grid type");
#else
  errMsg("Not yet supported...");
#endif
}

void Shape::applyToGridSmooth(GridBase *grid, Real sigma, Real shift, FlagGrid *respectFlags)
{
  Grid<Real> phi(grid->getParent());
  generateLevelset(phi);

#if NOPYTHON != 1
  if (grid->getType() & GridBase::TypeInt)
    ApplyShapeToGridSmooth<int>(
        (Grid<int> *)grid, phi, sigma, shift, _args.get<int>("value"), respectFlags);
  else if (grid->getType() & GridBase::TypeReal)
    ApplyShapeToGridSmooth<Real>(
        (Grid<Real> *)grid, phi, sigma, shift, _args.get<Real>("value"), respectFlags);
  else if (grid->getType() & GridBase::TypeVec3)
    ApplyShapeToGridSmooth<Vec3>(
        (Grid<Vec3> *)grid, phi, sigma, shift, _args.get<Vec3>("value"), respectFlags);
  else
    errMsg("Shape::applyToGridSmooth(): unknown grid type");
#else
  errMsg("Not yet supported...");
#endif
}

void Shape::collideMesh(Mesh &mesh)
{
  const Real margin = 0.2;

  Grid<Real> phi(getParent());
  Grid<Vec3> grad(getParent());
  generateLevelset(phi);
  GradientOp(grad, phi);

  const int num = mesh.numNodes();
  for (int i = 0; i < num; i++) {
    const Vec3 &p = mesh.nodes(i).pos;
    mesh.nodes(i).flags &= ~(Mesh::NfCollide | Mesh::NfMarked);
    if (!phi.isInBounds(p, 1))
      continue;

    for (int iter = 0; iter < 10; iter++) {
      const Real dist = phi.getInterpolated(p);
      if (dist < margin) {
        Vec3 n = grad.getInterpolated(p);
        normalize(n);
        mesh.nodes(i).pos += (margin - dist) * n;
        mesh.nodes(i).flags |= Mesh::NfCollide | Mesh::NfMarked;
      }
      else
        break;
    }
  }
}

//******************************************************************************
// Derived shape class members

Box::Box(FluidSolver *parent, Vec3 center, Vec3 p0, Vec3 p1, Vec3 size) : Shape(parent)
{
  mType = TypeBox;
  if (center.isValid() && size.isValid()) {
    mP0 = center - size;
    mP1 = center + size;
  }
  else if (p0.isValid() && p1.isValid()) {
    mP0 = p0;
    mP1 = p1;
  }
  else
    errMsg("Box: specify either p0,p1 or size,center");
}

bool Box::isInside(const Vec3 &pos) const
{
  return (pos.x >= mP0.x && pos.y >= mP0.y && pos.z >= mP0.z && pos.x <= mP1.x && pos.y <= mP1.y &&
          pos.z <= mP1.z);
}

void Box::generateMesh(Mesh *mesh)
{
  const int quadidx[24] = {0, 4, 6, 2, 3, 7, 5, 1, 0, 1, 5, 4, 6, 7, 3, 2, 0, 2, 3, 1, 5, 7, 6, 4};
  const int nodebase = mesh->numNodes();
  int oldtri = mesh->numTris();
  for (int i = 0; i < 8; i++) {
    Node p;
    p.flags = 0;
    p.pos = mP0;
    if (i & 1)
      p.pos.x = mP1.x;
    if (i & 2)
      p.pos.y = mP1.y;
    if (i & 4)
      p.pos.z = mP1.z;
    mesh->addNode(p);
  }
  for (int i = 0; i < 6; i++) {
    mesh->addTri(Triangle(nodebase + quadidx[i * 4 + 0],
                          nodebase + quadidx[i * 4 + 1],
                          nodebase + quadidx[i * 4 + 3]));
    mesh->addTri(Triangle(nodebase + quadidx[i * 4 + 1],
                          nodebase + quadidx[i * 4 + 2],
                          nodebase + quadidx[i * 4 + 3]));
  }
  mesh->rebuildCorners(oldtri, -1);
  mesh->rebuildLookup(oldtri, -1);
}

//! Kernel: Analytic SDF for box shape
struct BoxSDF : public KernelBase {
  BoxSDF(Grid<Real> &phi, const Vec3 &p1, const Vec3 &p2)
      : KernelBase(&phi, 0), phi(phi), p1(p1), p2(p2)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, Grid<Real> &phi, const Vec3 &p1, const Vec3 &p2) const
  {
    const Vec3 p(i + 0.5, j + 0.5, k + 0.5);
    if (p.x <= p2.x && p.x >= p1.x && p.y <= p2.y && p.y >= p1.y && p.z <= p2.z && p.z >= p1.z) {
      // inside: minimal surface distance
      Real mx = max(p.x - p2.x, p1.x - p.x);
      Real my = max(p.y - p2.y, p1.y - p.y);
      Real mz = max(p.z - p2.z, p1.z - p.z);
      if (!phi.is3D())
        mz = mx;  // skip for 2d...
      phi(i, j, k) = max(mx, max(my, mz));
    }
    else if (p.y <= p2.y && p.y >= p1.y && p.z <= p2.z && p.z >= p1.z) {
      // outside plane X
      phi(i, j, k) = max(p.x - p2.x, p1.x - p.x);
    }
    else if (p.x <= p2.x && p.x >= p1.x && p.z <= p2.z && p.z >= p1.z) {
      // outside plane Y
      phi(i, j, k) = max(p.y - p2.y, p1.y - p.y);
    }
    else if (p.x <= p2.x && p.x >= p1.x && p.y <= p2.y && p.y >= p1.y) {
      // outside plane Z
      phi(i, j, k) = max(p.z - p2.z, p1.z - p.z);
    }
    else if (p.x > p1.x && p.x < p2.x) {
      // lines X
      Real m1 = sqrt(square(p1.y - p.y) + square(p1.z - p.z));
      Real m2 = sqrt(square(p2.y - p.y) + square(p1.z - p.z));
      Real m3 = sqrt(square(p1.y - p.y) + square(p2.z - p.z));
      Real m4 = sqrt(square(p2.y - p.y) + square(p2.z - p.z));
      phi(i, j, k) = min(m1, min(m2, min(m3, m4)));
    }
    else if (p.y > p1.y && p.y < p2.y) {
      // lines Y
      Real m1 = sqrt(square(p1.x - p.x) + square(p1.z - p.z));
      Real m2 = sqrt(square(p2.x - p.x) + square(p1.z - p.z));
      Real m3 = sqrt(square(p1.x - p.x) + square(p2.z - p.z));
      Real m4 = sqrt(square(p2.x - p.x) + square(p2.z - p.z));
      phi(i, j, k) = min(m1, min(m2, min(m3, m4)));
    }
    else if (p.z > p1.x && p.z < p2.z) {
      // lines Z
      Real m1 = sqrt(square(p1.y - p.y) + square(p1.x - p.x));
      Real m2 = sqrt(square(p2.y - p.y) + square(p1.x - p.x));
      Real m3 = sqrt(square(p1.y - p.y) + square(p2.x - p.x));
      Real m4 = sqrt(square(p2.y - p.y) + square(p2.x - p.x));
      phi(i, j, k) = min(m1, min(m2, min(m3, m4)));
    }
    else {
      // points
      Real m = norm(p - Vec3(p1.x, p1.y, p1.z));
      m = min(m, norm(p - Vec3(p1.x, p1.y, p2.z)));
      m = min(m, norm(p - Vec3(p1.x, p2.y, p1.z)));
      m = min(m, norm(p - Vec3(p1.x, p2.y, p2.z)));
      m = min(m, norm(p - Vec3(p2.x, p1.y, p1.z)));
      m = min(m, norm(p - Vec3(p2.x, p1.y, p2.z)));
      m = min(m, norm(p - Vec3(p2.x, p2.y, p1.z)));
      m = min(m, norm(p - Vec3(p2.x, p2.y, p2.z)));
      phi(i, j, k) = m;
    }
  }
  inline Grid<Real> &getArg0()
  {
    return phi;
  }
  typedef Grid<Real> type0;
  inline const Vec3 &getArg1()
  {
    return p1;
  }
  typedef Vec3 type1;
  inline const Vec3 &getArg2()
  {
    return p2;
  }
  typedef Vec3 type2;
  void runMessage()
  {
    debMsg("Executing kernel BoxSDF ", 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, phi, p1, p2);
    }
    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, phi, p1, p2);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  Grid<Real> &phi;
  const Vec3 &p1;
  const Vec3 &p2;
};
void Box::generateLevelset(Grid<Real> &phi)
{
  BoxSDF(phi, mP0, mP1);
}

Sphere::Sphere(FluidSolver *parent, Vec3 center, Real radius, Vec3 scale)
    : Shape(parent), mCenter(center), mScale(scale), mRadius(radius)
{
  mType = TypeSphere;
}

bool Sphere::isInside(const Vec3 &pos) const
{
  return normSquare((pos - mCenter) / mScale) <= mRadius * mRadius;
}

struct Tri {
  Vec3 t[3];
  int i[3];
  Tri(Vec3 a, Vec3 b, Vec3 c)
  {
    t[0] = a;
    t[1] = b;
    t[2] = c;
  }
};
void Sphere::generateMesh(Mesh *mesh)
{
  vector<Tri> tris;
  const int iterations = 3;
  int oldtri = mesh->numTris();

  // start with octahedron
  const Real d = sqrt(0.5);
  Vec3 p[6] = {Vec3(0, 1, 0),
               Vec3(0, -1, 0),
               Vec3(-d, 0, -d),
               Vec3(d, 0, -d),
               Vec3(d, 0, d),
               Vec3(-d, 0, d)};
  tris.push_back(Tri(p[0], p[4], p[3]));
  tris.push_back(Tri(p[0], p[5], p[4]));
  tris.push_back(Tri(p[0], p[2], p[5]));
  tris.push_back(Tri(p[0], p[3], p[2]));
  tris.push_back(Tri(p[1], p[3], p[4]));
  tris.push_back(Tri(p[1], p[4], p[5]));
  tris.push_back(Tri(p[1], p[5], p[2]));
  tris.push_back(Tri(p[1], p[2], p[3]));

  // Bisect each edge and move to the surface of a unit sphere
  for (int it = 0; it < iterations; it++) {
    int ntold = tris.size();
    for (int i = 0; i < ntold; i++) {
      Vec3 pa = 0.5 * (tris[i].t[0] + tris[i].t[1]);
      Vec3 pb = 0.5 * (tris[i].t[1] + tris[i].t[2]);
      Vec3 pc = 0.5 * (tris[i].t[2] + tris[i].t[0]);
      normalize(pa);
      normalize(pb);
      normalize(pc);

      tris.push_back(Tri(tris[i].t[0], pa, pc));
      tris.push_back(Tri(pa, tris[i].t[1], pb));
      tris.push_back(Tri(pb, tris[i].t[2], pc));
      tris[i].t[0] = pa;
      tris[i].t[1] = pb;
      tris[i].t[2] = pc;
    }
  }

  // index + scale
  vector<Vec3> nodes;
  for (size_t i = 0; i < tris.size(); i++) {
    for (int t = 0; t < 3; t++) {
      Vec3 p = mCenter + tris[i].t[t] * mRadius * mScale;
      // vector already there ?
      int idx = nodes.size();
      for (size_t j = 0; j < nodes.size(); j++) {
        if (p == nodes[j]) {
          idx = j;
          break;
        }
      }
      if (idx == (int)nodes.size())
        nodes.push_back(p);
      tris[i].i[t] = idx;
    }
  }

  // add the to mesh
  const int ni = mesh->numNodes();
  for (size_t i = 0; i < nodes.size(); i++) {
    mesh->addNode(Node(nodes[i]));
  }
  for (size_t t = 0; t < tris.size(); t++)
    mesh->addTri(Triangle(tris[t].i[0] + ni, tris[t].i[1] + ni, tris[t].i[2] + ni));

  mesh->rebuildCorners(oldtri, -1);
  mesh->rebuildLookup(oldtri, -1);
}

struct SphereSDF : public KernelBase {
  SphereSDF(Grid<Real> &phi, Vec3 center, Real radius, Vec3 scale)
      : KernelBase(&phi, 0), phi(phi), center(center), radius(radius), scale(scale)
  {
    runMessage();
    run();
  }
  inline void op(int i, int j, int k, Grid<Real> &phi, Vec3 center, Real radius, Vec3 scale) const
  {
    phi(i, j, k) = norm((Vec3(i + 0.5, j + 0.5, k + 0.5) - center) / scale) - radius;
  }
  inline Grid<Real> &getArg0()
  {
    return phi;
  }
  typedef Grid<Real> type0;
  inline Vec3 &getArg1()
  {
    return center;
  }
  typedef Vec3 type1;
  inline Real &getArg2()
  {
    return radius;
  }
  typedef Real type2;
  inline Vec3 &getArg3()
  {
    return scale;
  }
  typedef Vec3 type3;
  void runMessage()
  {
    debMsg("Executing kernel SphereSDF ", 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, phi, center, radius, scale);
    }
    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, phi, center, radius, scale);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  Grid<Real> &phi;
  Vec3 center;
  Real radius;
  Vec3 scale;
};
void Sphere::generateLevelset(Grid<Real> &phi)
{
  SphereSDF(phi, mCenter, mRadius, mScale);
}

Cylinder::Cylinder(FluidSolver *parent, Vec3 center, Real radius, Vec3 z)
    : Shape(parent), mCenter(center), mRadius(radius)
{
  mType = TypeCylinder;
  mZDir = z;
  mZ = normalize(mZDir);
}

bool Cylinder::isInside(const Vec3 &pos) const
{
  Real z = dot(pos - mCenter, mZDir);
  if (fabs(z) > mZ)
    return false;
  Real r2 = normSquare(pos - mCenter) - square(z);
  return r2 < square(mRadius);
}

void Cylinder::generateMesh(Mesh *mesh)
{
  // generate coordinate system
  Vec3 x = getOrthogonalVector(mZDir) * mRadius;
  Vec3 y = cross(x, mZDir);
  Vec3 z = mZDir * mZ;
  int oldtri = mesh->numTris();

  // construct node ring
  const int N = 20;
  const int base = mesh->numNodes();
  for (int i = 0; i < N; i++) {
    const Real phi = 2.0 * M_PI * (Real)i / (Real)N;
    Vec3 r = x * cos(phi) + y * sin(phi) + mCenter;
    mesh->addNode(Node(r + z));
    mesh->addNode(Node(r - z));
  }
  // top/bottom center
  mesh->addNode(Node(mCenter + z));
  mesh->addNode(Node(mCenter - z));

  // connect with tris
  for (int i = 0; i < N; i++) {
    int cur = base + 2 * i;
    int next = base + 2 * ((i + 1) % N);
    // outside
    mesh->addTri(Triangle(cur, next, cur + 1));
    mesh->addTri(Triangle(next, next + 1, cur + 1));
    // upper / lower
    mesh->addTri(Triangle(cur, base + 2 * N, next));
    mesh->addTri(Triangle(cur + 1, next + 1, base + 2 * N + 1));
  }

  mesh->rebuildCorners(oldtri, -1);
  mesh->rebuildLookup(oldtri, -1);
}

struct CylinderSDF : public KernelBase {
  CylinderSDF(Grid<Real> &phi, Vec3 center, Real radius, Vec3 zaxis, Real maxz)
      : KernelBase(&phi, 0), phi(phi), center(center), radius(radius), zaxis(zaxis), maxz(maxz)
  {
    runMessage();
    run();
  }
  inline void op(
      int i, int j, int k, Grid<Real> &phi, Vec3 center, Real radius, Vec3 zaxis, Real maxz) const
  {
    Vec3 p = Vec3(i + 0.5, j + 0.5, k + 0.5) - center;
    Real z = fabs(dot(p, zaxis));
    Real r = sqrt(normSquare(p) - z * z);
    if (z < maxz) {
      // cylinder z area
      if (r < radius)
        phi(i, j, k) = max(r - radius, z - maxz);
      else
        phi(i, j, k) = r - radius;
    }
    else if (r < radius) {
      // cylinder top area
      phi(i, j, k) = fabs(z - maxz);
    }
    else {
      // edge
      phi(i, j, k) = sqrt(square(z - maxz) + square(r - radius));
    }
  }
  inline Grid<Real> &getArg0()
  {
    return phi;
  }
  typedef Grid<Real> type0;
  inline Vec3 &getArg1()
  {
    return center;
  }
  typedef Vec3 type1;
  inline Real &getArg2()
  {
    return radius;
  }
  typedef Real type2;
  inline Vec3 &getArg3()
  {
    return zaxis;
  }
  typedef Vec3 type3;
  inline Real &getArg4()
  {
    return maxz;
  }
  typedef Real type4;
  void runMessage()
  {
    debMsg("Executing kernel CylinderSDF ", 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, phi, center, radius, zaxis, maxz);
    }
    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, phi, center, radius, zaxis, maxz);
    }
  }
  void run()
  {
    if (maxZ > 1)
      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
    else
      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
  }
  Grid<Real> &phi;
  Vec3 center;
  Real radius;
  Vec3 zaxis;
  Real maxz;
};
void Cylinder::generateLevelset(Grid<Real> &phi)
{
  CylinderSDF(phi, mCenter, mRadius, mZDir, mZ);
}

Slope::Slope(FluidSolver *parent, Real anglexy, Real angleyz, Real origin, Vec3 gs)
    : Shape(parent), mAnglexy(anglexy), mAngleyz(angleyz), mOrigin(origin), mGs(gs)
{
  mType = TypeSlope;
}

void Slope::generateMesh(Mesh *mesh)
{

  const int oldtri = mesh->numTris();

  Vec3 v1(0., mOrigin, 0.);
  mesh->addNode(Node(v1));

  Real dy1 = mGs.z * std::tan(mAngleyz);
  Vec3 v2(0., mOrigin - dy1, mGs.z);
  mesh->addNode(Node(v2));

  Real dy2 = mGs.x * std::tan(mAnglexy);
  Vec3 v3(mGs.x, v2.y - dy2, mGs.z);
  mesh->addNode(Node(v3));

  Vec3 v4(mGs.x, mOrigin - dy2, 0.);
  mesh->addNode(Node(v4));

  mesh->addTri(Triangle(0, 1, 2));
  mesh->addTri(Triangle(2, 3, 0));

  mesh->rebuildCorners(oldtri, -1);
  mesh->rebuildLookup(oldtri, -1);
}

bool Slope::isInside(const Vec3 &pos) const
{

  const Real alpha = -mAnglexy * M_PI / 180.;
  const Real beta = -mAngleyz * M_PI / 180.;

  Vec3 n(0, 1, 0);

  n.x = std::sin(alpha) * std::cos(beta);
  n.y = std::cos(alpha) * std::cos(beta);
  n.z = std::sin(beta);

  normalize(n);

  const Real fac = std::sqrt(n.x * n.x + n.y * n.y + n.z * n.z);

  return ((n.x * (double)pos.x + n.y * (double)pos.y + n.z * (double)pos.z - mOrigin) / fac) <= 0.;
}

struct SlopeSDF : public KernelBase {
  SlopeSDF(const Vec3 &n, Grid<Real> &phiObs, const Real &fac, const Real &origin)
      : KernelBase(&phiObs, 0), n(n), phiObs(phiObs), fac(fac), origin(origin)
  {
    runMessage();
    run();
  }
  inline void op(int i,
                 int j,
                 int k,
                 const Vec3 &n,
                 Grid<Real> &phiObs,
                 const Real &fac,
                 const Real &origin) const
  {

    phiObs(i, j, k) = (n.x * (double)i + n.y * (double)j + n.z * (double)k - origin) * fac;
  }
  inline const Vec3 &getArg0()
  {
    return n;
  }
  typedef Vec3 type0;
  inline Grid<Real> &getArg1()
  {
    return phiObs;
  }
  typedef Grid<Real> type1;
  inline const Real &getArg2()
  {
    return fac;
  }
  typedef Real type2;
  inline const Real &getArg3()
  {
    return origin;
  }
  typedef Real type3;
  void runMessage()
  {
    debMsg("Executing kernel SlopeSDF ", 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, n, phiObs, fac, origin);
    }
    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, n, phiObs, fac, origin);
    }
  }
  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 Vec3 &n;
  Grid<Real> &phiObs;
  const Real &fac;
  const Real &origin;
};

void Slope::generateLevelset(Grid<Real> &phi)
{

  const Real alpha = -mAnglexy * M_PI / 180.;
  const Real beta = -mAngleyz * M_PI / 180.;

  Vec3 n(0, 1, 0);

  n.x = std::sin(alpha) * std::cos(beta);
  n.y = std::cos(alpha) * std::cos(beta);
  n.z = std::sin(beta);

  normalize(n);

  const Real fac = 1. / std::sqrt(n.x * n.x + n.y * n.y + n.z * n.z);

  SlopeSDF(n, phi, fac, mOrigin);
}

}  // namespace Manta