diff options
Diffstat (limited to 'extern/mantaflow/preprocessed/mesh.cpp')
-rw-r--r-- | extern/mantaflow/preprocessed/mesh.cpp | 2733 |
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 |