// 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
 *
 * Meshes
 *
 * note: this is only a temporary solution, details are bound to change
 *        long term goal is integration with Split&Merge code by Wojtan et al.
 *
 ******************************************************************************/

#include "mesh.h"
#include "integrator.h"
#include "mantaio.h"
#include "kernel.h"
#include "shapes.h"
#include "noisefield.h"
//#include "grid.h"
#include <stack>
#include <cstring>

using namespace std;
namespace Manta {

Mesh::Mesh(FluidSolver *parent) : PbClass(parent)
{
}

Mesh::~Mesh()
{
  for (IndexInt i = 0; i < (IndexInt)mMeshData.size(); ++i)
    mMeshData[i]->setMesh(NULL);

  if (mFreeMdata) {
    for (IndexInt i = 0; i < (IndexInt)mMeshData.size(); ++i)
      delete mMeshData[i];
  }
}

Mesh *Mesh::clone()
{
  Mesh *nm = new Mesh(mParent);
  *nm = *this;
  nm->setName(getName());
  return nm;
}

void Mesh::deregister(MeshDataBase *mdata)
{
  bool done = false;
  // remove pointer from mesh data list
  for (IndexInt i = 0; i < (IndexInt)mMeshData.size(); ++i) {
    if (mMeshData[i] == mdata) {
      if (i < (IndexInt)mMeshData.size() - 1)
        mMeshData[i] = mMeshData[mMeshData.size() - 1];
      mMeshData.pop_back();
      done = true;
    }
  }
  if (!done)
    errMsg("Invalid pointer given, not registered!");
}

// create and attach a new mdata field to this mesh
PbClass *Mesh::create(PbType t, PbTypeVec T, const string &name)
{
#if NOPYTHON != 1
  _args.add("nocheck", true);
  if (t.str() == "")
    errMsg("Specify mesh data type to create");
  // debMsg( "Mdata creating '"<< t.str <<" with size "<< this->getSizeSlow(), 5 );

  PbClass *pyObj = PbClass::createPyObject(t.str() + T.str(), name, _args, this->getParent());

  MeshDataBase *mdata = dynamic_cast<MeshDataBase *>(pyObj);
  if (!mdata) {
    errMsg(
        "Unable to get mesh data pointer from newly created object. Only create MeshData type "
        "with a Mesh.creat() call, eg, MdataReal, MdataVec3 etc.");
    delete pyObj;
    return NULL;
  }
  else {
    this->registerMdata(mdata);
  }

  // directly init size of new mdata field:
  mdata->resize(this->getSizeSlow());
#else
  PbClass *pyObj = NULL;
#endif
  return pyObj;
}

void Mesh::registerMdata(MeshDataBase *mdata)
{
  mdata->setMesh(this);
  mMeshData.push_back(mdata);

  if (mdata->getType() == MeshDataBase::TypeReal) {
    MeshDataImpl<Real> *pd = dynamic_cast<MeshDataImpl<Real> *>(mdata);
    if (!pd)
      errMsg("Invalid mdata object posing as real!");
    this->registerMdataReal(pd);
  }
  else if (mdata->getType() == MeshDataBase::TypeInt) {
    MeshDataImpl<int> *pd = dynamic_cast<MeshDataImpl<int> *>(mdata);
    if (!pd)
      errMsg("Invalid mdata object posing as int!");
    this->registerMdataInt(pd);
  }
  else if (mdata->getType() == MeshDataBase::TypeVec3) {
    MeshDataImpl<Vec3> *pd = dynamic_cast<MeshDataImpl<Vec3> *>(mdata);
    if (!pd)
      errMsg("Invalid mdata object posing as vec3!");
    this->registerMdataVec3(pd);
  }
}
void Mesh::registerMdataReal(MeshDataImpl<Real> *pd)
{
  mMdataReal.push_back(pd);
}
void Mesh::registerMdataVec3(MeshDataImpl<Vec3> *pd)
{
  mMdataVec3.push_back(pd);
}
void Mesh::registerMdataInt(MeshDataImpl<int> *pd)
{
  mMdataInt.push_back(pd);
}

void Mesh::addAllMdata()
{
  for (IndexInt i = 0; i < (IndexInt)mMeshData.size(); ++i) {
    mMeshData[i]->addEntry();
  }
}

Real Mesh::computeCenterOfMass(Vec3 &cm) const
{

  // use double precision for summation, otherwise too much error accumulation
  double vol = 0;
  Vector3D<double> cmd(0.0);
  for (size_t tri = 0; tri < mTris.size(); tri++) {
    Vector3D<double> p1(toVec3d(getNode(tri, 0)));
    Vector3D<double> p2(toVec3d(getNode(tri, 1)));
    Vector3D<double> p3(toVec3d(getNode(tri, 2)));

    double cvol = dot(cross(p1, p2), p3) / 6.0;
    cmd += (p1 + p2 + p3) * (cvol / 4.0);
    vol += cvol;
  }
  if (vol != 0.0)
    cmd /= vol;

  cm = toVec3(cmd);
  return (Real)vol;
}

void Mesh::clear()
{
  mNodes.clear();
  mTris.clear();
  mCorners.clear();
  m1RingLookup.clear();
  for (size_t i = 0; i < mNodeChannels.size(); i++)
    mNodeChannels[i]->resize(0);
  for (size_t i = 0; i < mTriChannels.size(); i++)
    mTriChannels[i]->resize(0);

  // clear mdata fields as well
  for (size_t i = 0; i < mMdataReal.size(); i++)
    mMdataReal[i]->resize(0);
  for (size_t i = 0; i < mMdataVec3.size(); i++)
    mMdataVec3[i]->resize(0);
  for (size_t i = 0; i < mMdataInt.size(); i++)
    mMdataInt[i]->resize(0);
}

Mesh &Mesh::operator=(const Mesh &o)
{
  // wipe current data
  clear();
  if (mNodeChannels.size() != o.mNodeChannels.size() ||
      mTriChannels.size() != o.mTriChannels.size())
    errMsg("can't copy mesh, channels not identical");
  mNodeChannels.clear();
  mTriChannels.clear();

  // copy corner, nodes, tris
  mCorners = o.mCorners;
  mNodes = o.mNodes;
  mTris = o.mTris;
  m1RingLookup = o.m1RingLookup;

  // copy channels
  for (size_t i = 0; i < mNodeChannels.size(); i++)
    mNodeChannels[i] = o.mNodeChannels[i];
  for (size_t i = 0; i < o.mTriChannels.size(); i++)
    mTriChannels[i] = o.mTriChannels[i];

  return *this;
}

int Mesh::load(string name, bool append)
{
  if (name.find_last_of('.') == string::npos)
    errMsg("file '" + name + "' does not have an extension");
  string ext = name.substr(name.find_last_of('.'));
  if (ext == ".gz")  // assume bobj gz
    return readBobjFile(name, this, append);
  else if (ext == ".obj")
    return readObjFile(name, this, append);
  else
    errMsg("file '" + name + "' filetype not supported");

  // dont always rebuild...
  // rebuildCorners();
  // rebuildLookup();
  return 0;
}

int Mesh::save(string name)
{
  if (name.find_last_of('.') == string::npos)
    errMsg("file '" + name + "' does not have an extension");
  string ext = name.substr(name.find_last_of('.'));
  if (ext == ".obj")
    return writeObjFile(name, this);
  else if (ext == ".gz")
    return writeBobjFile(name, this);
  else
    errMsg("file '" + name + "' filetype not supported");
  return 0;
}

void Mesh::fromShape(Shape &shape, bool append)
{
  if (!append)
    clear();
  shape.generateMesh(this);
}

void Mesh::resizeTris(int numTris)
{
  mTris.resize(numTris);
  rebuildChannels();
}
void Mesh::resizeNodes(int numNodes)
{
  mNodes.resize(numNodes);
  rebuildChannels();
}

//! do a quick check whether a rebuild is necessary, and if yes do rebuild
void Mesh::rebuildQuickCheck()
{
  if (mCorners.size() != 3 * mTris.size())
    rebuildCorners();
  if (m1RingLookup.size() != mNodes.size())
    rebuildLookup();
}

void Mesh::rebuildCorners(int from, int to)
{
  mCorners.resize(3 * mTris.size());
  if (to < 0)
    to = mTris.size();

  // fill in basic info
  for (int tri = from; tri < to; tri++) {
    for (int c = 0; c < 3; c++) {
      const int idx = tri * 3 + c;
      mCorners[idx].tri = tri;
      mCorners[idx].node = mTris[tri].c[c];
      mCorners[idx].next = 3 * tri + ((c + 1) % 3);
      mCorners[idx].prev = 3 * tri + ((c + 2) % 3);
      mCorners[idx].opposite = -1;
    }
  }

  // set opposite info
  int maxc = to * 3;
  for (int c = from * 3; c < maxc; c++) {
    int next = mCorners[mCorners[c].next].node;
    int prev = mCorners[mCorners[c].prev].node;

    // find corner with same next/prev nodes
    for (int c2 = c + 1; c2 < maxc; c2++) {
      int next2 = mCorners[mCorners[c2].next].node;
      if (next2 != next && next2 != prev)
        continue;
      int prev2 = mCorners[mCorners[c2].prev].node;
      if (prev2 != next && prev2 != prev)
        continue;

      // found
      mCorners[c].opposite = c2;
      mCorners[c2].opposite = c;
      break;
    }
    if (mCorners[c].opposite < 0) {
      // didn't find opposite
      errMsg("can't rebuild corners, index without an opposite");
    }
  }

  rebuildChannels();
}

void Mesh::rebuildLookup(int from, int to)
{
  if (from == 0 && to < 0)
    m1RingLookup.clear();
  m1RingLookup.resize(mNodes.size());
  if (to < 0)
    to = mTris.size();
  from *= 3;
  to *= 3;
  for (int i = from; i < to; i++) {
    const int node = mCorners[i].node;
    m1RingLookup[node].nodes.insert(mCorners[mCorners[i].next].node);
    m1RingLookup[node].nodes.insert(mCorners[mCorners[i].prev].node);
    m1RingLookup[node].tris.insert(mCorners[i].tri);
  }
}

void Mesh::rebuildChannels()
{
  for (size_t i = 0; i < mTriChannels.size(); i++)
    mTriChannels[i]->resize(mTris.size());
  for (size_t i = 0; i < mNodeChannels.size(); i++)
    mNodeChannels[i]->resize(mNodes.size());
}

struct _KnAdvectMeshInGrid : public KernelBase {
  _KnAdvectMeshInGrid(const KernelBase &base,
                      vector<Node> &nodes,
                      const FlagGrid &flags,
                      const MACGrid &vel,
                      const Real dt,
                      vector<Vec3> &u)
      : KernelBase(base), nodes(nodes), flags(flags), vel(vel), dt(dt), u(u)
  {
  }
  inline void op(IndexInt idx,
                 vector<Node> &nodes,
                 const FlagGrid &flags,
                 const MACGrid &vel,
                 const Real dt,
                 vector<Vec3> &u) const
  {
    if (nodes[idx].flags & Mesh::NfFixed)
      u[idx] = 0.0;
    else if (!flags.isInBounds(nodes[idx].pos, 1))
      u[idx] = 0.0;
    else
      u[idx] = vel.getInterpolated(nodes[idx].pos) * dt;
  }
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, nodes, flags, vel, dt, u);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  vector<Node> &nodes;
  const FlagGrid &flags;
  const MACGrid &vel;
  const Real dt;
  vector<Vec3> &u;
};
struct KnAdvectMeshInGrid : public KernelBase {
  KnAdvectMeshInGrid(vector<Node> &nodes, const FlagGrid &flags, const MACGrid &vel, const Real dt)
      : KernelBase(nodes.size()),
        _inner(KernelBase(nodes.size()), nodes, flags, vel, dt, u),
        nodes(nodes),
        flags(flags),
        vel(vel),
        dt(dt),
        u((size))
  {
    runMessage();
    run();
  }
  void run()
  {
    _inner.run();
  }
  inline operator vector<Vec3>()
  {
    return u;
  }
  inline vector<Vec3> &getRet()
  {
    return u;
  }
  inline vector<Node> &getArg0()
  {
    return nodes;
  }
  typedef vector<Node> type0;
  inline const FlagGrid &getArg1()
  {
    return flags;
  }
  typedef FlagGrid type1;
  inline const MACGrid &getArg2()
  {
    return vel;
  }
  typedef MACGrid type2;
  inline const Real &getArg3()
  {
    return dt;
  }
  typedef Real type3;
  void runMessage()
  {
    debMsg("Executing kernel KnAdvectMeshInGrid ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  _KnAdvectMeshInGrid _inner;
  vector<Node> &nodes;
  const FlagGrid &flags;
  const MACGrid &vel;
  const Real dt;
  vector<Vec3> u;
};

// advection plugin
void Mesh::advectInGrid(FlagGrid &flags, MACGrid &vel, int integrationMode)
{
  KnAdvectMeshInGrid kernel(mNodes, flags, vel, getParent()->getDt());
  integratePointSet(kernel, integrationMode);
}

void Mesh::scale(Vec3 s)
{
  for (size_t i = 0; i < mNodes.size(); i++)
    mNodes[i].pos *= s;
}

void Mesh::offset(Vec3 o)
{
  for (size_t i = 0; i < mNodes.size(); i++)
    mNodes[i].pos += o;
}

void Mesh::rotate(Vec3 thetas)
{
  // rotation thetas are in radians (e.g. pi is equal to 180 degrees)
  auto rotate = [&](Real theta, unsigned int first_axis, unsigned int second_axis) {
    if (theta == 0.0f)
      return;

    Real sin_t = sin(theta);
    Real cos_t = cos(theta);

    Real sin_sign = first_axis == 0u && second_axis == 2u ? -1.0f : 1.0f;
    sin_t *= sin_sign;

    size_t length = mNodes.size();
    for (size_t n = 0; n < length; ++n) {
      Vec3 &node = mNodes[n].pos;
      Real first_axis_val = node[first_axis];
      Real second_axis_val = node[second_axis];
      node[first_axis] = first_axis_val * cos_t - second_axis_val * sin_t;
      node[second_axis] = second_axis_val * cos_t + first_axis_val * sin_t;
    }
  };

  // rotate x
  rotate(thetas[0], 1u, 2u);
  // rotate y
  rotate(thetas[1], 0u, 2u);
  // rotate z
  rotate(thetas[2], 0u, 1u);
}

void Mesh::computeVelocity(Mesh &oldMesh, MACGrid &vel)
{
  // Early return if sizes do not match
  if (oldMesh.mNodes.size() != mNodes.size())
    return;

  // temp grid
  Grid<Vec3> veloMeanCounter(getParent());
  veloMeanCounter.setConst(0.0f);

  bool bIs2D = getParent()->is2D();

  // calculate velocities from previous to current frame (per vertex)
  for (size_t i = 0; i < mNodes.size(); ++i) {
    // skip vertices that are not needed for 2D
    if (bIs2D && (mNodes[i].pos.z < -0.5f || mNodes[i].pos.z > 0.5f))
      continue;

    Vec3 velo = mNodes[i].pos - oldMesh.mNodes[i].pos;
    vel.setInterpolated(mNodes[i].pos, velo, &(veloMeanCounter[0]));
  }

  // discretize the vertex velocities by averaging them on the grid
  vel.safeDivide(veloMeanCounter);
}

void Mesh::removeTri(int tri)
{
  // delete triangles by overwriting them with elements from the end of the array.
  if (tri != (int)mTris.size() - 1) {
    // if this is the last element, and it is marked for deletion,
    // don't waste cycles transfering data to itself,
    // and DEFINITELY don't transfer .opposite data to other, untainted triangles.

    // old corners hold indices on the end of the corners array
    // new corners holds indices in the new spot in the middle of the array
    Corner *oldcorners[3];
    Corner *newcorners[3];
    int oldtri = mTris.size() - 1;
    for (int c = 0; c < 3; c++) {
      oldcorners[c] = &corners(oldtri, c);
      newcorners[c] = &corners(tri, c);
    }

    // move the position of the triangle
    mTris[tri] = mTris[oldtri];

    // 1) update c.node, c.opposite (c.next and c.prev should be fine as they are)
    for (int c = 0; c < 3; c++) {
      newcorners[c]->node = mTris[tri].c[c];
      newcorners[c]->opposite = oldcorners[c]->opposite;
    }

    //  2) c.opposite.opposite = c
    for (int c = 0; c < 3; c++) {
      if (newcorners[c]->opposite >= 0)
        mCorners[newcorners[c]->opposite].opposite = 3 * tri + c;
    }

    // update tri lookup
    for (int c = 0; c < 3; c++) {
      int node = mTris[tri].c[c];
      m1RingLookup[node].tris.erase(oldtri);
      m1RingLookup[node].tris.insert(tri);
    }
  }

  // transfer tri props
  for (size_t p = 0; p < mTriChannels.size(); p++)
    mTriChannels[p]->remove(tri);

  // pop the triangle and corners out of the vector
  mTris.pop_back();
  mCorners.resize(mTris.size() * 3);
}

void Mesh::removeNodes(const vector<int> &deletedNodes)
{
  // After we delete the nodes that are marked for removal,
  // the size of mNodes will be the current size - the size of the deleted array.
  // We are going to move the elements at the end of the array
  // (everything with an index >= newsize)
  // to the deleted spots.
  // We have to map all references to the last few nodes to their new locations.
  int newsize = (int)(mNodes.size() - deletedNodes.size());

  vector<int> new_index(deletedNodes.size());
  int di, ni;
  for (ni = 0; ni < (int)new_index.size(); ni++)
    new_index[ni] = 0;
  for (di = 0; di < (int)deletedNodes.size(); di++) {
    if (deletedNodes[di] >= newsize)
      new_index[deletedNodes[di] - newsize] = -1;  // tag this node as invalid
  }
  for (di = 0, ni = 0; ni < (int)new_index.size(); ni++, di++) {
    // we need to find a valid node to move
    // we marked invalid nodes in the earlier loop with a (-1),
    // so pick anything but those
    while (ni < (int)new_index.size() && new_index[ni] == -1)
      ni++;

    if (ni >= (int)new_index.size())
      break;

    // next we need to find a valid spot to move the node to.
    // we iterate through deleted[] until we find a valid spot
    while (di < (int)new_index.size() && deletedNodes[di] >= newsize)
      di++;

    // now we assign the valid node to the valid spot
    new_index[ni] = deletedNodes[di];
  }

  // Now we have a map of valid indices.
  // we move node[newsize+i] to location new_index[i].
  // We ignore the nodes with a -1 index, because they should not be moved.
  for (int i = 0; i < (int)new_index.size(); i++) {
    if (new_index[i] != -1)
      mNodes[new_index[i]] = mNodes[newsize + i];
  }
  mNodes.resize(newsize);

  // handle vertex properties
  for (size_t i = 0; i < mNodeChannels.size(); i++)
    mNodeChannels[i]->renumber(new_index, newsize);

  // finally, we reconnect everything that used to point to this vertex.
  for (size_t tri = 0, n = 0; tri < mTris.size(); tri++) {
    for (int c = 0; c < 3; c++, n++) {
      if (mCorners[n].node >= newsize) {
        int newindex = new_index[mCorners[n].node - newsize];
        mCorners[n].node = newindex;
        mTris[mCorners[n].tri].c[c] = newindex;
      }
    }
  }

  // renumber 1-ring
  for (int i = 0; i < (int)new_index.size(); i++) {
    if (new_index[i] != -1) {
      m1RingLookup[new_index[i]].nodes.swap(m1RingLookup[newsize + i].nodes);
      m1RingLookup[new_index[i]].tris.swap(m1RingLookup[newsize + i].tris);
    }
  }
  m1RingLookup.resize(newsize);
  vector<int> reStack(new_index.size());
  for (int i = 0; i < newsize; i++) {
    set<int> &cs = m1RingLookup[i].nodes;
    int reNum = 0;
    // find all nodes > newsize
    set<int>::reverse_iterator itend = cs.rend();
    for (set<int>::reverse_iterator it = cs.rbegin(); it != itend; ++it) {
      if (*it < newsize)
        break;
      reStack[reNum++] = *it;
    }
    // kill them and insert shifted values
    if (reNum > 0) {
      cs.erase(cs.find(reStack[reNum - 1]), cs.end());
      for (int j = 0; j < reNum; j++) {
        cs.insert(new_index[reStack[j] - newsize]);
#ifdef DEBUG
        if (new_index[reStack[j] - newsize] == -1)
          errMsg("invalid node present in 1-ring set");
#endif
      }
    }
  }
}

void Mesh::mergeNode(int node, int delnode)
{
  set<int> &ring = m1RingLookup[delnode].nodes;
  for (set<int>::iterator it = ring.begin(); it != ring.end(); ++it) {
    m1RingLookup[*it].nodes.erase(delnode);
    if (*it != node) {
      m1RingLookup[*it].nodes.insert(node);
      m1RingLookup[node].nodes.insert(*it);
    }
  }
  set<int> &ringt = m1RingLookup[delnode].tris;
  for (set<int>::iterator it = ringt.begin(); it != ringt.end(); ++it) {
    const int t = *it;
    for (int c = 0; c < 3; c++) {
      if (mCorners[3 * t + c].node == delnode) {
        mCorners[3 * t + c].node = node;
        mTris[t].c[c] = node;
      }
    }
    m1RingLookup[node].tris.insert(t);
  }
  for (size_t i = 0; i < mNodeChannels.size(); i++) {
    // weight is fixed to 1/2 for now
    mNodeChannels[i]->mergeWith(node, delnode, 0.5);
  }
}

void Mesh::removeTriFromLookup(int tri)
{
  for (int c = 0; c < 3; c++) {
    int node = mTris[tri].c[c];
    m1RingLookup[node].tris.erase(tri);
  }
}

void Mesh::addCorner(Corner a)
{
  mCorners.push_back(a);
}

int Mesh::addTri(Triangle a)
{
  mTris.push_back(a);
  for (int c = 0; c < 3; c++) {
    int node = a.c[c];
    int nextnode = a.c[(c + 1) % 3];
    if ((int)m1RingLookup.size() <= node)
      m1RingLookup.resize(node + 1);
    if ((int)m1RingLookup.size() <= nextnode)
      m1RingLookup.resize(nextnode + 1);
    m1RingLookup[node].nodes.insert(nextnode);
    m1RingLookup[nextnode].nodes.insert(node);
    m1RingLookup[node].tris.insert(mTris.size() - 1);
  }
  return mTris.size() - 1;
}

int Mesh::addNode(Node a)
{
  mNodes.push_back(a);
  if (m1RingLookup.size() < mNodes.size())
    m1RingLookup.resize(mNodes.size());

  // if mdata exists, add zero init for every node
  addAllMdata();

  return mNodes.size() - 1;
}

void Mesh::computeVertexNormals()
{
  for (size_t i = 0; i < mNodes.size(); i++) {
    mNodes[i].normal = 0.0;
  }
  for (size_t t = 0; t < mTris.size(); t++) {
    Vec3 p0 = getNode(t, 0), p1 = getNode(t, 1), p2 = getNode(t, 2);
    Vec3 n0 = p0 - p1, n1 = p1 - p2, n2 = p2 - p0;
    Real l0 = normSquare(n0), l1 = normSquare(n1), l2 = normSquare(n2);

    Vec3 nm = cross(n0, n1);

    mNodes[mTris[t].c[0]].normal += nm * (1.0 / (l0 * l2));
    mNodes[mTris[t].c[1]].normal += nm * (1.0 / (l0 * l1));
    mNodes[mTris[t].c[2]].normal += nm * (1.0 / (l1 * l2));
  }
  for (size_t i = 0; i < mNodes.size(); i++) {
    normalize(mNodes[i].normal);
  }
}

void Mesh::fastNodeLookupRebuild(int corner)
{
  int node = mCorners[corner].node;
  m1RingLookup[node].nodes.clear();
  m1RingLookup[node].tris.clear();
  int start = mCorners[corner].prev;
  int current = start;
  do {
    m1RingLookup[node].nodes.insert(mCorners[current].node);
    m1RingLookup[node].tris.insert(mCorners[current].tri);
    current = mCorners[mCorners[current].opposite].next;
    if (current < 0)
      errMsg("Can't use fastNodeLookupRebuild on incomplete surfaces");
  } while (current != start);
}

void Mesh::sanityCheck(bool strict, vector<int> *deletedNodes, map<int, bool> *taintedTris)
{
  const int nodes = numNodes(), tris = numTris(), corners = 3 * tris;
  for (size_t i = 0; i < mNodeChannels.size(); i++) {
    if (mNodeChannels[i]->size() != nodes)
      errMsg("Node channel size mismatch");
  }
  for (size_t i = 0; i < mTriChannels.size(); i++) {
    if (mTriChannels[i]->size() != tris)
      errMsg("Tri channel size mismatch");
  }
  if ((int)m1RingLookup.size() != nodes)
    errMsg("1Ring size wrong");
  for (size_t t = 0; t < mTris.size(); t++) {
    if (taintedTris && taintedTris->find(t) != taintedTris->end())
      continue;
    for (int c = 0; c < 3; c++) {
      int corner = t * 3 + c;
      int node = mTris[t].c[c];
      int next = mTris[t].c[(c + 1) % 3];
      int prev = mTris[t].c[(c + 2) % 3];
      int rnext = mCorners[corner].next;
      int rprev = mCorners[corner].prev;
      int ro = mCorners[corner].opposite;
      if (node < 0 || node >= nodes || next < 0 || next >= nodes || prev < 0 || prev >= nodes)
        errMsg("invalid node entry");
      if (mCorners[corner].node != node || mCorners[corner].tri != (int)t)
        errMsg("invalid basic corner entry");
      if (rnext < 0 || rnext >= corners || rprev < 0 || rprev >= corners || ro >= corners)
        errMsg("invalid corner links");
      if (mCorners[rnext].node != next || mCorners[rprev].node != prev)
        errMsg("invalid corner next/prev");
      if (strict && ro < 0)
        errMsg("opposite missing");
      if (mCorners[ro].opposite != corner)
        errMsg("invalid opposite ref");
      set<int> &rnodes = m1RingLookup[node].nodes;
      set<int> &rtris = m1RingLookup[node].tris;
      if (rnodes.find(next) == rnodes.end() || rnodes.find(prev) == rnodes.end()) {
        debMsg("Tri " << t << " " << node << " " << next << " " << prev, 1);
        for (set<int>::iterator it = rnodes.begin(); it != rnodes.end(); ++it)
          debMsg(*it, 1);
        errMsg("node missing in 1ring");
      }
      if (rtris.find(t) == rtris.end()) {
        debMsg("Tri " << t << " " << node, 1);
        errMsg("tri missing in 1ring");
      }
    }
  }
  for (int n = 0; n < nodes; n++) {
    bool docheck = true;
    if (deletedNodes)
      for (size_t e = 0; e < deletedNodes->size(); e++)
        if ((*deletedNodes)[e] == n)
          docheck = false;
    ;

    if (docheck) {
      set<int> &sn = m1RingLookup[n].nodes;
      set<int> &st = m1RingLookup[n].tris;
      set<int> sn2;

      for (set<int>::iterator it = st.begin(); it != st.end(); ++it) {
        bool found = false;
        for (int c = 0; c < 3; c++) {
          if (mTris[*it].c[c] == n)
            found = true;
          else
            sn2.insert(mTris[*it].c[c]);
        }
        if (!found) {
          cout << *it << " " << n << endl;
          for (int c = 0; c < 3; c++)
            cout << mTris[*it].c[c] << endl;
          errMsg("invalid triangle in 1ring");
        }
        if (taintedTris && taintedTris->find(*it) != taintedTris->end()) {
          cout << *it << endl;
          errMsg("tainted tri still is use");
        }
      }
      if (sn.size() != sn2.size())
        errMsg("invalid nodes in 1ring");
      for (set<int>::iterator it = sn.begin(), it2 = sn2.begin(); it != sn.end(); ++it, ++it2) {
        if (*it != *it2) {
          cout << "Node " << n << ": " << *it << " vs " << *it2 << endl;
          errMsg("node ring mismatch");
        }
      }
    }
  }
}

//*****************************************************************************
// rasterization

void meshSDF(Mesh &mesh, LevelsetGrid &levelset, Real sigma, Real cutoff = 0.);

//! helper vec3 array container
struct CVec3Ptr {
  Real *x, *y, *z;
  inline Vec3 get(int i) const
  {
    return Vec3(x[i], y[i], z[i]);
  };
  inline void set(int i, const Vec3 &v)
  {
    x[i] = v.x;
    y[i] = v.y;
    z[i] = v.z;
  };
};
//! helper vec3 array, for CUDA compatibility, remove at some point
struct CVec3Array {
  CVec3Array(int sz)
  {
    x.resize(sz);
    y.resize(sz);
    z.resize(sz);
  }
  CVec3Array(const std::vector<Vec3> &v)
  {
    x.resize(v.size());
    y.resize(v.size());
    z.resize(v.size());
    for (size_t i = 0; i < v.size(); i++) {
      x[i] = v[i].x;
      y[i] = v[i].y;
      z[i] = v[i].z;
    }
  }
  CVec3Ptr data()
  {
    CVec3Ptr a = {x.data(), y.data(), z.data()};
    return a;
  }
  inline const Vec3 operator[](int idx) const
  {
    return Vec3((Real)x[idx], (Real)y[idx], (Real)z[idx]);
  }
  inline void set(int idx, const Vec3 &v)
  {
    x[idx] = v.x;
    y[idx] = v.y;
    z[idx] = v.z;
  }
  inline int size()
  {
    return x.size();
  }
  std::vector<Real> x, y, z;
};

// void SDFKernel(const int* partStart, const int* partLen, CVec3Ptr pos, CVec3Ptr normal, Real*
// sdf, Vec3i gridRes, int intRadius, Real safeRadius2, Real cutoff2, Real isigma2);
//! helper for rasterization
static void SDFKernel(Grid<int> &partStart,
                      Grid<int> &partLen,
                      CVec3Ptr pos,
                      CVec3Ptr normal,
                      LevelsetGrid &sdf,
                      Vec3i gridRes,
                      int intRadius,
                      Real safeRadius2,
                      Real cutoff2,
                      Real isigma2)
{
  for (int cnt_x(0); cnt_x < gridRes[0]; ++cnt_x) {
    for (int cnt_y(0); cnt_y < gridRes[1]; ++cnt_y) {
      for (int cnt_z(0); cnt_z < gridRes[2]; ++cnt_z) {
        // cell index, center
        Vec3i cell = Vec3i(cnt_x, cnt_y, cnt_z);
        if (cell.x >= gridRes.x || cell.y >= gridRes.y || cell.z >= gridRes.z)
          return;
        Vec3 cpos = Vec3(cell.x + 0.5f, cell.y + 0.5f, cell.z + 0.5f);
        Real sum = 0.0f;
        Real dist = 0.0f;

        // query cells within block radius
        Vec3i minBlock = Vec3i(
            max(cell.x - intRadius, 0), max(cell.y - intRadius, 0), max(cell.z - intRadius, 0));
        Vec3i maxBlock = Vec3i(min(cell.x + intRadius, gridRes.x - 1),
                               min(cell.y + intRadius, gridRes.y - 1),
                               min(cell.z + intRadius, gridRes.z - 1));
        for (int i = minBlock.x; i <= maxBlock.x; i++)
          for (int j = minBlock.y; j <= maxBlock.y; j++)
            for (int k = minBlock.z; k <= maxBlock.z; k++) {
              // test if block is within radius
              Vec3 d = Vec3(cell.x - i, cell.y - j, cell.z - k);
              Real normSqr = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
              if (normSqr > safeRadius2)
                continue;

              // find source cell, and divide it into thread blocks
              int block = i + gridRes.x * (j + gridRes.y * k);
              int slen = partLen[block];
              if (slen == 0)
                continue;
              int start = partStart[block];

              // process sources
              for (int s = 0; s < slen; s++) {

                // actual sdf kernel
                Vec3 r = cpos - pos.get(start + s);
                Real normSqr = r[0] * r[0] + r[1] * r[1] + r[2] * r[2];
                Real r2 = normSqr;
                if (r2 < cutoff2) {
                  Real w = expf(-r2 * isigma2);
                  sum += w;
                  dist += dot(normal.get(start + s), r) * w;
                }
              }
            }
        // writeback
        if (sum > 0.0f) {
          // sdf[cell.x + gridRes.x * (cell.y + gridRes.y * cell.z)] = dist / sum;
          sdf(cell.x, cell.y, cell.z) = dist / sum;
        }
      }
    }
  }
}

static inline IndexInt _cIndex(const Vec3 &pos, const Vec3i &s)
{
  Vec3i p = toVec3i(pos);
  if (p.x < 0 || p.y < 0 || p.z < 0 || p.x >= s.x || p.y >= s.y || p.z >= s.z)
    return -1;
  return p.x + s.x * (p.y + s.y * p.z);
}

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

template<class T> struct ApplyMeshToGrid : public KernelBase {
  ApplyMeshToGrid(Grid<T> *grid, Grid<Real> &sdf, T value, FlagGrid *respectFlags)
      : KernelBase(grid, 0), grid(grid), sdf(sdf), value(value), respectFlags(respectFlags)
  {
    runMessage();
    run();
  }
  inline void op(
      int i, int j, int k, Grid<T> *grid, Grid<Real> &sdf, T value, FlagGrid *respectFlags) const
  {
    if (respectFlags && respectFlags->isObstacle(i, j, k))
      return;
    if (sdf(i, j, k) < 0) {
      (*grid)(i, j, k) = value;
    }
  }
  inline Grid<T> *getArg0()
  {
    return grid;
  }
  typedef Grid<T> type0;
  inline Grid<Real> &getArg1()
  {
    return sdf;
  }
  typedef Grid<Real> type1;
  inline T &getArg2()
  {
    return value;
  }
  typedef T type2;
  inline FlagGrid *getArg3()
  {
    return respectFlags;
  }
  typedef FlagGrid type3;
  void runMessage()
  {
    debMsg("Executing kernel ApplyMeshToGrid ", 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, sdf, 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, sdf, 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> &sdf;
  T value;
  FlagGrid *respectFlags;
};

void Mesh::applyMeshToGrid(GridBase *grid, FlagGrid *respectFlags, Real cutoff, Real meshSigma)
{
  FluidSolver dummy(grid->getSize());
  LevelsetGrid mesh_sdf(&dummy, false);
  meshSDF(*this, mesh_sdf, meshSigma, cutoff);  // meshSigma=2 fixed here

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

void Mesh::computeLevelset(LevelsetGrid &levelset, Real sigma, Real cutoff)
{
  meshSDF(*this, levelset, sigma, cutoff);
}

LevelsetGrid Mesh::getLevelset(Real sigma, Real cutoff)
{
  LevelsetGrid phi(getParent());
  meshSDF(*this, phi, sigma, cutoff);
  return phi;
}

void meshSDF(Mesh &mesh, LevelsetGrid &levelset, Real sigma, Real cutoff)
{
  if (cutoff < 0)
    cutoff = 2 * sigma;
  Real maxEdgeLength = 0.75;
  Real numSamplesPerCell = 0.75;

  Vec3i gridRes = levelset.getSize();
  Vec3 mult = toVec3(gridRes) / toVec3(mesh.getParent()->getGridSize());

  // prepare center values
  std::vector<Vec3> center;
  std::vector<Vec3> normals;
  short bigEdges(0);
  std::vector<Vec3> samplePoints;
  for (int i = 0; i < mesh.numTris(); i++) {
    center.push_back(Vec3(mesh.getFaceCenter(i) * mult));
    normals.push_back(mesh.getFaceNormal(i));
    // count big, stretched edges
    bigEdges = 0;
    for (short edge(0); edge < 3; ++edge) {
      if (norm(mesh.getEdge(i, edge)) > maxEdgeLength) {
        bigEdges += 1 << edge;
      }
    }
    if (bigEdges > 0) {
      samplePoints.clear();
      short iterA, pointA, iterB, pointB;
      int numSamples0 = norm(mesh.getEdge(i, 1)) * numSamplesPerCell;
      int numSamples1 = norm(mesh.getEdge(i, 2)) * numSamplesPerCell;
      int numSamples2 = norm(mesh.getEdge(i, 0)) * numSamplesPerCell;
      if (!(bigEdges & (1 << 0))) {
        // loop through 0,1
        iterA = numSamples1;
        pointA = 0;
        iterB = numSamples2;
        pointB = 1;
      }
      else if (!(bigEdges & (1 << 1))) {
        // loop through 1,2
        iterA = numSamples2;
        pointA = 1;
        iterB = numSamples0;
        pointB = 2;
      }
      else {
        // loop through 2,0
        iterA = numSamples0;
        pointA = 2;
        iterB = numSamples1;
        pointB = 0;
      }

      Real u(0.), v(0.), w(0.);  // barycentric uvw coords
      Vec3 samplePoint, normal;
      for (int sample0(0); sample0 < iterA; ++sample0) {
        u = Real(1. * sample0 / iterA);
        for (int sample1(0); sample1 < iterB; ++sample1) {
          v = Real(1. * sample1 / iterB);
          w = 1 - u - v;
          if (w < 0.)
            continue;
          samplePoint = mesh.getNode(i, pointA) * mult * u + mesh.getNode(i, pointB) * mult * v +
                        mesh.getNode(i, (3 - pointA - pointB)) * mult * w;
          samplePoints.push_back(samplePoint);
          normal = mesh.getFaceNormal(i);
          normals.push_back(normal);
        }
      }
      center.insert(center.end(), samplePoints.begin(), samplePoints.end());
    }
  }

  // prepare grid
  levelset.setConst(-cutoff);

  // 1. count sources per cell
  Grid<int> srcPerCell(levelset.getParent());
  for (size_t i = 0; i < center.size(); i++) {
    IndexInt idx = _cIndex(center[i], gridRes);
    if (idx >= 0)
      srcPerCell[idx]++;
  }

  // 2. create start index lookup
  Grid<int> srcCellStart(levelset.getParent());
  int cnt = 0;
  FOR_IJK(srcCellStart)
  {
    IndexInt idx = srcCellStart.index(i, j, k);
    srcCellStart[idx] = cnt;
    cnt += srcPerCell[idx];
  }

  // 3. reorder nodes
  CVec3Array reorderPos(center.size());
  CVec3Array reorderNormal(center.size());
  {
    Grid<int> curSrcCell(levelset.getParent());
    for (int i = 0; i < (int)center.size(); i++) {
      IndexInt idx = _cIndex(center[i], gridRes);
      if (idx < 0)
        continue;
      IndexInt idx2 = srcCellStart[idx] + curSrcCell[idx];
      reorderPos.set(idx2, center[i]);
      reorderNormal.set(idx2, normals[i]);
      curSrcCell[idx]++;
    }
  }

  // construct parameters
  Real safeRadius = cutoff + sqrt(3.0) * 0.5;
  Real safeRadius2 = safeRadius * safeRadius;
  Real cutoff2 = cutoff * cutoff;
  Real isigma2 = 1.0 / (sigma * sigma);
  int intRadius = (int)(cutoff + 0.5);

  SDFKernel(srcCellStart,
            srcPerCell,
            reorderPos.data(),
            reorderNormal.data(),
            levelset,
            gridRes,
            intRadius,
            safeRadius2,
            cutoff2,
            isigma2);

  // floodfill outside
  std::stack<Vec3i> outside;
  FOR_IJK(levelset)
  {
    if (levelset(i, j, k) >= cutoff - 1.0f)
      outside.push(Vec3i(i, j, k));
  }
  while (!outside.empty()) {
    Vec3i c = outside.top();
    outside.pop();
    levelset(c) = cutoff;
    if (c.x > 0 && levelset(c.x - 1, c.y, c.z) < 0)
      outside.push(Vec3i(c.x - 1, c.y, c.z));
    if (c.y > 0 && levelset(c.x, c.y - 1, c.z) < 0)
      outside.push(Vec3i(c.x, c.y - 1, c.z));
    if (c.z > 0 && levelset(c.x, c.y, c.z - 1) < 0)
      outside.push(Vec3i(c.x, c.y, c.z - 1));
    if (c.x < levelset.getSizeX() - 1 && levelset(c.x + 1, c.y, c.z) < 0)
      outside.push(Vec3i(c.x + 1, c.y, c.z));
    if (c.y < levelset.getSizeY() - 1 && levelset(c.x, c.y + 1, c.z) < 0)
      outside.push(Vec3i(c.x, c.y + 1, c.z));
    if (c.z < levelset.getSizeZ() - 1 && levelset(c.x, c.y, c.z + 1) < 0)
      outside.push(Vec3i(c.x, c.y, c.z + 1));
  };
}

// Blender data pointer accessors
std::string Mesh::getNodesDataPointer()
{
  std::ostringstream out;
  out << &mNodes;
  return out.str();
}
std::string Mesh::getTrisDataPointer()
{
  std::ostringstream out;
  out << &mTris;
  return out.str();
}

// mesh data

MeshDataBase::MeshDataBase(FluidSolver *parent) : PbClass(parent), mMesh(NULL)
{
}

MeshDataBase::~MeshDataBase()
{
  // notify parent of deletion
  if (mMesh)
    mMesh->deregister(this);
}

// actual data implementation

template<class T>
MeshDataImpl<T>::MeshDataImpl(FluidSolver *parent)
    : MeshDataBase(parent), mpGridSource(NULL), mGridSourceMAC(false)
{
}

template<class T>
MeshDataImpl<T>::MeshDataImpl(FluidSolver *parent, MeshDataImpl<T> *other)
    : MeshDataBase(parent), mpGridSource(NULL), mGridSourceMAC(false)
{
  this->mData = other->mData;
  setName(other->getName());
}

template<class T> MeshDataImpl<T>::~MeshDataImpl()
{
}

template<class T> IndexInt MeshDataImpl<T>::getSizeSlow() const
{
  return mData.size();
}
template<class T> void MeshDataImpl<T>::addEntry()
{
  // add zero'ed entry
  T tmp = T(0.);
  // for debugging, force init:
  // tmp = T(0.02 * mData.size()); // increasing
  // tmp = T(1.); // constant 1
  return mData.push_back(tmp);
}
template<class T> void MeshDataImpl<T>::resize(IndexInt s)
{
  mData.resize(s);
}
template<class T> void MeshDataImpl<T>::copyValueSlow(IndexInt from, IndexInt to)
{
  this->copyValue(from, to);
}
template<class T> MeshDataBase *MeshDataImpl<T>::clone()
{
  MeshDataImpl<T> *npd = new MeshDataImpl<T>(getParent(), this);
  return npd;
}

template<class T> void MeshDataImpl<T>::setSource(Grid<T> *grid, bool isMAC)
{
  mpGridSource = grid;
  mGridSourceMAC = isMAC;
  if (grid && isMAC)
    assertMsg(grid->getType() & GridBase::TypeMAC, "Given grid is not a valid MAC grid");
}

template<class T> void MeshDataImpl<T>::initNewValue(IndexInt idx, Vec3 pos)
{
  if (!mpGridSource)
    mData[idx] = 0;
  else {
    mData[idx] = mpGridSource->getInterpolated(pos);
  }
}

// special handling needed for velocities
template<> void MeshDataImpl<Vec3>::initNewValue(IndexInt idx, Vec3 pos)
{
  if (!mpGridSource)
    mData[idx] = 0;
  else {
    if (!mGridSourceMAC)
      mData[idx] = mpGridSource->getInterpolated(pos);
    else
      mData[idx] = ((MACGrid *)mpGridSource)->getInterpolated(pos);
  }
}

//! update additional mesh data
void Mesh::updateDataFields()
{
  for (size_t i = 0; i < mNodes.size(); ++i) {
    Vec3 pos = mNodes[i].pos;
    for (IndexInt md = 0; md < (IndexInt)mMdataReal.size(); ++md)
      mMdataReal[md]->initNewValue(i, pos);
    for (IndexInt md = 0; md < (IndexInt)mMdataVec3.size(); ++md)
      mMdataVec3[md]->initNewValue(i, pos);
    for (IndexInt md = 0; md < (IndexInt)mMdataInt.size(); ++md)
      mMdataInt[md]->initNewValue(i, pos);
  }
}

template<typename T> int MeshDataImpl<T>::load(string name)
{
  if (name.find_last_of('.') == string::npos)
    errMsg("file '" + name + "' does not have an extension");
  string ext = name.substr(name.find_last_of('.'));
  if (ext == ".uni")
    return readMdataUni<T>(name, this);
  else if (ext == ".raw")  // raw = uni for now
    return readMdataUni<T>(name, this);
  else
    errMsg("mesh data '" + name + "' filetype not supported for loading");
  return 0;
}

template<typename T> int MeshDataImpl<T>::save(string name)
{
  if (name.find_last_of('.') == string::npos)
    errMsg("file '" + name + "' does not have an extension");
  string ext = name.substr(name.find_last_of('.'));
  if (ext == ".uni")
    return writeMdataUni<T>(name, this);
  else if (ext == ".raw")  // raw = uni for now
    return writeMdataUni<T>(name, this);
  else
    errMsg("mesh data '" + name + "' filetype not supported for saving");
  return 0;
}

// specializations

template<> MeshDataBase::MdataType MeshDataImpl<Real>::getType() const
{
  return MeshDataBase::TypeReal;
}
template<> MeshDataBase::MdataType MeshDataImpl<int>::getType() const
{
  return MeshDataBase::TypeInt;
}
template<> MeshDataBase::MdataType MeshDataImpl<Vec3>::getType() const
{
  return MeshDataBase::TypeVec3;
}

template<class T> struct knSetMdataConst : public KernelBase {
  knSetMdataConst(MeshDataImpl<T> &mdata, T value)
      : KernelBase(mdata.size()), mdata(mdata), value(value)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &mdata, T value) const
  {
    mdata[idx] = value;
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return mdata;
  }
  typedef MeshDataImpl<T> type0;
  inline T &getArg1()
  {
    return value;
  }
  typedef T type1;
  void runMessage()
  {
    debMsg("Executing kernel knSetMdataConst ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, mdata, value);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &mdata;
  T value;
};

template<class T, class S> struct knMdataSet : public KernelBase {
  knMdataSet(MeshDataImpl<T> &me, const MeshDataImpl<S> &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const MeshDataImpl<S> &other) const
  {
    me[idx] += other[idx];
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<S> &getArg1()
  {
    return other;
  }
  typedef MeshDataImpl<S> type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataSet ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const MeshDataImpl<S> &other;
};
template<class T, class S> struct knMdataAdd : public KernelBase {
  knMdataAdd(MeshDataImpl<T> &me, const MeshDataImpl<S> &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const MeshDataImpl<S> &other) const
  {
    me[idx] += other[idx];
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<S> &getArg1()
  {
    return other;
  }
  typedef MeshDataImpl<S> type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataAdd ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const MeshDataImpl<S> &other;
};
template<class T, class S> struct knMdataSub : public KernelBase {
  knMdataSub(MeshDataImpl<T> &me, const MeshDataImpl<S> &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const MeshDataImpl<S> &other) const
  {
    me[idx] -= other[idx];
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<S> &getArg1()
  {
    return other;
  }
  typedef MeshDataImpl<S> type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataSub ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const MeshDataImpl<S> &other;
};
template<class T, class S> struct knMdataMult : public KernelBase {
  knMdataMult(MeshDataImpl<T> &me, const MeshDataImpl<S> &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const MeshDataImpl<S> &other) const
  {
    me[idx] *= other[idx];
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<S> &getArg1()
  {
    return other;
  }
  typedef MeshDataImpl<S> type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataMult ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const MeshDataImpl<S> &other;
};
template<class T, class S> struct knMdataDiv : public KernelBase {
  knMdataDiv(MeshDataImpl<T> &me, const MeshDataImpl<S> &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const MeshDataImpl<S> &other) const
  {
    me[idx] /= other[idx];
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<S> &getArg1()
  {
    return other;
  }
  typedef MeshDataImpl<S> type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataDiv ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const MeshDataImpl<S> &other;
};

template<class T, class S> struct knMdataSetScalar : public KernelBase {
  knMdataSetScalar(MeshDataImpl<T> &me, const S &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const S &other) const
  {
    me[idx] = other;
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const S &getArg1()
  {
    return other;
  }
  typedef S type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataSetScalar ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const S &other;
};
template<class T, class S> struct knMdataAddScalar : public KernelBase {
  knMdataAddScalar(MeshDataImpl<T> &me, const S &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const S &other) const
  {
    me[idx] += other;
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const S &getArg1()
  {
    return other;
  }
  typedef S type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataAddScalar ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const S &other;
};
template<class T, class S> struct knMdataMultScalar : public KernelBase {
  knMdataMultScalar(MeshDataImpl<T> &me, const S &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const S &other) const
  {
    me[idx] *= other;
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const S &getArg1()
  {
    return other;
  }
  typedef S type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataMultScalar ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const S &other;
};
template<class T, class S> struct knMdataScaledAdd : public KernelBase {
  knMdataScaledAdd(MeshDataImpl<T> &me, const MeshDataImpl<T> &other, const S &factor)
      : KernelBase(me.size()), me(me), other(other), factor(factor)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 MeshDataImpl<T> &me,
                 const MeshDataImpl<T> &other,
                 const S &factor) const
  {
    me[idx] += factor * other[idx];
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<T> &getArg1()
  {
    return other;
  }
  typedef MeshDataImpl<T> type1;
  inline const S &getArg2()
  {
    return factor;
  }
  typedef S type2;
  void runMessage()
  {
    debMsg("Executing kernel knMdataScaledAdd ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other, factor);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const MeshDataImpl<T> &other;
  const S &factor;
};

template<class T> struct knMdataSafeDiv : public KernelBase {
  knMdataSafeDiv(MeshDataImpl<T> &me, const MeshDataImpl<T> &other)
      : KernelBase(me.size()), me(me), other(other)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const MeshDataImpl<T> &other) const
  {
    me[idx] = safeDivide(me[idx], other[idx]);
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<T> &getArg1()
  {
    return other;
  }
  typedef MeshDataImpl<T> type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataSafeDiv ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const MeshDataImpl<T> &other;
};
template<class T> struct knMdataSetConst : public KernelBase {
  knMdataSetConst(MeshDataImpl<T> &mdata, T value)
      : KernelBase(mdata.size()), mdata(mdata), value(value)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &mdata, T value) const
  {
    mdata[idx] = value;
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return mdata;
  }
  typedef MeshDataImpl<T> type0;
  inline T &getArg1()
  {
    return value;
  }
  typedef T type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataSetConst ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, mdata, value);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &mdata;
  T value;
};

template<class T> struct knMdataClamp : public KernelBase {
  knMdataClamp(MeshDataImpl<T> &me, T min, T max)
      : KernelBase(me.size()), me(me), min(min), max(max)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, T min, T max) const
  {
    me[idx] = clamp(me[idx], min, max);
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline T &getArg1()
  {
    return min;
  }
  typedef T type1;
  inline T &getArg2()
  {
    return max;
  }
  typedef T type2;
  void runMessage()
  {
    debMsg("Executing kernel knMdataClamp ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, min, max);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  T min;
  T max;
};
template<class T> struct knMdataClampMin : public KernelBase {
  knMdataClampMin(MeshDataImpl<T> &me, const T vmin) : KernelBase(me.size()), me(me), vmin(vmin)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const T vmin) const
  {
    me[idx] = std::max(vmin, me[idx]);
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const T &getArg1()
  {
    return vmin;
  }
  typedef T type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataClampMin ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, vmin);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const T vmin;
};
template<class T> struct knMdataClampMax : public KernelBase {
  knMdataClampMax(MeshDataImpl<T> &me, const T vmax) : KernelBase(me.size()), me(me), vmax(vmax)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<T> &me, const T vmax) const
  {
    me[idx] = std::min(vmax, me[idx]);
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const T &getArg1()
  {
    return vmax;
  }
  typedef T type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataClampMax ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, vmax);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const T vmax;
};
struct knMdataClampMinVec3 : public KernelBase {
  knMdataClampMinVec3(MeshDataImpl<Vec3> &me, const Real vmin)
      : KernelBase(me.size()), me(me), vmin(vmin)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<Vec3> &me, const Real vmin) const
  {
    me[idx].x = std::max(vmin, me[idx].x);
    me[idx].y = std::max(vmin, me[idx].y);
    me[idx].z = std::max(vmin, me[idx].z);
  }
  inline MeshDataImpl<Vec3> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<Vec3> type0;
  inline const Real &getArg1()
  {
    return vmin;
  }
  typedef Real type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataClampMinVec3 ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, vmin);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<Vec3> &me;
  const Real vmin;
};
struct knMdataClampMaxVec3 : public KernelBase {
  knMdataClampMaxVec3(MeshDataImpl<Vec3> &me, const Real vmax)
      : KernelBase(me.size()), me(me), vmax(vmax)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, MeshDataImpl<Vec3> &me, const Real vmax) const
  {
    me[idx].x = std::min(vmax, me[idx].x);
    me[idx].y = std::min(vmax, me[idx].y);
    me[idx].z = std::min(vmax, me[idx].z);
  }
  inline MeshDataImpl<Vec3> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<Vec3> type0;
  inline const Real &getArg1()
  {
    return vmax;
  }
  typedef Real type1;
  void runMessage()
  {
    debMsg("Executing kernel knMdataClampMaxVec3 ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, vmax);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<Vec3> &me;
  const Real vmax;
};

// python operators

template<typename T> MeshDataImpl<T> &MeshDataImpl<T>::copyFrom(const MeshDataImpl<T> &a)
{
  assertMsg(a.mData.size() == mData.size(),
            "different mdata size " << a.mData.size() << " vs " << this->mData.size());
  memcpy(&mData[0], &a.mData[0], sizeof(T) * mData.size());
  return *this;
}

template<typename T> void MeshDataImpl<T>::setConst(T s)
{
  knMdataSetScalar<T, T> op(*this, s);
}

template<typename T> void MeshDataImpl<T>::setConstRange(T s, const int begin, const int end)
{
  for (int i = begin; i < end; ++i)
    (*this)[i] = s;
}

// special set by flag
template<class T, class S> struct knMdataSetScalarIntFlag : public KernelBase {
  knMdataSetScalarIntFlag(MeshDataImpl<T> &me,
                          const S &other,
                          const MeshDataImpl<int> &t,
                          const int itype)
      : KernelBase(me.size()), me(me), other(other), t(t), itype(itype)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 MeshDataImpl<T> &me,
                 const S &other,
                 const MeshDataImpl<int> &t,
                 const int itype) const
  {
    if (t[idx] & itype)
      me[idx] = other;
  }
  inline MeshDataImpl<T> &getArg0()
  {
    return me;
  }
  typedef MeshDataImpl<T> type0;
  inline const S &getArg1()
  {
    return other;
  }
  typedef S type1;
  inline const MeshDataImpl<int> &getArg2()
  {
    return t;
  }
  typedef MeshDataImpl<int> type2;
  inline const int &getArg3()
  {
    return itype;
  }
  typedef int type3;
  void runMessage()
  {
    debMsg("Executing kernel knMdataSetScalarIntFlag ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r) const
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, me, other, t, itype);
  }
  void run()
  {
    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  MeshDataImpl<T> &me;
  const S &other;
  const MeshDataImpl<int> &t;
  const int itype;
};
template<typename T>
void MeshDataImpl<T>::setConstIntFlag(T s, const MeshDataImpl<int> &t, const int itype)
{
  knMdataSetScalarIntFlag<T, T> op(*this, s, t, itype);
}

template<typename T> void MeshDataImpl<T>::add(const MeshDataImpl<T> &a)
{
  knMdataAdd<T, T> op(*this, a);
}
template<typename T> void MeshDataImpl<T>::sub(const MeshDataImpl<T> &a)
{
  knMdataSub<T, T> op(*this, a);
}

template<typename T> void MeshDataImpl<T>::addConst(T s)
{
  knMdataAddScalar<T, T> op(*this, s);
}

template<typename T> void MeshDataImpl<T>::addScaled(const MeshDataImpl<T> &a, const T &factor)
{
  knMdataScaledAdd<T, T> op(*this, a, factor);
}

template<typename T> void MeshDataImpl<T>::mult(const MeshDataImpl<T> &a)
{
  knMdataMult<T, T> op(*this, a);
}

template<typename T> void MeshDataImpl<T>::safeDiv(const MeshDataImpl<T> &a)
{
  knMdataSafeDiv<T> op(*this, a);
}

template<typename T> void MeshDataImpl<T>::multConst(T s)
{
  knMdataMultScalar<T, T> op(*this, s);
}

template<typename T> void MeshDataImpl<T>::clamp(Real vmin, Real vmax)
{
  knMdataClamp<T> op(*this, vmin, vmax);
}

template<typename T> void MeshDataImpl<T>::clampMin(Real vmin)
{
  knMdataClampMin<T> op(*this, vmin);
}
template<typename T> void MeshDataImpl<T>::clampMax(Real vmax)
{
  knMdataClampMax<T> op(*this, vmax);
}

template<> void MeshDataImpl<Vec3>::clampMin(Real vmin)
{
  knMdataClampMinVec3 op(*this, vmin);
}
template<> void MeshDataImpl<Vec3>::clampMax(Real vmax)
{
  knMdataClampMaxVec3 op(*this, vmax);
}

template<typename T> struct KnPtsSum : public KernelBase {
  KnPtsSum(const MeshDataImpl<T> &val, const MeshDataImpl<int> *t, const int itype)
      : KernelBase(val.size()), val(val), t(t), itype(itype), result(T(0.))
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx,
                 const MeshDataImpl<T> &val,
                 const MeshDataImpl<int> *t,
                 const int itype,
                 T &result)
  {
    if (t && !((*t)[idx] & itype))
      return;
    result += val[idx];
  }
  inline operator T()
  {
    return result;
  }
  inline T &getRet()
  {
    return result;
  }
  inline const MeshDataImpl<T> &getArg0()
  {
    return val;
  }
  typedef MeshDataImpl<T> type0;
  inline const MeshDataImpl<int> *getArg1()
  {
    return t;
  }
  typedef MeshDataImpl<int> type1;
  inline const int &getArg2()
  {
    return itype;
  }
  typedef int type2;
  void runMessage()
  {
    debMsg("Executing kernel KnPtsSum ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, val, t, itype, result);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  KnPtsSum(KnPtsSum &o, tbb::split)
      : KernelBase(o), val(o.val), t(o.t), itype(o.itype), result(T(0.))
  {
  }
  void join(const KnPtsSum &o)
  {
    result += o.result;
  }
  const MeshDataImpl<T> &val;
  const MeshDataImpl<int> *t;
  const int itype;
  T result;
};
template<typename T> struct KnPtsSumSquare : public KernelBase {
  KnPtsSumSquare(const MeshDataImpl<T> &val) : KernelBase(val.size()), val(val), result(0.)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, const MeshDataImpl<T> &val, Real &result)
  {
    result += normSquare(val[idx]);
  }
  inline operator Real()
  {
    return result;
  }
  inline Real &getRet()
  {
    return result;
  }
  inline const MeshDataImpl<T> &getArg0()
  {
    return val;
  }
  typedef MeshDataImpl<T> type0;
  void runMessage()
  {
    debMsg("Executing kernel KnPtsSumSquare ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, val, result);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  KnPtsSumSquare(KnPtsSumSquare &o, tbb::split) : KernelBase(o), val(o.val), result(0.)
  {
  }
  void join(const KnPtsSumSquare &o)
  {
    result += o.result;
  }
  const MeshDataImpl<T> &val;
  Real result;
};
template<typename T> struct KnPtsSumMagnitude : public KernelBase {
  KnPtsSumMagnitude(const MeshDataImpl<T> &val) : KernelBase(val.size()), val(val), result(0.)
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, const MeshDataImpl<T> &val, Real &result)
  {
    result += norm(val[idx]);
  }
  inline operator Real()
  {
    return result;
  }
  inline Real &getRet()
  {
    return result;
  }
  inline const MeshDataImpl<T> &getArg0()
  {
    return val;
  }
  typedef MeshDataImpl<T> type0;
  void runMessage()
  {
    debMsg("Executing kernel KnPtsSumMagnitude ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, val, result);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  KnPtsSumMagnitude(KnPtsSumMagnitude &o, tbb::split) : KernelBase(o), val(o.val), result(0.)
  {
  }
  void join(const KnPtsSumMagnitude &o)
  {
    result += o.result;
  }
  const MeshDataImpl<T> &val;
  Real result;
};

template<typename T> T MeshDataImpl<T>::sum(const MeshDataImpl<int> *t, const int itype) const
{
  return KnPtsSum<T>(*this, t, itype);
}
template<typename T> Real MeshDataImpl<T>::sumSquare() const
{
  return KnPtsSumSquare<T>(*this);
}
template<typename T> Real MeshDataImpl<T>::sumMagnitude() const
{
  return KnPtsSumMagnitude<T>(*this);
}

template<typename T>

struct CompMdata_Min : public KernelBase {
  CompMdata_Min(const MeshDataImpl<T> &val)
      : KernelBase(val.size()), val(val), minVal(std::numeric_limits<Real>::max())
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, const MeshDataImpl<T> &val, Real &minVal)
  {
    if (val[idx] < minVal)
      minVal = val[idx];
  }
  inline operator Real()
  {
    return minVal;
  }
  inline Real &getRet()
  {
    return minVal;
  }
  inline const MeshDataImpl<T> &getArg0()
  {
    return val;
  }
  typedef MeshDataImpl<T> type0;
  void runMessage()
  {
    debMsg("Executing kernel CompMdata_Min ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, val, minVal);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  CompMdata_Min(CompMdata_Min &o, tbb::split)
      : KernelBase(o), val(o.val), minVal(std::numeric_limits<Real>::max())
  {
  }
  void join(const CompMdata_Min &o)
  {
    minVal = min(minVal, o.minVal);
  }
  const MeshDataImpl<T> &val;
  Real minVal;
};

template<typename T>

struct CompMdata_Max : public KernelBase {
  CompMdata_Max(const MeshDataImpl<T> &val)
      : KernelBase(val.size()), val(val), maxVal(-std::numeric_limits<Real>::max())
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, const MeshDataImpl<T> &val, Real &maxVal)
  {
    if (val[idx] > maxVal)
      maxVal = val[idx];
  }
  inline operator Real()
  {
    return maxVal;
  }
  inline Real &getRet()
  {
    return maxVal;
  }
  inline const MeshDataImpl<T> &getArg0()
  {
    return val;
  }
  typedef MeshDataImpl<T> type0;
  void runMessage()
  {
    debMsg("Executing kernel CompMdata_Max ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, val, maxVal);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  CompMdata_Max(CompMdata_Max &o, tbb::split)
      : KernelBase(o), val(o.val), maxVal(-std::numeric_limits<Real>::max())
  {
  }
  void join(const CompMdata_Max &o)
  {
    maxVal = max(maxVal, o.maxVal);
  }
  const MeshDataImpl<T> &val;
  Real maxVal;
};

template<typename T> Real MeshDataImpl<T>::getMin()
{
  return CompMdata_Min<T>(*this);
}

template<typename T> Real MeshDataImpl<T>::getMaxAbs()
{
  Real amin = CompMdata_Min<T>(*this);
  Real amax = CompMdata_Max<T>(*this);
  return max(fabs(amin), fabs(amax));
}

template<typename T> Real MeshDataImpl<T>::getMax()
{
  return CompMdata_Max<T>(*this);
}

template<typename T>
void MeshDataImpl<T>::printMdata(IndexInt start, IndexInt stop, bool printIndex)
{
  std::ostringstream sstr;
  IndexInt s = (start > 0 ? start : 0);
  IndexInt e = (stop > 0 ? stop : (IndexInt)mData.size());
  s = Manta::clamp(s, (IndexInt)0, (IndexInt)mData.size());
  e = Manta::clamp(e, (IndexInt)0, (IndexInt)mData.size());

  for (IndexInt i = s; i < e; ++i) {
    if (printIndex)
      sstr << i << ": ";
    sstr << mData[i] << " "
         << "\n";
  }
  debMsg(sstr.str(), 1);
}
template<class T> std::string MeshDataImpl<T>::getDataPointer()
{
  std::ostringstream out;
  out << &mData;
  return out.str();
}

// specials for vec3
// work on length values, ie, always positive (in contrast to scalar versions above)

struct CompMdata_MinVec3 : public KernelBase {
  CompMdata_MinVec3(const MeshDataImpl<Vec3> &val)
      : KernelBase(val.size()), val(val), minVal(-std::numeric_limits<Real>::max())
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, const MeshDataImpl<Vec3> &val, Real &minVal)
  {
    const Real s = normSquare(val[idx]);
    if (s < minVal)
      minVal = s;
  }
  inline operator Real()
  {
    return minVal;
  }
  inline Real &getRet()
  {
    return minVal;
  }
  inline const MeshDataImpl<Vec3> &getArg0()
  {
    return val;
  }
  typedef MeshDataImpl<Vec3> type0;
  void runMessage()
  {
    debMsg("Executing kernel CompMdata_MinVec3 ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, val, minVal);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  CompMdata_MinVec3(CompMdata_MinVec3 &o, tbb::split)
      : KernelBase(o), val(o.val), minVal(-std::numeric_limits<Real>::max())
  {
  }
  void join(const CompMdata_MinVec3 &o)
  {
    minVal = min(minVal, o.minVal);
  }
  const MeshDataImpl<Vec3> &val;
  Real minVal;
};

struct CompMdata_MaxVec3 : public KernelBase {
  CompMdata_MaxVec3(const MeshDataImpl<Vec3> &val)
      : KernelBase(val.size()), val(val), maxVal(-std::numeric_limits<Real>::min())
  {
    runMessage();
    run();
  }
  inline void op(IndexInt idx, const MeshDataImpl<Vec3> &val, Real &maxVal)
  {
    const Real s = normSquare(val[idx]);
    if (s > maxVal)
      maxVal = s;
  }
  inline operator Real()
  {
    return maxVal;
  }
  inline Real &getRet()
  {
    return maxVal;
  }
  inline const MeshDataImpl<Vec3> &getArg0()
  {
    return val;
  }
  typedef MeshDataImpl<Vec3> type0;
  void runMessage()
  {
    debMsg("Executing kernel CompMdata_MaxVec3 ", 3);
    debMsg("Kernel range"
               << " size " << size << " ",
           4);
  };
  void operator()(const tbb::blocked_range<IndexInt> &__r)
  {
    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
      op(idx, val, maxVal);
  }
  void run()
  {
    tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this);
  }
  CompMdata_MaxVec3(CompMdata_MaxVec3 &o, tbb::split)
      : KernelBase(o), val(o.val), maxVal(-std::numeric_limits<Real>::min())
  {
  }
  void join(const CompMdata_MaxVec3 &o)
  {
    maxVal = max(maxVal, o.maxVal);
  }
  const MeshDataImpl<Vec3> &val;
  Real maxVal;
};

template<> Real MeshDataImpl<Vec3>::getMin()
{
  return sqrt(CompMdata_MinVec3(*this));
}

template<> Real MeshDataImpl<Vec3>::getMaxAbs()
{
  return sqrt(CompMdata_MaxVec3(*this));  // no minimum necessary here
}

template<> Real MeshDataImpl<Vec3>::getMax()
{
  return sqrt(CompMdata_MaxVec3(*this));
}

// explicit instantiation
template class MeshDataImpl<int>;
template class MeshDataImpl<Real>;
template class MeshDataImpl<Vec3>;

}  // namespace Manta