From e09d0c0d077cff79b55ce32ec5124d5faa73e2e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastia=CC=81n=20Barschkis?= Date: Wed, 25 Nov 2020 23:17:47 +0100 Subject: Fluid: Updated Mantaflow source files This update introduces two improvements from the Mantaflow repository: (1) Improved particle sampling: - Liquid and secondary particles are sampled more predictably. With all parameters being equal, baked particles will be computed at the exact same position during every bake. - Before, this was not guaranteed. (2) Sparse grid caching: - While saving grid data to disk, grids will from now on be saved in a sparse structure whenever possible (e.g. density, flame but not levelsets). - With the sparse optimization grid cells with a value under the 'Empty Space' value (already present in domain settings) will not be cached. - The main benefits of this optimization are: Smaller cache sizes and faster playback of simulation data in the viewport. - This optimization works 'out-of-the-box'. There is no option in the UI to enable it. - For now, only smoke simulation grids will take advantage of this optimization. --- extern/mantaflow/preprocessed/fileio/iovdb.cpp | 120 ++++++++++++++++----- extern/mantaflow/preprocessed/fileio/mantaio.cpp | 20 +++- extern/mantaflow/preprocessed/fileio/mantaio.h | 4 +- extern/mantaflow/preprocessed/gitinfo.h | 2 +- extern/mantaflow/preprocessed/grid.cpp | 29 ++--- extern/mantaflow/preprocessed/grid.h | 30 ++++-- extern/mantaflow/preprocessed/particle.cpp | 10 +- extern/mantaflow/preprocessed/particle.h | 40 ++++--- extern/mantaflow/preprocessed/plugin/flip.cpp | 16 +-- .../preprocessed/plugin/secondaryparticles.cpp | 65 +++++++---- 10 files changed, 239 insertions(+), 97 deletions(-) (limited to 'extern/mantaflow') diff --git a/extern/mantaflow/preprocessed/fileio/iovdb.cpp b/extern/mantaflow/preprocessed/fileio/iovdb.cpp index 287c78f0608..8aca19ade04 100644 --- a/extern/mantaflow/preprocessed/fileio/iovdb.cpp +++ b/extern/mantaflow/preprocessed/fileio/iovdb.cpp @@ -31,6 +31,8 @@ # include "openvdb/openvdb.h" # include # include +# include +# include #endif #define POSITION_NAME "P" @@ -45,15 +47,29 @@ namespace Manta { template void importVDB(typename GridType::Ptr from, Grid *to) { using ValueT = typename GridType::ValueType; - typename GridType::Accessor accessor = from->getAccessor(); - FOR_IJK(*to) - { - openvdb::Coord xyz(i, j, k); - ValueT vdbValue = accessor.getValue(xyz); - T toMantaValue; - convertFrom(vdbValue, &toMantaValue); - to->set(i, j, k, toMantaValue); + // Check if current grid is to be read as a sparse grid, active voxels (only) will be copied + if (to->saveSparse()) { + to->clear(); // Ensure that destination grid is empty before writing + for (typename GridType::ValueOnCIter iter = from->cbeginValueOn(); iter.test(); ++iter) { + ValueT vdbValue = *iter; + openvdb::Coord coord = iter.getCoord(); + T toMantaValue; + convertFrom(vdbValue, &toMantaValue); + to->set(coord.x(), coord.y(), coord.z(), toMantaValue); + } + } + // When importing all grid cells, using a grid accessor is usually faster than a value iterator + else { + typename GridType::Accessor accessor = from->getAccessor(); + FOR_IJK(*to) + { + openvdb::Coord xyz(i, j, k); + ValueT vdbValue = accessor.getValue(xyz); + T toMantaValue; + convertFrom(vdbValue, &toMantaValue); + to->set(i, j, k, toMantaValue); + } } } @@ -173,19 +189,43 @@ static void setGridOptions(typename GridType::Ptr grid, grid->setSaveFloatAsHalf(precision == PRECISION_MINI || precision == PRECISION_HALF); } -template typename GridType::Ptr exportVDB(Grid *from) +template +typename GridType::Ptr exportVDB(Grid *from, float clip, openvdb::FloatGrid::Ptr clipGrid) { using ValueT = typename GridType::ValueType; - typename GridType::Ptr to = GridType::create(); - typename GridType::Accessor accessor = to->getAccessor(); - - FOR_IJK(*from) - { - openvdb::Coord xyz(i, j, k); - T fromMantaValue = (*from)(i, j, k); - ValueT vdbValue; - convertTo(&vdbValue, fromMantaValue); - accessor.setValue(xyz, vdbValue); + typename GridType::Ptr to = GridType::create(ValueT(0)); + + // Copy data from grid by creating a vdb dense structure and then copying that into a vdb grid + // This is the fastest way to copy data for both dense and sparse grids -> if (true) + if (true) { + ValueT *data = (ValueT *)from->getData(); + openvdb::math::CoordBBox bbox( + openvdb::Coord(0), + openvdb::Coord(from->getSizeX() - 1, from->getSizeY() - 1, from->getSizeZ() - 1)); + openvdb::tools::Dense dense(bbox, data); + + // Trick: Set clip value to very small / negative value in order to copy all values of dense + // grids + float tmpClip = (from->saveSparse()) ? clip : -std::numeric_limits::max(); + // Copy from dense to sparse grid structure considering clip value + openvdb::tools::copyFromDense(dense, *to, ValueT(tmpClip)); + + // If present, use clip grid to trim down current vdb grid even more + if (from->saveSparse() && clipGrid && !clipGrid->empty()) { + to = openvdb::tools::clip(*to, *clipGrid); + } + } + // Alternatively, reading all grid cells with an accessor (slightly slower) is possible like this + else { + typename GridType::Accessor accessor = to->getAccessor(); + FOR_IJK(*from) + { + openvdb::Coord xyz(i, j, k); + T fromMantaValue = (*from)(i, j, k); + ValueT vdbValue; + convertTo(&vdbValue, fromMantaValue); + accessor.setValue(xyz, vdbValue); + } } return to; } @@ -346,7 +386,9 @@ int writeObjectsVDB(const string &filename, float worldSize, bool skipDeletedParts, int compression, - int precision) + int precision, + float clip, + const Grid *clipGrid) { openvdb::initialize(); openvdb::io::File file(filename); @@ -357,6 +399,19 @@ int writeObjectsVDB(const string &filename, std::vector pdbBuffer; + // Convert given clip grid to vdb clip grid + openvdb::FloatGrid::Ptr vdbClipGrid = nullptr; + if (clipGrid) { + vdbClipGrid = openvdb::FloatGrid::create(); + Real *data = (Real *)clipGrid->getData(); + openvdb::math::CoordBBox bbox(openvdb::Coord(0), + openvdb::Coord(clipGrid->getSizeX() - 1, + clipGrid->getSizeY() - 1, + clipGrid->getSizeZ() - 1)); + openvdb::tools::Dense dense(bbox, data); + openvdb::tools::copyFromDense(dense, *vdbClipGrid, clip); + } + for (std::vector::iterator iter = objects->begin(); iter != objects->end(); ++iter) { openvdb::GridClass gClass = openvdb::GRID_UNKNOWN; openvdb::GridBase::Ptr vdbGrid; @@ -368,10 +423,14 @@ int writeObjectsVDB(const string &filename, if (GridBase *mantaGrid = dynamic_cast(*iter)) { + if (clipGrid) { + assertMsg(clipGrid->getSize() == mantaGrid->getSize(), + "writeObjectsVDB: Clip grid and exported grid must have the same size"); + } if (mantaGrid->getType() & GridBase::TypeInt) { debMsg("Writing int grid '" << mantaGrid->getName() << "' to vdb file " << filename, 1); Grid *mantaIntGrid = (Grid *)mantaGrid; - vdbGrid = exportVDB(mantaIntGrid); + vdbGrid = exportVDB(mantaIntGrid, clip, vdbClipGrid); gridsVDB.push_back(vdbGrid); } else if (mantaGrid->getType() & GridBase::TypeReal) { @@ -379,7 +438,9 @@ int writeObjectsVDB(const string &filename, gClass = (mantaGrid->getType() & GridBase::TypeLevelset) ? openvdb::GRID_LEVEL_SET : openvdb::GRID_FOG_VOLUME; Grid *mantaRealGrid = (Grid *)mantaGrid; - vdbGrid = exportVDB(mantaRealGrid); + // Only supply clip grid if real grid is not equal to the clip grid + openvdb::FloatGrid::Ptr tmpClipGrid = (mantaRealGrid == clipGrid) ? nullptr : vdbClipGrid; + vdbGrid = exportVDB(mantaRealGrid, clip, tmpClipGrid); gridsVDB.push_back(vdbGrid); } else if (mantaGrid->getType() & GridBase::TypeVec3) { @@ -387,7 +448,7 @@ int writeObjectsVDB(const string &filename, gClass = (mantaGrid->getType() & GridBase::TypeMAC) ? openvdb::GRID_STAGGERED : openvdb::GRID_UNKNOWN; Grid *mantaVec3Grid = (Grid *)mantaGrid; - vdbGrid = exportVDB(mantaVec3Grid); + vdbGrid = exportVDB(mantaVec3Grid, clip, vdbClipGrid); gridsVDB.push_back(vdbGrid); } else { @@ -620,9 +681,12 @@ template void importVDB(openvdb::Int32Grid::Ptr from, G template void importVDB(openvdb::FloatGrid::Ptr from, Grid *to); template void importVDB(openvdb::Vec3SGrid::Ptr from, Grid *to); -template openvdb::Int32Grid::Ptr exportVDB(Grid *from); -template openvdb::FloatGrid::Ptr exportVDB(Grid *from); -template openvdb::Vec3SGrid::Ptr exportVDB(Grid *from); +template openvdb::Int32Grid::Ptr exportVDB( + Grid *from, float clip = 1e-4, openvdb::FloatGrid::Ptr clipGrid = nullptr); +template openvdb::FloatGrid::Ptr exportVDB( + Grid *from, float clip = 1e-4, openvdb::FloatGrid::Ptr clipGrid = nullptr); +template openvdb::Vec3SGrid::Ptr exportVDB( + Grid *from, float clip = 1e-4, openvdb::FloatGrid::Ptr clipGrid = nullptr); openvdb::points::PointDataGrid::Ptr exportVDB(BasicParticleSystem *from, std::vector &fromPData, @@ -652,7 +716,9 @@ int writeObjectsVDB(const string &filename, float worldSize, bool skipDeletedParts, int compression, - int precision) + int precision, + float clip, + const Grid *clipGrid) { errMsg("Cannot save to .vdb file. Mantaflow has not been built with OpenVDB support."); return 0; diff --git a/extern/mantaflow/preprocessed/fileio/mantaio.cpp b/extern/mantaflow/preprocessed/fileio/mantaio.cpp index 3048572261a..fe29890ec11 100644 --- a/extern/mantaflow/preprocessed/fileio/mantaio.cpp +++ b/extern/mantaflow/preprocessed/fileio/mantaio.cpp @@ -83,7 +83,9 @@ int save(const string &name, bool skipDeletedParts = false, int compression = COMPRESSION_ZIP, bool precisionHalf = true, - int precision = PRECISION_HALF) + int precision = PRECISION_HALF, + float clip = 1e-4, + const Grid *clipGrid = nullptr) { if (!precisionHalf) { @@ -102,7 +104,8 @@ int save(const string &name, else if (ext == ".vol") return writeGridsVol(name, &objects); if (ext == ".vdb") - return writeObjectsVDB(name, &objects, worldSize, skipDeletedParts, compression, precision); + return writeObjectsVDB( + name, &objects, worldSize, skipDeletedParts, compression, precision, clip, clipGrid); else if (ext == ".npz") return writeGridsNumpy(name, &objects); else if (ext == ".txt") @@ -129,8 +132,17 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds) int compression = _args.getOpt("compression", 4, COMPRESSION_ZIP, &_lock); bool precisionHalf = _args.getOpt("precisionHalf", 5, true, &_lock); int precision = _args.getOpt("precision", 6, PRECISION_HALF, &_lock); - _retval = toPy( - save(name, objects, worldSize, skipDeletedParts, compression, precisionHalf, precision)); + float clip = _args.getOpt("clip", 7, 1e-4, &_lock); + const Grid *clipGrid = _args.getPtrOpt>("clipGrid", 8, nullptr, &_lock); + _retval = toPy(save(name, + objects, + worldSize, + skipDeletedParts, + compression, + precisionHalf, + precision, + clip, + clipGrid)); _args.check(); } pbFinalizePlugin(parent, "save", !noTiming); diff --git a/extern/mantaflow/preprocessed/fileio/mantaio.h b/extern/mantaflow/preprocessed/fileio/mantaio.h index d49f4b1369d..088d43556e1 100644 --- a/extern/mantaflow/preprocessed/fileio/mantaio.h +++ b/extern/mantaflow/preprocessed/fileio/mantaio.h @@ -75,7 +75,9 @@ int writeObjectsVDB(const std::string &filename, float scale = 1.0, bool skipDeletedParts = false, int compression = COMPRESSION_ZIP, - int precision = PRECISION_HALF); + int precision = PRECISION_HALF, + float clip = 1e-4, + const Grid *clipGrid = nullptr); int readObjectsVDB(const std::string &filename, std::vector *objects, float scale = 1.0); diff --git a/extern/mantaflow/preprocessed/gitinfo.h b/extern/mantaflow/preprocessed/gitinfo.h index c7e62bf1248..d7a69564a3b 100644 --- a/extern/mantaflow/preprocessed/gitinfo.h +++ b/extern/mantaflow/preprocessed/gitinfo.h @@ -1,3 +1,3 @@ -#define MANTA_GIT_VERSION "commit 73990a8a5b876e2b136a646258f714c08c5828da" +#define MANTA_GIT_VERSION "commit bb7cde47b6e04fa62815c70775dc70f02065599f" diff --git a/extern/mantaflow/preprocessed/grid.cpp b/extern/mantaflow/preprocessed/grid.cpp index cf8e4635462..61672129f37 100644 --- a/extern/mantaflow/preprocessed/grid.cpp +++ b/extern/mantaflow/preprocessed/grid.cpp @@ -60,7 +60,7 @@ template<> inline GridBase::GridType typeList() } template -Grid::Grid(FluidSolver *parent, bool show) : GridBase(parent), externalData(false) +Grid::Grid(FluidSolver *parent, bool show, bool sparse) : GridBase(parent), mExternalData(false) { mType = typeList(); mSize = parent->getGridSize(); @@ -70,22 +70,23 @@ Grid::Grid(FluidSolver *parent, bool show) : GridBase(parent), externalData(f mDx = 1.0 / mSize.max(); clear(); setHidden(!show); + +#if OPENVDB == 1 + mSaveSparse = sparse; +#else + if (sparse) + debMsg("Cannot enable sparse save option without OpenVDB", 1); + mSaveSparse = false; +#endif } template -Grid::Grid(FluidSolver *parent, T *data, bool show) - : GridBase(parent), mData(data), externalData(true) +Grid::Grid(FluidSolver *parent, T *data, bool show) : Grid::Grid(parent, show) { - mType = typeList(); - mSize = parent->getGridSize(); - - mStrideZ = parent->is2D() ? 0 : (mSize.x * mSize.y); - mDx = 1.0 / mSize.max(); - - setHidden(!show); + mData = data; } -template Grid::Grid(const Grid &a) : GridBase(a.getParent()), externalData(false) +template Grid::Grid(const Grid &a) : GridBase(a.getParent()), mExternalData(false) { mSize = a.mSize; mType = a.mType; @@ -98,7 +99,7 @@ template Grid::Grid(const Grid &a) : GridBase(a.getParent()), ext template Grid::~Grid() { - if (!externalData) { + if (!mExternalData) { mParent->freeGridPointer(mData); } } @@ -114,8 +115,8 @@ template void Grid::swap(Grid &other) other.getSizeZ() != getSizeZ()) errMsg("Grid::swap(): Grid dimensions mismatch."); - if (externalData || other.externalData) - errMsg("Grid::swap(): Cannot swap if one grid stores externalData."); + if (mExternalData || other.mExternalData) + errMsg("Grid::swap(): Cannot swap if one grid stores mExternalData."); T *dswap = other.mData; other.mData = mData; diff --git a/extern/mantaflow/preprocessed/grid.h b/extern/mantaflow/preprocessed/grid.h index 1f3dc2789ac..cf942a19e9a 100644 --- a/extern/mantaflow/preprocessed/grid.h +++ b/extern/mantaflow/preprocessed/grid.h @@ -402,7 +402,7 @@ class GridBase : public PbClass { template class Grid : public GridBase { public: //! init new grid, values are set to zero - Grid(FluidSolver *parent, bool show = true); + Grid(FluidSolver *parent, bool show = true, bool sparse = false); static int _W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds) { PbClass *obj = Pb::objFromPy(_self); @@ -416,7 +416,8 @@ template class Grid : public GridBase { ArgLocker _lock; FluidSolver *parent = _args.getPtr("parent", 0, &_lock); bool show = _args.getOpt("show", 1, true, &_lock); - obj = new Grid(parent, show); + bool sparse = _args.getOpt("sparse", 2, false, &_lock); + obj = new Grid(parent, show, sparse); obj->registerObject(_self, &_args); _args.check(); } @@ -581,6 +582,16 @@ template class Grid : public GridBase { DEBUG_ONLY(checkIndex(idx)); return mData[idx]; } + //! raw data access + inline T *getData() const + { + return mData; + } + //! query if this grid should be saved as a sparse grid + inline bool saveSparse() + { + return mSaveSparse; + } //! set data inline void set(int i, int j, int k, T &val) @@ -1290,7 +1301,8 @@ template class Grid : public GridBase { protected: T *mData; - bool externalData; // True if mData is managed outside of the Fluidsolver + bool mExternalData; // True if mData is managed outside of the Fluidsolver + bool mSaveSparse; // True if this grid may be cached in a sparse structure public: PbArgs _args; } @@ -1302,7 +1314,8 @@ template class Grid : public GridBase { //! Special function for staggered grids class MACGrid : public Grid { public: - MACGrid(FluidSolver *parent, bool show = true) : Grid(parent, show) + MACGrid(FluidSolver *parent, bool show = true, bool sparse = false) + : Grid(parent, show, sparse) { mType = (GridType)(TypeMAC | TypeVec3); } @@ -1319,7 +1332,8 @@ class MACGrid : public Grid { ArgLocker _lock; FluidSolver *parent = _args.getPtr("parent", 0, &_lock); bool show = _args.getOpt("show", 1, true, &_lock); - obj = new MACGrid(parent, show); + bool sparse = _args.getOpt("sparse", 2, false, &_lock); + obj = new MACGrid(parent, show, sparse); obj->registerObject(_self, &_args); _args.check(); } @@ -1425,7 +1439,8 @@ class MACGrid : public Grid { //! Special functions for FlagGrid class FlagGrid : public Grid { public: - FlagGrid(FluidSolver *parent, int dim = 3, bool show = true) : Grid(parent, show) + FlagGrid(FluidSolver *parent, int dim = 3, bool show = true, bool sparse = false) + : Grid(parent, show, sparse) { mType = (GridType)(TypeFlags | TypeInt); } @@ -1443,7 +1458,8 @@ class FlagGrid : public Grid { FluidSolver *parent = _args.getPtr("parent", 0, &_lock); int dim = _args.getOpt("dim", 1, 3, &_lock); bool show = _args.getOpt("show", 2, true, &_lock); - obj = new FlagGrid(parent, dim, show); + bool sparse = _args.getOpt("sparse", 3, false, &_lock); + obj = new FlagGrid(parent, dim, show, sparse); obj->registerObject(_self, &_args); _args.check(); } diff --git a/extern/mantaflow/preprocessed/particle.cpp b/extern/mantaflow/preprocessed/particle.cpp index 05c561e2c60..ad1c344d307 100644 --- a/extern/mantaflow/preprocessed/particle.cpp +++ b/extern/mantaflow/preprocessed/particle.cpp @@ -28,9 +28,15 @@ using namespace std; namespace Manta { -ParticleBase::ParticleBase(FluidSolver *parent) - : PbClass(parent), mMaxParticles(0), mAllowCompress(true), mFreePdata(false) +int ParticleBase::globalSeed = 9832; + +ParticleBase::ParticleBase(FluidSolver *parent, int fixedSeed) + : PbClass(parent), mMaxParticles(0), mAllowCompress(true), mFreePdata(false), mSeed(fixedSeed) { + // use global random seed if none is given + if (fixedSeed == -1) { + mSeed = globalSeed; + } } ParticleBase::~ParticleBase() diff --git a/extern/mantaflow/preprocessed/particle.h b/extern/mantaflow/preprocessed/particle.h index 7fed33c737d..7fcc7e5ca32 100644 --- a/extern/mantaflow/preprocessed/particle.h +++ b/extern/mantaflow/preprocessed/particle.h @@ -47,7 +47,7 @@ class ParticleBase : public PbClass { PINVALID = (1 << 30), // unused }; - ParticleBase(FluidSolver *parent); + ParticleBase(FluidSolver *parent, int fixedSeed = -1); static int _W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds) { PbClass *obj = Pb::objFromPy(_self); @@ -60,7 +60,8 @@ class ParticleBase : public PbClass { { ArgLocker _lock; FluidSolver *parent = _args.getPtr("parent", 0, &_lock); - obj = new ParticleBase(parent); + int fixedSeed = _args.getOpt("fixedSeed", 1, -1, &_lock); + obj = new ParticleBase(parent, fixedSeed); obj->registerObject(_self, &_args); _args.check(); } @@ -111,6 +112,11 @@ class ParticleBase : public PbClass { return; } + inline int getSeed() + { + return mSeed; + } + //! particle data functions //! create a particle data object @@ -192,9 +198,13 @@ class ParticleBase : public PbClass { //! per particle) std::vector *> mPdataReal; std::vector *> mPdataVec3; - std::vector *> - mPdataInt; //! indicate that pdata of this particle system is copied, and needs to be freed + std::vector *> mPdataInt; + //! indicate that pdata of this particle system is copied, and needs to be freed bool mFreePdata; + + //! custom seed for particle systems, used by plugins + int mSeed; //! fix global random seed storage, used mainly by functions in this class + static int globalSeed; public: PbArgs _args; } @@ -2285,15 +2295,14 @@ void ParticleSystem::advectInGrid(const FlagGrid &flags, } template struct KnProjectParticles : public KernelBase { - KnProjectParticles(ParticleSystem &part, Grid &gradient) - : KernelBase(part.size()), part(part), gradient(gradient) + KnProjectParticles(ParticleSystem &part, Grid &gradient, RandomStream &rand) + : KernelBase(part.size()), part(part), gradient(gradient), rand(rand) { runMessage(); run(); } - inline void op(IndexInt idx, ParticleSystem &part, Grid &gradient) + inline void op(IndexInt idx, ParticleSystem &part, Grid &gradient, RandomStream &rand) { - static RandomStream rand(3123984); const double jlen = 0.1; if (part.isActive(idx)) { @@ -2321,6 +2330,11 @@ template struct KnProjectParticles : public KernelBase { return gradient; } typedef Grid type1; + inline RandomStream &getArg2() + { + return rand; + } + typedef RandomStream type2; void runMessage() { debMsg("Executing kernel KnProjectParticles ", 3); @@ -2332,15 +2346,17 @@ template struct KnProjectParticles : public KernelBase { { const IndexInt _sz = size; for (IndexInt i = 0; i < _sz; i++) - op(i, part, gradient); + op(i, part, gradient, rand); } ParticleSystem ∂ Grid &gradient; + RandomStream &rand; }; template void ParticleSystem::projectOutside(Grid &gradient) { - KnProjectParticles(*this, gradient); + RandomStream rand(globalSeed); + KnProjectParticles(*this, gradient, rand); } template struct KnProjectOutOfBnd : public KernelBase { @@ -2531,14 +2547,14 @@ template void ParticleSystem::insertBufferedParticles() int insertFlag; Vec3 insertPos; - static RandomStream mRand(9832); + RandomStream rand(globalSeed); for (IndexInt i = 0; i < numNewParts; ++i) { // get random index in newBuffer vector // we are inserting particles randomly so that they are sampled uniformly in the fluid region // otherwise, regions of fluid can remain completely empty once mData.size() == maxParticles is // reached. - int randIndex = floor(mRand.getReal() * mNewBufferPos.size()); + int randIndex = floor(rand.getReal() * mNewBufferPos.size()); // get elements from new buffers with random index std::swap(mNewBufferPos[randIndex], mNewBufferPos.back()); diff --git a/extern/mantaflow/preprocessed/plugin/flip.cpp b/extern/mantaflow/preprocessed/plugin/flip.cpp index 1acdac1c094..eaa1ebe45d5 100644 --- a/extern/mantaflow/preprocessed/plugin/flip.cpp +++ b/extern/mantaflow/preprocessed/plugin/flip.cpp @@ -41,7 +41,7 @@ void sampleFlagsWithParticles(const FlagGrid &flags, const bool is3D = flags.is3D(); const Real jlen = randomness / discretization; const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization); - RandomStream mRand(9832); + RandomStream rand(parts.getSeed()); FOR_IJK_BND(flags, 0) { @@ -53,7 +53,7 @@ void sampleFlagsWithParticles(const FlagGrid &flags, for (int dj = 0; dj < discretization; dj++) for (int di = 0; di < discretization; di++) { Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk); - subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3()); + subpos += jlen * (Vec3(1, 1, 1) - 2.0 * rand.getVec3()); if (!is3D) subpos[2] = 0.5; parts.addBuffered(subpos); @@ -113,7 +113,7 @@ void sampleLevelsetWithParticles(const LevelsetGrid &phi, const bool is3D = phi.is3D(); const Real jlen = randomness / discretization; const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization); - RandomStream mRand(9832); + RandomStream rand(parts.getSeed()); if (reset) { parts.clear(); @@ -132,7 +132,7 @@ void sampleLevelsetWithParticles(const LevelsetGrid &phi, for (int dj = 0; dj < discretization; dj++) for (int di = 0; di < discretization; di++) { Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk); - subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3()); + subpos += jlen * (Vec3(1, 1, 1) - 2.0 * rand.getVec3()); if (!is3D) subpos[2] = 0.5; if (phi.getInterpolated(subpos) > 0.) @@ -205,7 +205,7 @@ void sampleShapeWithParticles(const Shape &shape, const bool is3D = flags.is3D(); const Real jlen = randomness / discretization; const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization); - RandomStream mRand(9832); + RandomStream rand(parts.getSeed()); if (reset) { parts.clear(); @@ -223,7 +223,7 @@ void sampleShapeWithParticles(const Shape &shape, for (int dj = 0; dj < discretization; dj++) for (int di = 0; di < discretization; di++) { Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk); - subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3()); + subpos += jlen * (Vec3(1, 1, 1) - 2.0 * rand.getVec3()); if (!is3D) subpos[2] = 0.5; if (exclude && exclude->getInterpolated(subpos) <= 0.) @@ -576,7 +576,7 @@ void adjustNumber(BasicParticleSystem &parts, } // seed new particles - RandomStream mRand(9832); + RandomStream rand(parts.getSeed()); FOR_IJK(tmp) { int cnt = tmp(i, j, k); @@ -593,7 +593,7 @@ void adjustNumber(BasicParticleSystem &parts, if (flags.isFluid(i, j, k) && cnt < minParticles) { for (int m = cnt; m < minParticles; m++) { - Vec3 pos = Vec3(i, j, k) + mRand.getVec3(); + Vec3 pos = Vec3(i, j, k) + rand.getVec3(); // Vec3 pos (i + 0.5, j + 0.5, k + 0.5); // cell center parts.addBuffered(pos); } diff --git a/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp b/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp index 82afe2d7ab2..8ebc239e8fc 100644 --- a/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp +++ b/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp @@ -479,7 +479,8 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { const Real k_ta, const Real k_wc, const Real dt, - const int itype = FlagGrid::TypeFluid) + const int itype, + RandomStream &rand) : KernelBase(&flags, 0), flags(flags), v(v), @@ -497,7 +498,8 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { k_ta(k_ta), k_wc(k_wc), dt(dt), - itype(itype) + itype(itype), + rand(rand) { runMessage(); run(); @@ -521,13 +523,13 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { const Real k_ta, const Real k_wc, const Real dt, - const int itype = FlagGrid::TypeFluid) + const int itype, + RandomStream &rand) { if (!(flags(i, j, k) & itype)) return; - static RandomStream mRand(9832); Real radius = 0.25; // diameter=0.5 => sampling with two cylinders in each dimension since cell size=1 for (Real x = i - radius; x <= i + radius; x += 2 * radius) { @@ -549,9 +551,9 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { cross(e1, dir)); // perpendicular to dir and e1, so e1 and e1 create reference plane for (int di = 0; di < n; di++) { - const Real r = radius * sqrt(mRand.getReal()); // distance to cylinder axis - const Real theta = mRand.getReal() * Real(2) * M_PI; // azimuth - const Real h = mRand.getReal() * norm(dt * vi); // distance to reference plane + const Real r = radius * sqrt(rand.getReal()); // distance to cylinder axis + const Real theta = rand.getReal() * Real(2) * M_PI; // azimuth + const Real h = rand.getReal() * norm(dt * vi); // distance to reference plane Vec3 xd = xi + r * cos(theta) * e1 + r * sin(theta) * e2 + h * getNormalized(vi); if (!flags.is3D()) xd.z = 0; @@ -561,7 +563,7 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { vi; // init velocity of new particle Real temp = (KE + TA + WC) / 3; l_sec[l_sec.size() - 1] = ((lMax - lMin) * temp) + lMin + - mRand.getReal() * 0.1; // init lifetime of new particle + rand.getReal() * 0.1; // init lifetime of new particle // init type of new particle if (neighborRatio(i, j, k) < c_s) { @@ -663,6 +665,11 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { return itype; } typedef int type16; + inline RandomStream &getArg17() + { + return rand; + } + typedef RandomStream type17; void runMessage() { debMsg("Executing kernel knFlipSampleSecondaryParticlesMoreCylinders ", 3); @@ -696,7 +703,8 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { k_ta, k_wc, dt, - itype); + itype, + rand); } const FlagGrid &flags; const MACGrid &v; @@ -715,6 +723,7 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase { const Real k_wc; const Real dt; const int itype; + RandomStream &rand; }; // adds secondary particles to &pts_sec for every fluid cell in &flags according to the potential @@ -738,7 +747,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase { const Real k_ta, const Real k_wc, const Real dt, - const int itype = FlagGrid::TypeFluid) + const int itype, + RandomStream &rand) : KernelBase(&flags, 0), flags(flags), v(v), @@ -756,7 +766,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase { k_ta(k_ta), k_wc(k_wc), dt(dt), - itype(itype) + itype(itype), + rand(rand) { runMessage(); run(); @@ -780,7 +791,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase { const Real k_ta, const Real k_wc, const Real dt, - const int itype = FlagGrid::TypeFluid) + const int itype, + RandomStream &rand) { if (!(flags(i, j, k) & itype)) @@ -793,9 +805,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase { const int n = KE * (k_ta * TA + k_wc * WC) * dt; // number of secondary particles if (n == 0) return; - static RandomStream mRand(9832); - Vec3 xi = Vec3(i, j, k) + mRand.getVec3(); // randomized offset uniform in cell + Vec3 xi = Vec3(i, j, k) + rand.getVec3(); // randomized offset uniform in cell Vec3 vi = v.getInterpolated(xi); Vec3 dir = dt * vi; // direction of movement of current particle Vec3 e1 = getNormalized(Vec3(dir.z, 0, -dir.x)); // perpendicular to dir @@ -803,9 +814,9 @@ struct knFlipSampleSecondaryParticles : public KernelBase { cross(e1, dir)); // perpendicular to dir and e1, so e1 and e1 create reference plane for (int di = 0; di < n; di++) { - const Real r = Real(0.5) * sqrt(mRand.getReal()); // distance to cylinder axis - const Real theta = mRand.getReal() * Real(2) * M_PI; // azimuth - const Real h = mRand.getReal() * norm(dt * vi); // distance to reference plane + const Real r = Real(0.5) * sqrt(rand.getReal()); // distance to cylinder axis + const Real theta = rand.getReal() * Real(2) * M_PI; // azimuth + const Real h = rand.getReal() * norm(dt * vi); // distance to reference plane Vec3 xd = xi + r * cos(theta) * e1 + r * sin(theta) * e2 + h * getNormalized(vi); if (!flags.is3D()) xd.z = 0; @@ -815,7 +826,7 @@ struct knFlipSampleSecondaryParticles : public KernelBase { vi; // init velocity of new particle Real temp = (KE + TA + WC) / 3; l_sec[l_sec.size() - 1] = ((lMax - lMin) * temp) + lMin + - mRand.getReal() * 0.1; // init lifetime of new particle + rand.getReal() * 0.1; // init lifetime of new particle // init type of new particle if (neighborRatio(i, j, k) < c_s) { @@ -914,6 +925,11 @@ struct knFlipSampleSecondaryParticles : public KernelBase { return itype; } typedef int type16; + inline RandomStream &getArg17() + { + return rand; + } + typedef RandomStream type17; void runMessage() { debMsg("Executing kernel knFlipSampleSecondaryParticles ", 3); @@ -947,7 +963,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase { k_ta, k_wc, dt, - itype); + itype, + rand); } const FlagGrid &flags; const MACGrid &v; @@ -966,6 +983,7 @@ struct knFlipSampleSecondaryParticles : public KernelBase { const Real k_wc; const Real dt; const int itype; + RandomStream &rand; }; void flipSampleSecondaryParticles(const std::string mode, @@ -992,6 +1010,9 @@ void flipSampleSecondaryParticles(const std::string mode, if (dt <= 0) timestep = flags.getParent()->getDt(); + /* Every particle needs to get a different random offset. */ + RandomStream rand(pts_sec.getSeed()); + if (mode == "single") { knFlipSampleSecondaryParticles(flags, v, @@ -1009,7 +1030,8 @@ void flipSampleSecondaryParticles(const std::string mode, k_ta, k_wc, timestep, - itype); + itype, + rand); } else if (mode == "multiple") { knFlipSampleSecondaryParticlesMoreCylinders(flags, @@ -1028,7 +1050,8 @@ void flipSampleSecondaryParticles(const std::string mode, k_ta, k_wc, timestep, - itype); + itype, + rand); } else { throw std::invalid_argument("Unknown mode: use \"single\" or \"multiple\" instead!"); -- cgit v1.2.3