Page MenuHome
Differential D3850 Diff 17619 intern/mantaflow/intern/manta_develop/preprocessed/tbb/mesh.cpp

Changeset View

Standalone View

View Options

intern/mantaflow/intern/manta_develop/preprocessed/tbb/mesh.cpp

  • This file was added.
// 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)
{
Vec3 velo = mNodes[i].pos - oldMesh.mNodes[i].pos;
if(bIs2D && (mNodes[i].pos.z > -0.5f && mNodes[i].pos.z < 0.5f))
{
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) {
FluidSolver dummy(grid->getSize());
LevelsetGrid mesh_sdf(&dummy, false);
meshSDF(*this, mesh_sdf, 2., 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