Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/extern/mantaflow/preprocessed/mesh.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'extern/mantaflow/preprocessed/mesh.cpp')
-rw-r--r--extern/mantaflow/preprocessed/mesh.cpp2733
1 files changed, 2733 insertions, 0 deletions
diff --git a/extern/mantaflow/preprocessed/mesh.cpp b/extern/mantaflow/preprocessed/mesh.cpp
new file mode 100644
index 00000000000..d93c2ac04c0
--- /dev/null
+++ b/extern/mantaflow/preprocessed/mesh.cpp
@@ -0,0 +1,2733 @@
+
+
+// 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;
+}
+
+void 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
+ readBobjFile(name, this, append);
+ else if (ext == ".obj")
+ readObjFile(name, this, append);
+ else
+ errMsg("file '" + name + "' filetype not supported");
+
+ // dont always rebuild...
+ // rebuildCorners();
+ // rebuildLookup();
+}
+
+void 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")
+ writeObjFile(name, this);
+ else if (ext == ".gz")
+ writeBobjFile(name, this);
+ else
+ errMsg("file '" + name + "' filetype not supported");
+}
+
+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 (isMAC)
+ assertMsg(dynamic_cast<MACGrid *>(grid) != NULL, "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, mNodes[i].pos);
+ for (IndexInt md = 0; md < (IndexInt)mMdataVec3.size(); ++md)
+ mMdataVec3[md]->initNewValue(i, mNodes[i].pos);
+ for (IndexInt md = 0; md < (IndexInt)mMdataInt.size(); ++md)
+ mMdataInt[md]->initNewValue(i, mNodes[i].pos);
+ }
+}
+
+template<typename T> void 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")
+ readMdataUni<T>(name, this);
+ else if (ext == ".raw") // raw = uni for now
+ readMdataUni<T>(name, this);
+ else
+ errMsg("mesh data '" + name + "' filetype not supported for loading");
+}
+
+template<typename T> void 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")
+ writeMdataUni<T>(name, this);
+ else if (ext == ".raw") // raw = uni for now
+ writeMdataUni<T>(name, this);
+ else
+ errMsg("mesh data '" + name + "' filetype not supported for saving");
+}
+
+// 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
generated by cgit v1.2.3 (git 2.25.1) at 2024-08-26 07:44:24 +0300