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

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/extern/mantaflow/preprocessed/particle.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'extern/mantaflow/preprocessed/particle.cpp')
-rw-r--r--extern/mantaflow/preprocessed/particle.cpp1620
1 files changed, 1620 insertions, 0 deletions
diff --git a/extern/mantaflow/preprocessed/particle.cpp b/extern/mantaflow/preprocessed/particle.cpp
new file mode 100644
index 00000000000..478f1417109
--- /dev/null
+++ b/extern/mantaflow/preprocessed/particle.cpp
@@ -0,0 +1,1620 @@
+
+
+// DO NOT EDIT !
+// This file is generated using the MantaFlow preprocessor (prep generate).
+
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2013 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
+ *
+ * Particle data functionality
+ *
+ ******************************************************************************/
+
+#include <fstream>
+#include <cstring>
+#if NO_ZLIB != 1
+# include <zlib.h>
+#endif
+#include "particle.h"
+#include "levelset.h"
+#include "mantaio.h"
+
+using namespace std;
+namespace Manta {
+
+ParticleBase::ParticleBase(FluidSolver *parent)
+ : PbClass(parent), mAllowCompress(true), mFreePdata(false)
+{
+}
+
+ParticleBase::~ParticleBase()
+{
+ // make sure data fields now parent system is deleted
+ for (IndexInt i = 0; i < (IndexInt)mPartData.size(); ++i)
+ mPartData[i]->setParticleSys(NULL);
+
+ if (mFreePdata) {
+ for (IndexInt i = 0; i < (IndexInt)mPartData.size(); ++i)
+ delete mPartData[i];
+ }
+}
+
+std::string ParticleBase::infoString() const
+{
+ return "ParticleSystem " + mName + " <no info>";
+}
+
+void ParticleBase::cloneParticleData(ParticleBase *nm)
+{
+ // clone additional data , and make sure the copied particle system deletes it
+ nm->mFreePdata = true;
+ for (IndexInt i = 0; i < (IndexInt)mPartData.size(); ++i) {
+ ParticleDataBase *pdata = mPartData[i]->clone();
+ nm->registerPdata(pdata);
+ }
+}
+
+void ParticleBase::deregister(ParticleDataBase *pdata)
+{
+ bool done = false;
+ // remove pointer from particle data list
+ for (IndexInt i = 0; i < (IndexInt)mPartData.size(); ++i) {
+ if (mPartData[i] == pdata) {
+ if (i < (IndexInt)mPartData.size() - 1)
+ mPartData[i] = mPartData[mPartData.size() - 1];
+ mPartData.pop_back();
+ done = true;
+ }
+ }
+ if (!done)
+ errMsg("Invalid pointer given, not registered!");
+}
+
+// create and attach a new pdata field to this particle system
+PbClass *ParticleBase::create(PbType t, PbTypeVec T, const string &name)
+{
+#if NOPYTHON != 1
+ _args.add("nocheck", true);
+ if (t.str() == "")
+ errMsg("Specify particle data type to create");
+ // debMsg( "Pdata creating '"<< t.str <<" with size "<< this->getSizeSlow(), 5 );
+
+ PbClass *pyObj = PbClass::createPyObject(t.str() + T.str(), name, _args, this->getParent());
+
+ ParticleDataBase *pdata = dynamic_cast<ParticleDataBase *>(pyObj);
+ if (!pdata) {
+ errMsg(
+ "Unable to get particle data pointer from newly created object. Only create ParticleData "
+ "type with a ParticleSys.creat() call, eg, PdataReal, PdataVec3 etc.");
+ delete pyObj;
+ return NULL;
+ }
+ else {
+ this->registerPdata(pdata);
+ }
+
+ // directly init size of new pdata field:
+ pdata->resize(this->getSizeSlow());
+#else
+ PbClass *pyObj = NULL;
+#endif
+ return pyObj;
+}
+
+void ParticleBase::registerPdata(ParticleDataBase *pdata)
+{
+ pdata->setParticleSys(this);
+ mPartData.push_back(pdata);
+
+ if (pdata->getType() == ParticleDataBase::TypeReal) {
+ ParticleDataImpl<Real> *pd = dynamic_cast<ParticleDataImpl<Real> *>(pdata);
+ if (!pd)
+ errMsg("Invalid pdata object posing as real!");
+ this->registerPdataReal(pd);
+ }
+ else if (pdata->getType() == ParticleDataBase::TypeInt) {
+ ParticleDataImpl<int> *pd = dynamic_cast<ParticleDataImpl<int> *>(pdata);
+ if (!pd)
+ errMsg("Invalid pdata object posing as int!");
+ this->registerPdataInt(pd);
+ }
+ else if (pdata->getType() == ParticleDataBase::TypeVec3) {
+ ParticleDataImpl<Vec3> *pd = dynamic_cast<ParticleDataImpl<Vec3> *>(pdata);
+ if (!pd)
+ errMsg("Invalid pdata object posing as vec3!");
+ this->registerPdataVec3(pd);
+ }
+}
+void ParticleBase::registerPdataReal(ParticleDataImpl<Real> *pd)
+{
+ mPdataReal.push_back(pd);
+}
+void ParticleBase::registerPdataVec3(ParticleDataImpl<Vec3> *pd)
+{
+ mPdataVec3.push_back(pd);
+}
+void ParticleBase::registerPdataInt(ParticleDataImpl<int> *pd)
+{
+ mPdataInt.push_back(pd);
+}
+
+void ParticleBase::addAllPdata()
+{
+ for (IndexInt i = 0; i < (IndexInt)mPartData.size(); ++i) {
+ mPartData[i]->addEntry();
+ }
+}
+
+BasicParticleSystem::BasicParticleSystem(FluidSolver *parent)
+ : ParticleSystem<BasicParticleData>(parent)
+{
+ this->mAllowCompress = false;
+}
+
+// file io
+
+void BasicParticleSystem::writeParticlesText(const string name) const
+{
+ ofstream ofs(name.c_str());
+ if (!ofs.good())
+ errMsg("can't open file!");
+ ofs << this->size() << ", pdata: " << mPartData.size() << " (" << mPdataInt.size() << ","
+ << mPdataReal.size() << "," << mPdataVec3.size() << ") \n";
+ for (IndexInt i = 0; i < this->size(); ++i) {
+ ofs << i << ": " << this->getPos(i) << " , " << this->getStatus(i) << ". ";
+ for (IndexInt pd = 0; pd < (IndexInt)mPdataInt.size(); ++pd)
+ ofs << mPdataInt[pd]->get(i) << " ";
+ for (IndexInt pd = 0; pd < (IndexInt)mPdataReal.size(); ++pd)
+ ofs << mPdataReal[pd]->get(i) << " ";
+ for (IndexInt pd = 0; pd < (IndexInt)mPdataVec3.size(); ++pd)
+ ofs << mPdataVec3[pd]->get(i) << " ";
+ ofs << "\n";
+ }
+ ofs.close();
+}
+
+void BasicParticleSystem::writeParticlesRawPositionsGz(const string name) const
+{
+#if NO_ZLIB != 1
+ gzFile gzf = gzopen(name.c_str(), "wb1");
+ if (!gzf)
+ errMsg("can't open file " << name);
+ for (IndexInt i = 0; i < this->size(); ++i) {
+ Vector3D<float> p = toVec3f(this->getPos(i));
+ gzwrite(gzf, &p, sizeof(float) * 3);
+ }
+ gzclose(gzf);
+#else
+ cout << "file format not supported without zlib" << endl;
+#endif
+}
+
+void BasicParticleSystem::writeParticlesRawVelocityGz(const string name) const
+{
+#if NO_ZLIB != 1
+ gzFile gzf = gzopen(name.c_str(), "wb1");
+ if (!gzf)
+ errMsg("can't open file " << name);
+ if (mPdataVec3.size() < 1)
+ errMsg("no vec3 particle data channel found!");
+ // note , assuming particle data vec3 0 is velocity! make optional...
+ for (IndexInt i = 0; i < this->size(); ++i) {
+ Vector3D<float> p = toVec3f(mPdataVec3[0]->get(i));
+ gzwrite(gzf, &p, sizeof(float) * 3);
+ }
+ gzclose(gzf);
+#else
+ cout << "file format not supported without zlib" << endl;
+#endif
+}
+
+void BasicParticleSystem::load(const 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")
+ readParticlesUni(name, this);
+ else if (ext == ".raw") // raw = uni for now
+ readParticlesUni(name, this);
+ else
+ errMsg("particle '" + name + "' filetype not supported for loading");
+}
+
+void BasicParticleSystem::save(const string name) const
+{
+ 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 == ".txt")
+ this->writeParticlesText(name);
+ else if (ext == ".uni")
+ writeParticlesUni(name, this);
+ else if (ext == ".raw") // raw = uni for now
+ writeParticlesUni(name, this);
+ // raw data formats, very basic for simple data transfer to other programs
+ else if (ext == ".posgz")
+ this->writeParticlesRawPositionsGz(name);
+ else if (ext == ".velgz")
+ this->writeParticlesRawVelocityGz(name);
+ else
+ errMsg("particle '" + name + "' filetype not supported for saving");
+}
+
+void BasicParticleSystem::printParts(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].pos << " " << mData[i].flag << "\n";
+ }
+ debMsg(sstr.str(), 1);
+}
+
+std::string BasicParticleSystem::getDataPointer()
+{
+ std::ostringstream out;
+ out << &mData;
+ return out.str();
+}
+
+void BasicParticleSystem::readParticles(BasicParticleSystem *from)
+{
+ // re-allocate all data
+ this->resizeAll(from->size());
+ assertMsg(from->size() == this->size(), "particle size doesn't match");
+
+ for (int i = 0; i < this->size(); ++i) {
+ (*this)[i].pos = (*from)[i].pos;
+ (*this)[i].flag = (*from)[i].flag;
+ }
+ this->transformPositions(from->getParent()->getGridSize(), this->getParent()->getGridSize());
+}
+
+// particle data
+
+ParticleDataBase::ParticleDataBase(FluidSolver *parent) : PbClass(parent), mpParticleSys(NULL)
+{
+}
+
+ParticleDataBase::~ParticleDataBase()
+{
+ // notify parent of deletion
+ if (mpParticleSys)
+ mpParticleSys->deregister(this);
+}
+
+// actual data implementation
+
+template<class T>
+ParticleDataImpl<T>::ParticleDataImpl(FluidSolver *parent)
+ : ParticleDataBase(parent), mpGridSource(NULL), mGridSourceMAC(false)
+{
+}
+
+template<class T>
+ParticleDataImpl<T>::ParticleDataImpl(FluidSolver *parent, ParticleDataImpl<T> *other)
+ : ParticleDataBase(parent), mpGridSource(NULL), mGridSourceMAC(false)
+{
+ this->mData = other->mData;
+ setName(other->getName());
+}
+
+template<class T> ParticleDataImpl<T>::~ParticleDataImpl()
+{
+}
+
+template<class T> IndexInt ParticleDataImpl<T>::getSizeSlow() const
+{
+ return mData.size();
+}
+template<class T> void ParticleDataImpl<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 ParticleDataImpl<T>::resize(IndexInt s)
+{
+ mData.resize(s);
+}
+template<class T> void ParticleDataImpl<T>::copyValueSlow(IndexInt from, IndexInt to)
+{
+ this->copyValue(from, to);
+}
+template<class T> ParticleDataBase *ParticleDataImpl<T>::clone()
+{
+ ParticleDataImpl<T> *npd = new ParticleDataImpl<T>(getParent(), this);
+ return npd;
+}
+
+template<class T> void ParticleDataImpl<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 ParticleDataImpl<T>::initNewValue(IndexInt idx, Vec3 pos)
+{
+ if (!mpGridSource)
+ mData[idx] = 0;
+ else {
+ mData[idx] = mpGridSource->getInterpolated(pos);
+ }
+}
+// special handling needed for velocities
+template<> void ParticleDataImpl<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);
+ }
+}
+
+template<typename T> void ParticleDataImpl<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")
+ readPdataUni<T>(name, this);
+ else if (ext == ".raw") // raw = uni for now
+ readPdataUni<T>(name, this);
+ else
+ errMsg("particle data '" + name + "' filetype not supported for loading");
+}
+
+template<typename T> void ParticleDataImpl<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")
+ writePdataUni<T>(name, this);
+ else if (ext == ".raw") // raw = uni for now
+ writePdataUni<T>(name, this);
+ else
+ errMsg("particle data '" + name + "' filetype not supported for saving");
+}
+
+// specializations
+
+template<> ParticleDataBase::PdataType ParticleDataImpl<Real>::getType() const
+{
+ return ParticleDataBase::TypeReal;
+}
+template<> ParticleDataBase::PdataType ParticleDataImpl<int>::getType() const
+{
+ return ParticleDataBase::TypeInt;
+}
+template<> ParticleDataBase::PdataType ParticleDataImpl<Vec3>::getType() const
+{
+ return ParticleDataBase::TypeVec3;
+}
+
+// note, we need a flag value for functions such as advection
+// ideally, this value should never be modified
+int ParticleIndexData::flag = 0;
+Vec3 ParticleIndexData::pos = Vec3(0., 0., 0.);
+
+template<class T, class S> struct knPdataAdd : public KernelBase {
+ knPdataAdd(ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other) const
+ {
+ me[idx] += other[idx];
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const ParticleDataImpl<S> &getArg1()
+ {
+ return other;
+ }
+ typedef ParticleDataImpl<S> type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataAdd ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const ParticleDataImpl<S> &other;
+};
+template<class T, class S> struct knPdataSub : public KernelBase {
+ knPdataSub(ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other) const
+ {
+ me[idx] -= other[idx];
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const ParticleDataImpl<S> &getArg1()
+ {
+ return other;
+ }
+ typedef ParticleDataImpl<S> type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataSub ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const ParticleDataImpl<S> &other;
+};
+template<class T, class S> struct knPdataMult : public KernelBase {
+ knPdataMult(ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other) const
+ {
+ me[idx] *= other[idx];
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const ParticleDataImpl<S> &getArg1()
+ {
+ return other;
+ }
+ typedef ParticleDataImpl<S> type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataMult ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const ParticleDataImpl<S> &other;
+};
+template<class T, class S> struct knPdataDiv : public KernelBase {
+ knPdataDiv(ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const ParticleDataImpl<S> &other) const
+ {
+ me[idx] /= other[idx];
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const ParticleDataImpl<S> &getArg1()
+ {
+ return other;
+ }
+ typedef ParticleDataImpl<S> type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataDiv ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const ParticleDataImpl<S> &other;
+};
+template<class T> struct knPdataSafeDiv : public KernelBase {
+ knPdataSafeDiv(ParticleDataImpl<T> &me, const ParticleDataImpl<T> &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const ParticleDataImpl<T> &other) const
+ {
+ me[idx] = safeDivide(me[idx], other[idx]);
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const ParticleDataImpl<T> &getArg1()
+ {
+ return other;
+ }
+ typedef ParticleDataImpl<T> type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataSafeDiv ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const ParticleDataImpl<T> &other;
+};
+
+template<class T, class S> struct knPdataSetScalar : public KernelBase {
+ knPdataSetScalar(ParticleDataImpl<T> &me, const S &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const S &other) const
+ {
+ me[idx] = other;
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const S &getArg1()
+ {
+ return other;
+ }
+ typedef S type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataSetScalar ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const S &other;
+};
+template<class T, class S> struct knPdataAddScalar : public KernelBase {
+ knPdataAddScalar(ParticleDataImpl<T> &me, const S &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const S &other) const
+ {
+ me[idx] += other;
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const S &getArg1()
+ {
+ return other;
+ }
+ typedef S type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataAddScalar ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const S &other;
+};
+template<class T, class S> struct knPdataMultScalar : public KernelBase {
+ knPdataMultScalar(ParticleDataImpl<T> &me, const S &other)
+ : KernelBase(me.size()), me(me), other(other)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const S &other) const
+ {
+ me[idx] *= other;
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const S &getArg1()
+ {
+ return other;
+ }
+ typedef S type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataMultScalar ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const S &other;
+};
+template<class T, class S> struct knPdataScaledAdd : public KernelBase {
+ knPdataScaledAdd(ParticleDataImpl<T> &me, const ParticleDataImpl<T> &other, const S &factor)
+ : KernelBase(me.size()), me(me), other(other), factor(factor)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx,
+ ParticleDataImpl<T> &me,
+ const ParticleDataImpl<T> &other,
+ const S &factor) const
+ {
+ me[idx] += factor * other[idx];
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const ParticleDataImpl<T> &getArg1()
+ {
+ return other;
+ }
+ typedef ParticleDataImpl<T> type1;
+ inline const S &getArg2()
+ {
+ return factor;
+ }
+ typedef S type2;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataScaledAdd ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const ParticleDataImpl<T> &other;
+ const S &factor;
+};
+
+template<class T> struct knPdataClamp : public KernelBase {
+ knPdataClamp(ParticleDataImpl<T> &me, const T vmin, const T vmax)
+ : KernelBase(me.size()), me(me), vmin(vmin), vmax(vmax)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const T vmin, const T vmax) const
+ {
+ me[idx] = clamp(me[idx], vmin, vmax);
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const T &getArg1()
+ {
+ return vmin;
+ }
+ typedef T type1;
+ inline const T &getArg2()
+ {
+ return vmax;
+ }
+ typedef T type2;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataClamp ", 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, vmax);
+ }
+ void run()
+ {
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+ }
+ ParticleDataImpl<T> &me;
+ const T vmin;
+ const T vmax;
+};
+template<class T> struct knPdataClampMin : public KernelBase {
+ knPdataClampMin(ParticleDataImpl<T> &me, const T vmin)
+ : KernelBase(me.size()), me(me), vmin(vmin)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const T vmin) const
+ {
+ me[idx] = std::max(vmin, me[idx]);
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const T &getArg1()
+ {
+ return vmin;
+ }
+ typedef T type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataClampMin ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const T vmin;
+};
+template<class T> struct knPdataClampMax : public KernelBase {
+ knPdataClampMax(ParticleDataImpl<T> &me, const T vmax)
+ : KernelBase(me.size()), me(me), vmax(vmax)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<T> &me, const T vmax) const
+ {
+ me[idx] = std::min(vmax, me[idx]);
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const T &getArg1()
+ {
+ return vmax;
+ }
+ typedef T type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataClampMax ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const T vmax;
+};
+
+struct knPdataClampMinVec3 : public KernelBase {
+ knPdataClampMinVec3(ParticleDataImpl<Vec3> &me, const Real vmin)
+ : KernelBase(me.size()), me(me), vmin(vmin)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<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 ParticleDataImpl<Vec3> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<Vec3> type0;
+ inline const Real &getArg1()
+ {
+ return vmin;
+ }
+ typedef Real type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataClampMinVec3 ", 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);
+ }
+ ParticleDataImpl<Vec3> &me;
+ const Real vmin;
+};
+
+struct knPdataClampMaxVec3 : public KernelBase {
+ knPdataClampMaxVec3(ParticleDataImpl<Vec3> &me, const Real vmax)
+ : KernelBase(me.size()), me(me), vmax(vmax)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, ParticleDataImpl<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 ParticleDataImpl<Vec3> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<Vec3> type0;
+ inline const Real &getArg1()
+ {
+ return vmax;
+ }
+ typedef Real type1;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataClampMaxVec3 ", 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);
+ }
+ ParticleDataImpl<Vec3> &me;
+ const Real vmax;
+};
+
+// python operators
+
+template<typename T>
+ParticleDataImpl<T> &ParticleDataImpl<T>::copyFrom(const ParticleDataImpl<T> &a)
+{
+ assertMsg(a.mData.size() == mData.size(),
+ "different pdata size " << a.mData.size() << " vs " << this->mData.size());
+ mData = a.mData;
+ return *this;
+}
+
+template<typename T> void ParticleDataImpl<T>::setConst(const T &s)
+{
+ knPdataSetScalar<T, T> op(*this, s);
+}
+
+template<typename T>
+void ParticleDataImpl<T>::setConstRange(const 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 knPdataSetScalarIntFlag : public KernelBase {
+ knPdataSetScalarIntFlag(ParticleDataImpl<T> &me,
+ const S &other,
+ const ParticleDataImpl<int> &t,
+ const int itype)
+ : KernelBase(me.size()), me(me), other(other), t(t), itype(itype)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx,
+ ParticleDataImpl<T> &me,
+ const S &other,
+ const ParticleDataImpl<int> &t,
+ const int itype) const
+ {
+ if (t[idx] & itype)
+ me[idx] = other;
+ }
+ inline ParticleDataImpl<T> &getArg0()
+ {
+ return me;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const S &getArg1()
+ {
+ return other;
+ }
+ typedef S type1;
+ inline const ParticleDataImpl<int> &getArg2()
+ {
+ return t;
+ }
+ typedef ParticleDataImpl<int> type2;
+ inline const int &getArg3()
+ {
+ return itype;
+ }
+ typedef int type3;
+ void runMessage()
+ {
+ debMsg("Executing kernel knPdataSetScalarIntFlag ", 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);
+ }
+ ParticleDataImpl<T> &me;
+ const S &other;
+ const ParticleDataImpl<int> &t;
+ const int itype;
+};
+template<typename T>
+void ParticleDataImpl<T>::setConstIntFlag(const T &s,
+ const ParticleDataImpl<int> &t,
+ const int itype)
+{
+ knPdataSetScalarIntFlag<T, T> op(*this, s, t, itype);
+}
+
+template<typename T> void ParticleDataImpl<T>::add(const ParticleDataImpl<T> &a)
+{
+ knPdataAdd<T, T> op(*this, a);
+}
+template<typename T> void ParticleDataImpl<T>::sub(const ParticleDataImpl<T> &a)
+{
+ knPdataSub<T, T> op(*this, a);
+}
+
+template<typename T> void ParticleDataImpl<T>::addConst(const T &s)
+{
+ knPdataAddScalar<T, T> op(*this, s);
+}
+
+template<typename T>
+void ParticleDataImpl<T>::addScaled(const ParticleDataImpl<T> &a, const T &factor)
+{
+ knPdataScaledAdd<T, T> op(*this, a, factor);
+}
+
+template<typename T> void ParticleDataImpl<T>::mult(const ParticleDataImpl<T> &a)
+{
+ knPdataMult<T, T> op(*this, a);
+}
+
+template<typename T> void ParticleDataImpl<T>::safeDiv(const ParticleDataImpl<T> &a)
+{
+ knPdataSafeDiv<T> op(*this, a);
+}
+
+template<typename T> void ParticleDataImpl<T>::multConst(const T &s)
+{
+ knPdataMultScalar<T, T> op(*this, s);
+}
+
+template<typename T> void ParticleDataImpl<T>::clamp(const Real vmin, const Real vmax)
+{
+ knPdataClamp<T> op(*this, vmin, vmax);
+}
+
+template<typename T> void ParticleDataImpl<T>::clampMin(const Real vmin)
+{
+ knPdataClampMin<T> op(*this, vmin);
+}
+template<typename T> void ParticleDataImpl<T>::clampMax(const Real vmax)
+{
+ knPdataClampMax<T> op(*this, vmax);
+}
+
+template<> void ParticleDataImpl<Vec3>::clampMin(const Real vmin)
+{
+ knPdataClampMinVec3 op(*this, vmin);
+}
+template<> void ParticleDataImpl<Vec3>::clampMax(const Real vmax)
+{
+ knPdataClampMaxVec3 op(*this, vmax);
+}
+
+template<typename T> struct KnPtsSum : public KernelBase {
+ KnPtsSum(const ParticleDataImpl<T> &val, const ParticleDataImpl<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 ParticleDataImpl<T> &val,
+ const ParticleDataImpl<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 ParticleDataImpl<T> &getArg0()
+ {
+ return val;
+ }
+ typedef ParticleDataImpl<T> type0;
+ inline const ParticleDataImpl<int> *getArg1()
+ {
+ return t;
+ }
+ typedef ParticleDataImpl<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 ParticleDataImpl<T> &val;
+ const ParticleDataImpl<int> *t;
+ const int itype;
+ T result;
+};
+template<typename T> struct KnPtsSumSquare : public KernelBase {
+ KnPtsSumSquare(const ParticleDataImpl<T> &val) : KernelBase(val.size()), val(val), result(0.)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, const ParticleDataImpl<T> &val, Real &result)
+ {
+ result += normSquare(val[idx]);
+ }
+ inline operator Real()
+ {
+ return result;
+ }
+ inline Real &getRet()
+ {
+ return result;
+ }
+ inline const ParticleDataImpl<T> &getArg0()
+ {
+ return val;
+ }
+ typedef ParticleDataImpl<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 ParticleDataImpl<T> &val;
+ Real result;
+};
+template<typename T> struct KnPtsSumMagnitude : public KernelBase {
+ KnPtsSumMagnitude(const ParticleDataImpl<T> &val) : KernelBase(val.size()), val(val), result(0.)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, const ParticleDataImpl<T> &val, Real &result)
+ {
+ result += norm(val[idx]);
+ }
+ inline operator Real()
+ {
+ return result;
+ }
+ inline Real &getRet()
+ {
+ return result;
+ }
+ inline const ParticleDataImpl<T> &getArg0()
+ {
+ return val;
+ }
+ typedef ParticleDataImpl<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 ParticleDataImpl<T> &val;
+ Real result;
+};
+
+template<typename T>
+T ParticleDataImpl<T>::sum(const ParticleDataImpl<int> *t, const int itype) const
+{
+ return KnPtsSum<T>(*this, t, itype);
+}
+template<typename T> Real ParticleDataImpl<T>::sumSquare() const
+{
+ return KnPtsSumSquare<T>(*this);
+}
+template<typename T> Real ParticleDataImpl<T>::sumMagnitude() const
+{
+ return KnPtsSumMagnitude<T>(*this);
+}
+
+template<typename T>
+
+struct CompPdata_Min : public KernelBase {
+ CompPdata_Min(const ParticleDataImpl<T> &val)
+ : KernelBase(val.size()), val(val), minVal(std::numeric_limits<Real>::max())
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, const ParticleDataImpl<T> &val, Real &minVal)
+ {
+ if (val[idx] < minVal)
+ minVal = val[idx];
+ }
+ inline operator Real()
+ {
+ return minVal;
+ }
+ inline Real &getRet()
+ {
+ return minVal;
+ }
+ inline const ParticleDataImpl<T> &getArg0()
+ {
+ return val;
+ }
+ typedef ParticleDataImpl<T> type0;
+ void runMessage()
+ {
+ debMsg("Executing kernel CompPdata_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);
+ }
+ CompPdata_Min(CompPdata_Min &o, tbb::split)
+ : KernelBase(o), val(o.val), minVal(std::numeric_limits<Real>::max())
+ {
+ }
+ void join(const CompPdata_Min &o)
+ {
+ minVal = min(minVal, o.minVal);
+ }
+ const ParticleDataImpl<T> &val;
+ Real minVal;
+};
+
+template<typename T>
+
+struct CompPdata_Max : public KernelBase {
+ CompPdata_Max(const ParticleDataImpl<T> &val)
+ : KernelBase(val.size()), val(val), maxVal(-std::numeric_limits<Real>::max())
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, const ParticleDataImpl<T> &val, Real &maxVal)
+ {
+ if (val[idx] > maxVal)
+ maxVal = val[idx];
+ }
+ inline operator Real()
+ {
+ return maxVal;
+ }
+ inline Real &getRet()
+ {
+ return maxVal;
+ }
+ inline const ParticleDataImpl<T> &getArg0()
+ {
+ return val;
+ }
+ typedef ParticleDataImpl<T> type0;
+ void runMessage()
+ {
+ debMsg("Executing kernel CompPdata_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);
+ }
+ CompPdata_Max(CompPdata_Max &o, tbb::split)
+ : KernelBase(o), val(o.val), maxVal(-std::numeric_limits<Real>::max())
+ {
+ }
+ void join(const CompPdata_Max &o)
+ {
+ maxVal = max(maxVal, o.maxVal);
+ }
+ const ParticleDataImpl<T> &val;
+ Real maxVal;
+};
+
+template<typename T> Real ParticleDataImpl<T>::getMin() const
+{
+ return CompPdata_Min<T>(*this);
+}
+
+template<typename T> Real ParticleDataImpl<T>::getMaxAbs() const
+{
+ Real amin = CompPdata_Min<T>(*this);
+ Real amax = CompPdata_Max<T>(*this);
+ return max(fabs(amin), fabs(amax));
+}
+
+template<typename T> Real ParticleDataImpl<T>::getMax() const
+{
+ return CompPdata_Max<T>(*this);
+}
+
+template<typename T>
+void ParticleDataImpl<T>::printPdata(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 ParticleDataImpl<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 CompPdata_MinVec3 : public KernelBase {
+ CompPdata_MinVec3(const ParticleDataImpl<Vec3> &val)
+ : KernelBase(val.size()), val(val), minVal(std::numeric_limits<Real>::max())
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, const ParticleDataImpl<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 ParticleDataImpl<Vec3> &getArg0()
+ {
+ return val;
+ }
+ typedef ParticleDataImpl<Vec3> type0;
+ void runMessage()
+ {
+ debMsg("Executing kernel CompPdata_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);
+ }
+ CompPdata_MinVec3(CompPdata_MinVec3 &o, tbb::split)
+ : KernelBase(o), val(o.val), minVal(std::numeric_limits<Real>::max())
+ {
+ }
+ void join(const CompPdata_MinVec3 &o)
+ {
+ minVal = min(minVal, o.minVal);
+ }
+ const ParticleDataImpl<Vec3> &val;
+ Real minVal;
+};
+
+struct CompPdata_MaxVec3 : public KernelBase {
+ CompPdata_MaxVec3(const ParticleDataImpl<Vec3> &val)
+ : KernelBase(val.size()), val(val), maxVal(-std::numeric_limits<Real>::max())
+ {
+ runMessage();
+ run();
+ }
+ inline void op(IndexInt idx, const ParticleDataImpl<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 ParticleDataImpl<Vec3> &getArg0()
+ {
+ return val;
+ }
+ typedef ParticleDataImpl<Vec3> type0;
+ void runMessage()
+ {
+ debMsg("Executing kernel CompPdata_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);
+ }
+ CompPdata_MaxVec3(CompPdata_MaxVec3 &o, tbb::split)
+ : KernelBase(o), val(o.val), maxVal(-std::numeric_limits<Real>::max())
+ {
+ }
+ void join(const CompPdata_MaxVec3 &o)
+ {
+ maxVal = max(maxVal, o.maxVal);
+ }
+ const ParticleDataImpl<Vec3> &val;
+ Real maxVal;
+};
+
+template<> Real ParticleDataImpl<Vec3>::getMin() const
+{
+ return sqrt(CompPdata_MinVec3(*this));
+}
+
+template<> Real ParticleDataImpl<Vec3>::getMaxAbs() const
+{
+ return sqrt(CompPdata_MaxVec3(*this)); // no minimum necessary here
+}
+
+template<> Real ParticleDataImpl<Vec3>::getMax() const
+{
+ return sqrt(CompPdata_MaxVec3(*this));
+}
+
+// explicit instantiation
+template class ParticleDataImpl<int>;
+template class ParticleDataImpl<Real>;
+template class ParticleDataImpl<Vec3>;
+
+} // namespace Manta
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-30 00:04:45 +0300