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
diff options
context:
space:
mode:
authorSebastián Barschkis <sebbas@sebbas.org>2020-12-23 15:54:47 +0300
committerSebastián Barschkis <sebbas@sebbas.org>2020-12-23 17:48:38 +0300
commit635694c0ff8fc5c9828bf920ecb81bb9bf792a82 (patch)
treebf2c3166130710a5a460753f5e6eedd9348c08ea
parent5cfda8e7f74da959a5c1081d45e7608bef95191d (diff)
Fluid: Added new viscosity solver
Mainly updated the Mantaflow version. It includes the new viscosity solver plugin based on the method from 'Accurate Viscous Free Surfaces for Buckling, Coiling, and Rotating Liquids' (Batty & Bridson). In the UI, this update adds a new 'Viscosity' section to the fluid modifier UI (liquid domains only). For now, there is a single 'strength' value to control the viscosity of liquids.
Diffstat
-rw-r--r--extern/mantaflow/CMakeLists.txt1
-rw-r--r--extern/mantaflow/helper/util/rcmatrix.h2
-rw-r--r--extern/mantaflow/preprocessed/conjugategrad.cpp94
-rw-r--r--extern/mantaflow/preprocessed/conjugategrad.h466
-rw-r--r--extern/mantaflow/preprocessed/general.h35
-rw-r--r--extern/mantaflow/preprocessed/gitinfo.h2
-rw-r--r--extern/mantaflow/preprocessed/plugin/extforces.cpp25
-rw-r--r--extern/mantaflow/preprocessed/plugin/pressure.cpp14
-rw-r--r--extern/mantaflow/preprocessed/plugin/viscosity.cpp1430
-rw-r--r--extern/mantaflow/preprocessed/plugin/vortexplugins.cpp4
-rw-r--r--extern/mantaflow/preprocessed/plugin/waves.cpp13
-rw-r--r--extern/mantaflow/preprocessed/registration.cpp2
-rw-r--r--intern/mantaflow/intern/MANTA_main.cpp27
-rw-r--r--intern/mantaflow/intern/MANTA_main.h2
-rw-r--r--intern/mantaflow/intern/strings/fluid_script.h11
-rw-r--r--intern/mantaflow/intern/strings/liquid_script.h17
-rw-r--r--release/scripts/startup/bl_ui/properties_physics_fluid.py41
-rw-r--r--source/blender/blenkernel/BKE_blender_version.h2
-rw-r--r--source/blender/blenkernel/intern/fluid.c3
-rw-r--r--source/blender/blenloader/intern/versioning_290.c40
-rw-r--r--source/blender/blenloader/intern/versioning_userdef.c9
-rw-r--r--source/blender/makesdna/DNA_fluid_defaults.h1
-rw-r--r--source/blender/makesdna/DNA_fluid_types.h63
-rw-r--r--source/blender/makesrna/intern/rna_fluid.c16
24 files changed, 2148 insertions, 172 deletions
diff --git a/extern/mantaflow/CMakeLists.txt b/extern/mantaflow/CMakeLists.txt
index ccf272650e3..82bf95a9742 100644
--- a/extern/mantaflow/CMakeLists.txt
+++ b/extern/mantaflow/CMakeLists.txt
@@ -200,6 +200,7 @@ set(SRC
${MANTA_PP}/plugin/ptsplugins.cpp
${MANTA_PP}/plugin/secondaryparticles.cpp
${MANTA_PP}/plugin/surfaceturbulence.cpp
+ ${MANTA_PP}/plugin/viscosity.cpp
${MANTA_PP}/plugin/vortexplugins.cpp
${MANTA_PP}/plugin/waveletturbulence.cpp
${MANTA_PP}/plugin/waves.cpp
diff --git a/extern/mantaflow/helper/util/rcmatrix.h b/extern/mantaflow/helper/util/rcmatrix.h
index f1f0efe6416..330fd1f64f7 100644
--- a/extern/mantaflow/helper/util/rcmatrix.h
+++ b/extern/mantaflow/helper/util/rcmatrix.h
@@ -1035,7 +1035,7 @@ template<class N, class T> struct RCFixedMatrix {
typedef RCMatrix<int, Real> Matrix;
typedef RCFixedMatrix<int, Real> FixedMatrix;
-} // namespace Manta
+}
#undef parallel_for
#undef parallel_end
diff --git a/extern/mantaflow/preprocessed/conjugategrad.cpp b/extern/mantaflow/preprocessed/conjugategrad.cpp
index da82a412a97..bdcceb29520 100644
--- a/extern/mantaflow/preprocessed/conjugategrad.cpp
+++ b/extern/mantaflow/preprocessed/conjugategrad.cpp
@@ -397,7 +397,7 @@ struct UpdateSearchVec : public KernelBase {
};
//*****************************************************************************
-// CG class
+// CG class
template<class APPLYMAT>
GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
@@ -406,10 +406,8 @@ GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
Grid<Real> &search,
const FlagGrid &flags,
Grid<Real> &tmp,
- Grid<Real> *pA0,
- Grid<Real> *pAi,
- Grid<Real> *pAj,
- Grid<Real> *pAk)
+ std::vector<Grid<Real> *> matrixAVec,
+ std::vector<Grid<Real> *> rhsVec)
: GridCgInterface(),
mInited(false),
mIterations(0),
@@ -419,10 +417,8 @@ GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
mSearch(search),
mFlags(flags),
mTmp(tmp),
- mpA0(pA0),
- mpAi(pAi),
- mpAj(pAj),
- mpAk(pAk),
+ mMatrixA(matrixAVec),
+ mVecRhs(rhsVec),
mPcMethod(PC_None),
mpPCA0(nullptr),
mpPCAi(nullptr),
@@ -445,19 +441,37 @@ template<class APPLYMAT> void GridCg<APPLYMAT>::doInit()
if (mPcMethod == PC_ICP) {
assertMsg(mDst.is3D(), "ICP only supports 3D grids so far");
- InitPreconditionIncompCholesky(
- mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
- ApplyPreconditionIncompCholesky(
- mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
+ InitPreconditionIncompCholesky(mFlags,
+ *mpPCA0,
+ *mpPCAi,
+ *mpPCAj,
+ *mpPCAk,
+ *mMatrixA[0],
+ *mMatrixA[1],
+ *mMatrixA[2],
+ *mMatrixA[3]);
+ ApplyPreconditionIncompCholesky(mTmp,
+ mResidual,
+ mFlags,
+ *mpPCA0,
+ *mpPCAi,
+ *mpPCAj,
+ *mpPCAk,
+ *mMatrixA[0],
+ *mMatrixA[1],
+ *mMatrixA[2],
+ *mMatrixA[3]);
}
else if (mPcMethod == PC_mICP) {
assertMsg(mDst.is3D(), "mICP only supports 3D grids so far");
- InitPreconditionModifiedIncompCholesky2(mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
+ InitPreconditionModifiedIncompCholesky2(
+ mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]);
ApplyPreconditionModifiedIncompCholesky2(
- mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
+ mTmp, mResidual, mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]);
}
else if (mPcMethod == PC_MGP) {
- InitPreconditionMultigrid(mMG, *mpA0, *mpAi, *mpAj, *mpAk, mAccuracy);
+ InitPreconditionMultigrid(
+ mMG, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3], mAccuracy);
ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
}
else {
@@ -465,7 +479,6 @@ template<class APPLYMAT> void GridCg<APPLYMAT>::doInit()
}
mSearch.copyFrom(mTmp);
-
mSigma = GridDotProduct(mTmp, mResidual);
}
@@ -480,7 +493,7 @@ template<class APPLYMAT> bool GridCg<APPLYMAT>::iterate()
// this could reinterpret the mpA pointers (not so clean right now)
// tmp = applyMat(search)
- APPLYMAT(mFlags, mTmp, mSearch, *mpA0, *mpAi, *mpAj, *mpAk);
+ APPLYMAT(mFlags, mTmp, mSearch, mMatrixA, mVecRhs);
// alpha = sigma/dot(tmp, search)
Real dp = GridDotProduct(mTmp, mSearch);
@@ -492,11 +505,20 @@ template<class APPLYMAT> bool GridCg<APPLYMAT>::iterate()
gridScaledAdd<Real, Real>(mResidual, mTmp, -alpha); // residual += tmp * -alpha
if (mPcMethod == PC_ICP)
- ApplyPreconditionIncompCholesky(
- mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
+ ApplyPreconditionIncompCholesky(mTmp,
+ mResidual,
+ mFlags,
+ *mpPCA0,
+ *mpPCAi,
+ *mpPCAj,
+ *mpPCAk,
+ *mMatrixA[0],
+ *mMatrixA[1],
+ *mMatrixA[2],
+ *mMatrixA[3]);
else if (mPcMethod == PC_mICP)
ApplyPreconditionModifiedIncompCholesky2(
- mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
+ mTmp, mResidual, mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]);
else if (mPcMethod == PC_MGP)
ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
else
@@ -584,13 +606,15 @@ void GridCg<APPLYMAT>::setMGPreconditioner(PreconditionType method, GridMg *MG)
assertMsg(method == PC_MGP, "GridCg<APPLYMAT>::setMGPreconditioner: Invalid method specified.");
mPcMethod = method;
-
mMG = MG;
}
// explicit instantiation
template class GridCg<ApplyMatrix>;
template class GridCg<ApplyMatrix2D>;
+template class GridCg<ApplyMatrixViscosityU>;
+template class GridCg<ApplyMatrixViscosityV>;
+template class GridCg<ApplyMatrixViscosityW>;
//*****************************************************************************
// diffusion for real and vec grids, e.g. for viscosity
@@ -638,10 +662,15 @@ void cgSolveDiffusion(const FlagGrid &flags,
if (grid.getType() & GridBase::TypeReal) {
Grid<Real> &u = ((Grid<Real> &)grid);
rhs.copyFrom(u);
- if (flags.is3D())
- gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
- else
- gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+ vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
+
+ if (flags.is3D()) {
+ matA.push_back(&Ak);
+ gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, matA);
+ }
+ else {
+ gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, matA);
+ }
gcg->setAccuracy(cgAccuracy);
gcg->solve(maxIter);
@@ -653,12 +682,17 @@ void cgSolveDiffusion(const FlagGrid &flags,
else if ((grid.getType() & GridBase::TypeVec3) || (grid.getType() & GridBase::TypeMAC)) {
Grid<Vec3> &vec = ((Grid<Vec3> &)grid);
Grid<Real> u(parent);
+ vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
// core solve is same as for a regular real grid
- if (flags.is3D())
- gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
- else
- gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+ if (flags.is3D()) {
+ matA.push_back(&Ak);
+ gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, matA);
+ }
+ else {
+ gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, matA);
+ }
+
gcg->setAccuracy(cgAccuracy);
// diffuse every component separately
diff --git a/extern/mantaflow/preprocessed/conjugategrad.h b/extern/mantaflow/preprocessed/conjugategrad.h
index 58ccff28179..4974aa6a4d6 100644
--- a/extern/mantaflow/preprocessed/conjugategrad.h
+++ b/extern/mantaflow/preprocessed/conjugategrad.h
@@ -78,13 +78,9 @@ template<class APPLYMAT> class GridCg : public GridCgInterface {
Grid<Real> &search,
const FlagGrid &flags,
Grid<Real> &tmp,
- Grid<Real> *A0,
- Grid<Real> *pAi,
- Grid<Real> *pAj,
- Grid<Real> *pAk);
- ~GridCg()
- {
- }
+ std::vector<Grid<Real> *> matrixAVec,
+ std::vector<Grid<Real> *> rhsVec = {});
+ ~GridCg(){};
void doInit();
bool iterate();
@@ -133,7 +129,10 @@ template<class APPLYMAT> class GridCg : public GridCgInterface {
const FlagGrid &mFlags;
Grid<Real> &mTmp;
- Grid<Real> *mpA0, *mpAi, *mpAj, *mpAk;
+ //! shape of A matrix defined here (e.g. diagonal, positive neighbor cells, etc)
+ std::vector<Grid<Real> *> mMatrixA;
+ //! shape of rhs vector defined here (e.g. 1 rhs for regular fluids solve, 3 rhs for viscosity)
+ std::vector<Grid<Real> *> mVecRhs;
PreconditionType mPcMethod;
//! preconditioning grids
@@ -154,11 +153,9 @@ struct ApplyMatrix : public KernelBase {
ApplyMatrix(const FlagGrid &flags,
Grid<Real> &dst,
const Grid<Real> &src,
- Grid<Real> &A0,
- Grid<Real> &Ai,
- Grid<Real> &Aj,
- Grid<Real> &Ak)
- : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak)
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs)
+ : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
{
runMessage();
run();
@@ -167,11 +164,18 @@ struct ApplyMatrix : public KernelBase {
const FlagGrid &flags,
Grid<Real> &dst,
const Grid<Real> &src,
- Grid<Real> &A0,
- Grid<Real> &Ai,
- Grid<Real> &Aj,
- Grid<Real> &Ak) const
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs) const
{
+ unusedParameter(vecRhs); // Not needed in this matrix application
+
+ if (matrixA.size() != 4)
+ errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+ Grid<Real> &A0 = *matrixA[0];
+ Grid<Real> &Ai = *matrixA[1];
+ Grid<Real> &Aj = *matrixA[2];
+ Grid<Real> &Ak = *matrixA[3];
+
if (!flags.isFluid(idx)) {
dst[idx] = src[idx];
return;
@@ -196,26 +200,16 @@ struct ApplyMatrix : public KernelBase {
return src;
}
typedef Grid<Real> type2;
- inline Grid<Real> &getArg3()
- {
- return A0;
- }
- typedef Grid<Real> type3;
- inline Grid<Real> &getArg4()
- {
- return Ai;
- }
- typedef Grid<Real> type4;
- inline Grid<Real> &getArg5()
+ inline const std::vector<Grid<Real> *> &getArg3()
{
- return Aj;
+ return matrixA;
}
- typedef Grid<Real> type5;
- inline Grid<Real> &getArg6()
+ typedef std::vector<Grid<Real> *> type3;
+ inline const std::vector<Grid<Real> *> &getArg4()
{
- return Ak;
+ return vecRhs;
}
- typedef Grid<Real> type6;
+ typedef std::vector<Grid<Real> *> type4;
void runMessage()
{
debMsg("Executing kernel ApplyMatrix ", 3);
@@ -226,7 +220,7 @@ struct ApplyMatrix : public KernelBase {
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
- op(idx, flags, dst, src, A0, Ai, Aj, Ak);
+ op(idx, flags, dst, src, matrixA, vecRhs);
}
void run()
{
@@ -235,10 +229,8 @@ struct ApplyMatrix : public KernelBase {
const FlagGrid &flags;
Grid<Real> &dst;
const Grid<Real> &src;
- Grid<Real> &A0;
- Grid<Real> &Ai;
- Grid<Real> &Aj;
- Grid<Real> &Ak;
+ const std::vector<Grid<Real> *> matrixA;
+ const std::vector<Grid<Real> *> vecRhs;
};
//! Kernel: Apply symmetric stored Matrix. 2D version
@@ -247,11 +239,9 @@ struct ApplyMatrix2D : public KernelBase {
ApplyMatrix2D(const FlagGrid &flags,
Grid<Real> &dst,
const Grid<Real> &src,
- Grid<Real> &A0,
- Grid<Real> &Ai,
- Grid<Real> &Aj,
- Grid<Real> &Ak)
- : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak)
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs)
+ : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
{
runMessage();
run();
@@ -260,12 +250,16 @@ struct ApplyMatrix2D : public KernelBase {
const FlagGrid &flags,
Grid<Real> &dst,
const Grid<Real> &src,
- Grid<Real> &A0,
- Grid<Real> &Ai,
- Grid<Real> &Aj,
- Grid<Real> &Ak) const
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs) const
{
- unusedParameter(Ak); // only there for parameter compatibility with ApplyMatrix
+ unusedParameter(vecRhs); // Not needed in this matrix application
+
+ if (matrixA.size() != 3)
+ errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+ Grid<Real> &A0 = *matrixA[0];
+ Grid<Real> &Ai = *matrixA[1];
+ Grid<Real> &Aj = *matrixA[2];
if (!flags.isFluid(idx)) {
dst[idx] = src[idx];
@@ -290,51 +284,387 @@ struct ApplyMatrix2D : public KernelBase {
return src;
}
typedef Grid<Real> type2;
- inline Grid<Real> &getArg3()
+ inline const std::vector<Grid<Real> *> &getArg3()
{
- return A0;
+ return matrixA;
}
- typedef Grid<Real> type3;
- inline Grid<Real> &getArg4()
+ typedef std::vector<Grid<Real> *> type3;
+ inline const std::vector<Grid<Real> *> &getArg4()
{
- return Ai;
+ return vecRhs;
}
- typedef Grid<Real> type4;
- inline Grid<Real> &getArg5()
+ typedef std::vector<Grid<Real> *> type4;
+ void runMessage()
{
- return Aj;
+ debMsg("Executing kernel ApplyMatrix2D ", 3);
+ debMsg("Kernel range"
+ << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+ 4);
+ };
+ void operator()(const tbb::blocked_range<IndexInt> &__r) const
+ {
+ for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+ op(idx, flags, dst, src, matrixA, vecRhs);
}
- typedef Grid<Real> type5;
- inline Grid<Real> &getArg6()
+ void run()
{
- return Ak;
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+ }
+ const FlagGrid &flags;
+ Grid<Real> &dst;
+ const Grid<Real> &src;
+ const std::vector<Grid<Real> *> matrixA;
+ const std::vector<Grid<Real> *> vecRhs;
+};
+
+struct ApplyMatrixViscosityU : public KernelBase {
+ ApplyMatrixViscosityU(const FlagGrid &flags,
+ Grid<Real> &dst,
+ const Grid<Real> &src,
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs)
+ : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
+ {
+ runMessage();
+ run();
}
- typedef Grid<Real> type6;
+ inline void op(int i,
+ int j,
+ int k,
+ const FlagGrid &flags,
+ Grid<Real> &dst,
+ const Grid<Real> &src,
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs) const
+ {
+ if (matrixA.size() != 15)
+ errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+ Grid<Real> &A0 = *matrixA[0];
+ Grid<Real> &Aplusi = *matrixA[1];
+ Grid<Real> &Aplusj = *matrixA[2];
+ Grid<Real> &Aplusk = *matrixA[3];
+ Grid<Real> &Aminusi = *matrixA[4];
+ Grid<Real> &Aminusj = *matrixA[5];
+ Grid<Real> &Aminusk = *matrixA[6];
+
+ if (vecRhs.size() != 2)
+ errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step");
+ Grid<Real> &srcV = *vecRhs[0];
+ Grid<Real> &srcW = *vecRhs[1];
+
+ dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
+ src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
+ src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
+ src(i, j, k - 1) * Aminusk(i, j, k);
+
+ dst(i, j, k) += srcV(i, j + 1, k) * (*matrixA[7])(i, j, k) +
+ srcV(i - 1, j + 1, k) * (*matrixA[8])(i, j, k) +
+ srcV(i, j, k) * (*matrixA[9])(i, j, k) +
+ srcV(i - 1, j, k) * (*matrixA[10])(i, j, k) +
+ srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) +
+ srcW(i - 1, j, k + 1) * (*matrixA[12])(i, j, k) +
+ srcW(i, j, k) * (*matrixA[13])(i, j, k) +
+ srcW(i - 1, j, k) * (*matrixA[14])(i, j, k);
+ }
+ inline const FlagGrid &getArg0()
+ {
+ return flags;
+ }
+ typedef FlagGrid type0;
+ inline Grid<Real> &getArg1()
+ {
+ return dst;
+ }
+ typedef Grid<Real> type1;
+ inline const Grid<Real> &getArg2()
+ {
+ return src;
+ }
+ typedef Grid<Real> type2;
+ inline const std::vector<Grid<Real> *> &getArg3()
+ {
+ return matrixA;
+ }
+ typedef std::vector<Grid<Real> *> type3;
+ inline const std::vector<Grid<Real> *> &getArg4()
+ {
+ return vecRhs;
+ }
+ typedef std::vector<Grid<Real> *> type4;
void runMessage()
{
- debMsg("Executing kernel ApplyMatrix2D ", 3);
+ debMsg("Executing kernel ApplyMatrixViscosityU ", 3);
debMsg("Kernel range"
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
- for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
- op(idx, flags, dst, src, A0, Ai, Aj, Ak);
+ const int _maxX = maxX;
+ const int _maxY = maxY;
+ if (maxZ > 1) {
+ for (int k = __r.begin(); k != (int)__r.end(); k++)
+ for (int j = 1; j < _maxY; j++)
+ for (int i = 1; i < _maxX; i++)
+ op(i, j, k, flags, dst, src, matrixA, vecRhs);
+ }
+ else {
+ const int k = 0;
+ for (int j = __r.begin(); j != (int)__r.end(); j++)
+ for (int i = 1; i < _maxX; i++)
+ op(i, j, k, flags, dst, src, matrixA, vecRhs);
+ }
}
void run()
{
- tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+ if (maxZ > 1)
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+ else
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
}
const FlagGrid &flags;
Grid<Real> &dst;
const Grid<Real> &src;
- Grid<Real> &A0;
- Grid<Real> &Ai;
- Grid<Real> &Aj;
- Grid<Real> &Ak;
+ const std::vector<Grid<Real> *> matrixA;
+ const std::vector<Grid<Real> *> vecRhs;
};
+struct ApplyMatrixViscosityV : public KernelBase {
+ ApplyMatrixViscosityV(const FlagGrid &flags,
+ Grid<Real> &dst,
+ const Grid<Real> &src,
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs)
+ : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(int i,
+ int j,
+ int k,
+ const FlagGrid &flags,
+ Grid<Real> &dst,
+ const Grid<Real> &src,
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs) const
+ {
+ if (matrixA.size() != 15)
+ errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+ Grid<Real> &A0 = *matrixA[0];
+ Grid<Real> &Aplusi = *matrixA[1];
+ Grid<Real> &Aplusj = *matrixA[2];
+ Grid<Real> &Aplusk = *matrixA[3];
+ Grid<Real> &Aminusi = *matrixA[4];
+ Grid<Real> &Aminusj = *matrixA[5];
+ Grid<Real> &Aminusk = *matrixA[6];
+
+ if (vecRhs.size() != 2)
+ errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step");
+ Grid<Real> &srcU = *vecRhs[0];
+ Grid<Real> &srcW = *vecRhs[1];
+
+ dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
+ src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
+ src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
+ src(i, j, k - 1) * Aminusk(i, j, k);
+
+ dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) +
+ srcU(i + 1, j - 1, k) * (*matrixA[8])(i, j, k) +
+ srcU(i, j, k) * (*matrixA[9])(i, j, k) +
+ srcU(i, j - 1, k) * (*matrixA[10])(i, j, k) +
+ srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) +
+ srcW(i, j - 1, k + 1) * (*matrixA[12])(i, j, k) +
+ srcW(i, j, k) * (*matrixA[13])(i, j, k) +
+ srcW(i, j - 1, k) * (*matrixA[14])(i, j, k);
+ }
+ inline const FlagGrid &getArg0()
+ {
+ return flags;
+ }
+ typedef FlagGrid type0;
+ inline Grid<Real> &getArg1()
+ {
+ return dst;
+ }
+ typedef Grid<Real> type1;
+ inline const Grid<Real> &getArg2()
+ {
+ return src;
+ }
+ typedef Grid<Real> type2;
+ inline const std::vector<Grid<Real> *> &getArg3()
+ {
+ return matrixA;
+ }
+ typedef std::vector<Grid<Real> *> type3;
+ inline const std::vector<Grid<Real> *> &getArg4()
+ {
+ return vecRhs;
+ }
+ typedef std::vector<Grid<Real> *> type4;
+ void runMessage()
+ {
+ debMsg("Executing kernel ApplyMatrixViscosityV ", 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 = 1; j < _maxY; j++)
+ for (int i = 1; i < _maxX; i++)
+ op(i, j, k, flags, dst, src, matrixA, vecRhs);
+ }
+ else {
+ const int k = 0;
+ for (int j = __r.begin(); j != (int)__r.end(); j++)
+ for (int i = 1; i < _maxX; i++)
+ op(i, j, k, flags, dst, src, matrixA, vecRhs);
+ }
+ }
+ void run()
+ {
+ if (maxZ > 1)
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+ else
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
+ }
+ const FlagGrid &flags;
+ Grid<Real> &dst;
+ const Grid<Real> &src;
+ const std::vector<Grid<Real> *> matrixA;
+ const std::vector<Grid<Real> *> vecRhs;
+};
+
+struct ApplyMatrixViscosityW : public KernelBase {
+ ApplyMatrixViscosityW(const FlagGrid &flags,
+ Grid<Real> &dst,
+ const Grid<Real> &src,
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs)
+ : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(int i,
+ int j,
+ int k,
+ const FlagGrid &flags,
+ Grid<Real> &dst,
+ const Grid<Real> &src,
+ const std::vector<Grid<Real> *> matrixA,
+ const std::vector<Grid<Real> *> vecRhs) const
+ {
+ if (matrixA.size() != 15)
+ errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+ Grid<Real> &A0 = *matrixA[0];
+ Grid<Real> &Aplusi = *matrixA[1];
+ Grid<Real> &Aplusj = *matrixA[2];
+ Grid<Real> &Aplusk = *matrixA[3];
+ Grid<Real> &Aminusi = *matrixA[4];
+ Grid<Real> &Aminusj = *matrixA[5];
+ Grid<Real> &Aminusk = *matrixA[6];
+
+ if (vecRhs.size() != 2)
+ errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step");
+ Grid<Real> &srcU = *vecRhs[0];
+ Grid<Real> &srcV = *vecRhs[1];
+
+ dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
+ src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
+ src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
+ src(i, j, k - 1) * Aminusk(i, j, k);
+
+ dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) +
+ srcU(i + 1, j, k - 1) * (*matrixA[8])(i, j, k) +
+ srcU(i, j, k) * (*matrixA[9])(i, j, k) +
+ srcU(i, j, k - 1) * (*matrixA[10])(i, j, k) +
+ srcV(i, j + 1, k) * (*matrixA[11])(i, j, k) +
+ srcV(i, j + 1, k - 1) * (*matrixA[12])(i, j, k) +
+ srcV(i, j, k) * (*matrixA[13])(i, j, k) +
+ srcV(i, j, k - 1) * (*matrixA[14])(i, j, k);
+ }
+ inline const FlagGrid &getArg0()
+ {
+ return flags;
+ }
+ typedef FlagGrid type0;
+ inline Grid<Real> &getArg1()
+ {
+ return dst;
+ }
+ typedef Grid<Real> type1;
+ inline const Grid<Real> &getArg2()
+ {
+ return src;
+ }
+ typedef Grid<Real> type2;
+ inline const std::vector<Grid<Real> *> &getArg3()
+ {
+ return matrixA;
+ }
+ typedef std::vector<Grid<Real> *> type3;
+ inline const std::vector<Grid<Real> *> &getArg4()
+ {
+ return vecRhs;
+ }
+ typedef std::vector<Grid<Real> *> type4;
+ void runMessage()
+ {
+ debMsg("Executing kernel ApplyMatrixViscosityW ", 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 = 1; j < _maxY; j++)
+ for (int i = 1; i < _maxX; i++)
+ op(i, j, k, flags, dst, src, matrixA, vecRhs);
+ }
+ else {
+ const int k = 0;
+ for (int j = __r.begin(); j != (int)__r.end(); j++)
+ for (int i = 1; i < _maxX; i++)
+ op(i, j, k, flags, dst, src, matrixA, vecRhs);
+ }
+ }
+ void run()
+ {
+ if (maxZ > 1)
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+ else
+ tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
+ }
+ const FlagGrid &flags;
+ Grid<Real> &dst;
+ const Grid<Real> &src;
+ const std::vector<Grid<Real> *> matrixA;
+ const std::vector<Grid<Real> *> vecRhs;
+};
+
+/* NOTE: Use this template for new matrix application kernels
+
+//! Template for matrix application kernels
+KERNEL()
+void ApplyMatrixTemplate (const FlagGrid& flags, Grid<Real>& dst, const Grid<Real>& src,
+ const std::vector<Grid<Real> *> matrixA, const std::vector<Grid<Real> *> vecRhs)
+{
+ // The kernel must define how to use the grids from the matrixA and vecRhs lists
+}
+
+*/
+
//! Kernel: Construct the matrix for the poisson equation
struct MakeLaplaceMatrix : public KernelBase {
diff --git a/extern/mantaflow/preprocessed/general.h b/extern/mantaflow/preprocessed/general.h
index 7a840517cef..50eac71e87e 100644
--- a/extern/mantaflow/preprocessed/general.h
+++ b/extern/mantaflow/preprocessed/general.h
@@ -42,7 +42,7 @@ inline void updateQtGui(bool full, int frame, float time, const std::string &cur
# ifdef _DEBUG
# define DEBUG 1
# endif // _DEBUG
-#endif // DEBUG
+#endif // DEBUG
// Standard exception
class Error : public std::exception {
@@ -242,6 +242,39 @@ inline bool c_isnan(float c)
return d != d;
}
+//! Swap so that a<b
+template<class T> inline void sort(T &a, T &b)
+{
+ if (a > b)
+ std::swap(a, b);
+}
+
+//! Swap so that a<b<c
+template<class T> inline void sort(T &a, T &b, T &c)
+{
+ if (a > b)
+ std::swap(a, b);
+ if (a > c)
+ std::swap(a, c);
+ if (b > c)
+ std::swap(b, c);
+}
+
+//! Swap so that a<b<c<d
+template<class T> inline void sort(T &a, T &b, T &c, T &d)
+{
+ if (a > b)
+ std::swap(a, b);
+ if (c > d)
+ std::swap(c, d);
+ if (a > c)
+ std::swap(a, c);
+ if (b > d)
+ std::swap(b, d);
+ if (b > c)
+ std::swap(b, c);
+}
+
} // namespace Manta
#endif
diff --git a/extern/mantaflow/preprocessed/gitinfo.h b/extern/mantaflow/preprocessed/gitinfo.h
index a0590e6a0b0..afc34a00d3d 100644
--- a/extern/mantaflow/preprocessed/gitinfo.h
+++ b/extern/mantaflow/preprocessed/gitinfo.h
@@ -1,3 +1,3 @@
-#define MANTA_GIT_VERSION "commit 327917cd59b03bef3a953b5f58fc1637b3a83e01"
+#define MANTA_GIT_VERSION "commit e2285cb9bc492987f728123be6cfc1fe11fe73d6"
diff --git a/extern/mantaflow/preprocessed/plugin/extforces.cpp b/extern/mantaflow/preprocessed/plugin/extforces.cpp
index 558008afe3b..88935fa7ae9 100644
--- a/extern/mantaflow/preprocessed/plugin/extforces.cpp
+++ b/extern/mantaflow/preprocessed/plugin/extforces.cpp
@@ -1135,26 +1135,27 @@ struct KnAddForceIfLower : public KernelBase {
if (!curFluid && !curEmpty)
return;
+ Real minVal, maxVal, sum;
if (flags.isFluid(i - 1, j, k) || (curFluid && flags.isEmpty(i - 1, j, k))) {
Real forceMACX = 0.5 * (force(i - 1, j, k).x + force(i, j, k).x);
- Real min = std::min(vel(i, j, k).x, forceMACX);
- Real max = std::max(vel(i, j, k).x, forceMACX);
- Real sum = vel(i, j, k).x + forceMACX;
- vel(i, j, k).x = (forceMACX > 0) ? std::min(sum, max) : std::max(sum, min);
+ minVal = min(vel(i, j, k).x, forceMACX);
+ maxVal = max(vel(i, j, k).x, forceMACX);
+ sum = vel(i, j, k).x + forceMACX;
+ vel(i, j, k).x = (forceMACX > 0) ? min(sum, maxVal) : max(sum, minVal);
}
if (flags.isFluid(i, j - 1, k) || (curFluid && flags.isEmpty(i, j - 1, k))) {
Real forceMACY = 0.5 * (force(i, j - 1, k).y + force(i, j, k).y);
- Real min = std::min(vel(i, j, k).y, forceMACY);
- Real max = std::max(vel(i, j, k).y, forceMACY);
- Real sum = vel(i, j, k).y + forceMACY;
- vel(i, j, k).y = (forceMACY > 0) ? std::min(sum, max) : std::max(sum, min);
+ minVal = min(vel(i, j, k).y, forceMACY);
+ maxVal = max(vel(i, j, k).y, forceMACY);
+ sum = vel(i, j, k).y + forceMACY;
+ vel(i, j, k).y = (forceMACY > 0) ? min(sum, maxVal) : max(sum, minVal);
}
if (vel.is3D() && (flags.isFluid(i, j, k - 1) || (curFluid && flags.isEmpty(i, j, k - 1)))) {
Real forceMACZ = 0.5 * (force(i, j, k - 1).z + force(i, j, k).z);
- Real min = std::min(vel(i, j, k).z, forceMACZ);
- Real max = std::max(vel(i, j, k).z, forceMACZ);
- Real sum = vel(i, j, k).z + forceMACZ;
- vel(i, j, k).z = (forceMACZ > 0) ? std::min(sum, max) : std::max(sum, min);
+ minVal = min(vel(i, j, k).z, forceMACZ);
+ maxVal = max(vel(i, j, k).z, forceMACZ);
+ sum = vel(i, j, k).z + forceMACZ;
+ vel(i, j, k).z = (forceMACZ > 0) ? min(sum, maxVal) : max(sum, minVal);
}
}
inline const FlagGrid &getArg0()
diff --git a/extern/mantaflow/preprocessed/plugin/pressure.cpp b/extern/mantaflow/preprocessed/plugin/pressure.cpp
index 1100a58db47..593aeb16859 100644
--- a/extern/mantaflow/preprocessed/plugin/pressure.cpp
+++ b/extern/mantaflow/preprocessed/plugin/pressure.cpp
@@ -1138,11 +1138,15 @@ void solvePressureSystem(Grid<Real> &rhs,
// note: the last factor increases the max iterations for 2d, which right now can't use a
// preconditioner
GridCgInterface *gcg;
- if (vel.is3D())
- gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
- else
- gcg = new GridCg<ApplyMatrix2D>(
- pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+ vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
+
+ if (vel.is3D()) {
+ matA.push_back(&Ak);
+ gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, matA);
+ }
+ else {
+ gcg = new GridCg<ApplyMatrix2D>(pressure, rhs, residual, search, flags, tmp, matA);
+ }
gcg->setAccuracy(cgAccuracy);
gcg->setUseL2Norm(useL2Norm);
diff --git a/extern/mantaflow/preprocessed/plugin/viscosity.cpp b/extern/mantaflow/preprocessed/plugin/viscosity.cpp
new file mode 100644
index 00000000000..5dc478b13dd
--- /dev/null
+++ b/extern/mantaflow/preprocessed/plugin/viscosity.cpp
@@ -0,0 +1,1430 @@
+
+
+// DO NOT EDIT !
+// This file is generated using the MantaFlow preprocessor (prep generate).
+
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2020 Sebastian Barschkis, 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
+ *
+ * Accurate Viscous Free Surfaces for Buckling, Coiling, and Rotating Liquids
+ * Batty et al., SCA 2008
+ *
+ ******************************************************************************/
+
+#include "conjugategrad.h"
+#include "general.h"
+#include "grid.h"
+#include "vectorbase.h"
+
+#include <chrono>
+
+#if OPENMP == 1 || TBB == 1
+# define ENABLE_PARALLEL 0
+#endif
+
+#if ENABLE_PARALLEL == 1
+# include <thread>
+# include <algorithm>
+
+static const int manta_num_threads = std::thread::hardware_concurrency();
+
+# define parallel_block \
+ { \
+ std::vector<std::thread> threads; \
+ {
+
+# define do_parallel threads.push_back( std::thread([&]() {
+# define do_end \
+ } ) );
+
+# define block_end \
+ } \
+ for (auto &thread : threads) { \
+ thread.join(); \
+ } \
+ }
+
+#endif
+
+#define FOR_INT_IJK(num) \
+ for (int k_off = 0; k_off < num; ++k_off) \
+ for (int j_off = 0; j_off < num; ++j_off) \
+ for (int i_off = 0; i_off < num; ++i_off)
+
+using namespace std;
+
+namespace Manta {
+
+//! Assumes phi0<0 and phi1>=0, phi2>=0, and phi3>=0 or vice versa.
+//! In particular, phi0 must not equal any of phi1, phi2 or phi3.
+static Real sortedTetFraction(Real phi0, Real phi1, Real phi2, Real phi3)
+{
+ return phi0 * phi0 * phi0 / ((phi0 - phi1) * (phi0 - phi2) * (phi0 - phi3));
+}
+
+//! Assumes phi0<0, phi1<0, and phi2>=0, and phi3>=0 or vice versa.
+//! In particular, phi0 and phi1 must not equal any of phi2 and phi3.
+static Real sortedPrismFraction(Real phi0, Real phi1, Real phi2, Real phi3)
+{
+ Real a = phi0 / (phi0 - phi2);
+ Real b = phi0 / (phi0 - phi3);
+ Real c = phi1 / (phi1 - phi3);
+ Real d = phi1 / (phi1 - phi2);
+ return a * b * (1 - d) + b * (1 - c) * d + c * d;
+}
+
+Real volumeFraction(Real phi0, Real phi1, Real phi2, Real phi3)
+{
+ sort(phi0, phi1, phi2, phi3);
+ if (phi3 <= 0)
+ return 1;
+ else if (phi2 <= 0)
+ return 1 - sortedTetFraction(phi3, phi2, phi1, phi0);
+ else if (phi1 <= 0)
+ return sortedPrismFraction(phi0, phi1, phi2, phi3);
+ else if (phi0 <= 0)
+ return sortedTetFraction(phi0, phi1, phi2, phi3);
+ else
+ return 0;
+}
+
+//! The average of the two possible decompositions of the cube into five tetrahedra.
+Real volumeFraction(Real phi000,
+ Real phi100,
+ Real phi010,
+ Real phi110,
+ Real phi001,
+ Real phi101,
+ Real phi011,
+ Real phi111)
+{
+ return (volumeFraction(phi000, phi001, phi101, phi011) +
+ volumeFraction(phi000, phi101, phi100, phi110) +
+ volumeFraction(phi000, phi010, phi011, phi110) +
+ volumeFraction(phi101, phi011, phi111, phi110) +
+ 2 * volumeFraction(phi000, phi011, phi101, phi110) +
+ volumeFraction(phi100, phi101, phi001, phi111) +
+ volumeFraction(phi100, phi001, phi000, phi010) +
+ volumeFraction(phi100, phi110, phi111, phi010) +
+ volumeFraction(phi001, phi111, phi011, phi010) +
+ 2 * volumeFraction(phi100, phi111, phi001, phi010)) /
+ 12;
+}
+
+//! Kernel loop over grid with 2x base resolution!
+
+struct KnEstimateVolumeFraction : public KernelBase {
+ KnEstimateVolumeFraction(Grid<Real> &volumes,
+ const Grid<Real> &phi,
+ const Vec3 &startCentre,
+ const Real dx)
+ : KernelBase(&volumes, 0), volumes(volumes), phi(phi), startCentre(startCentre), dx(dx)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(int i,
+ int j,
+ int k,
+ Grid<Real> &volumes,
+ const Grid<Real> &phi,
+ const Vec3 &startCentre,
+ const Real dx) const
+ {
+ const Vec3 centre = startCentre + Vec3(i, j, k) * 0.5;
+ const Real offset = 0.5 * dx;
+ const int order = 2;
+
+ Real phi000 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, -offset), order);
+ Real phi001 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, +offset), order);
+ Real phi010 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, -offset), order);
+ Real phi011 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, +offset), order);
+ Real phi100 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, -offset), order);
+ Real phi101 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, +offset), order);
+ Real phi110 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, -offset), order);
+ Real phi111 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, +offset), order);
+
+ volumes(i, j, k) = volumeFraction(
+ phi000, phi100, phi010, phi110, phi001, phi101, phi011, phi111);
+ }
+ inline Grid<Real> &getArg0()
+ {
+ return volumes;
+ }
+ typedef Grid<Real> type0;
+ inline const Grid<Real> &getArg1()
+ {
+ return phi;
+ }
+ typedef Grid<Real> type1;
+ inline const Vec3 &getArg2()
+ {
+ return startCentre;
+ }
+ typedef Vec3 type2;
+ inline const Real &getArg3()
+ {
+ return dx;
+ }
+ typedef Real type3;
+ void runMessage()
+ {
+ debMsg("Executing kernel KnEstimateVolumeFraction ", 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, volumes, phi, startCentre, dx);
+ }
+ 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, volumes, phi, startCentre, dx);
+ }
+ }
+ 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<Real> &volumes;
+ const Grid<Real> &phi;
+ const Vec3 &startCentre;
+ const Real dx;
+};
+
+struct KnUpdateVolumeGrid : public KernelBase {
+ KnUpdateVolumeGrid(Grid<Real> &cVolLiquid,
+ Grid<Real> &uVolLiquid,
+ Grid<Real> &vVolLiquid,
+ Grid<Real> &wVolLiquid,
+ Grid<Real> &exVolLiquid,
+ Grid<Real> &eyVolLiquid,
+ Grid<Real> &ezVolLiquid,
+ const Grid<Real> &src)
+ : KernelBase(&cVolLiquid, 0),
+ cVolLiquid(cVolLiquid),
+ uVolLiquid(uVolLiquid),
+ vVolLiquid(vVolLiquid),
+ wVolLiquid(wVolLiquid),
+ exVolLiquid(exVolLiquid),
+ eyVolLiquid(eyVolLiquid),
+ ezVolLiquid(ezVolLiquid),
+ src(src)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(int i,
+ int j,
+ int k,
+ Grid<Real> &cVolLiquid,
+ Grid<Real> &uVolLiquid,
+ Grid<Real> &vVolLiquid,
+ Grid<Real> &wVolLiquid,
+ Grid<Real> &exVolLiquid,
+ Grid<Real> &eyVolLiquid,
+ Grid<Real> &ezVolLiquid,
+ const Grid<Real> &src) const
+ {
+ // Work out c
+ cVolLiquid(i, j, k) = 0;
+ FOR_INT_IJK(2)
+ {
+ cVolLiquid(i, j, k) += src(2 * i + i_off, 2 * j + j_off, 2 * k + k_off);
+ }
+ cVolLiquid(i, j, k) /= 8;
+
+ // Work out u
+ if (i >= 1) {
+ uVolLiquid(i, j, k) = 0;
+ int base_i = 2 * i - 1;
+ int base_j = 2 * j;
+ int base_k = 2 * k;
+ FOR_INT_IJK(2)
+ {
+ uVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
+ }
+ uVolLiquid(i, j, k) /= 8;
+ }
+
+ // v
+ if (j >= 1) {
+ vVolLiquid(i, j, k) = 0;
+ int base_i = 2 * i;
+ int base_j = 2 * j - 1;
+ int base_k = 2 * k;
+ FOR_INT_IJK(2)
+ {
+ vVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
+ }
+ vVolLiquid(i, j, k) /= 8;
+ }
+
+ // w
+ if (k >= 1) {
+ wVolLiquid(i, j, k) = 0;
+ int base_i = 2 * i;
+ int base_j = 2 * j;
+ int base_k = 2 * k - 1;
+ FOR_INT_IJK(2)
+ {
+ wVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
+ }
+ wVolLiquid(i, j, k) /= 8;
+ }
+
+ // e-x
+ if (j >= 1 && k >= 1) {
+ exVolLiquid(i, j, k) = 0;
+ int base_i = 2 * i;
+ int base_j = 2 * j - 1;
+ int base_k = 2 * k - 1;
+ FOR_INT_IJK(2)
+ {
+ exVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
+ }
+ exVolLiquid(i, j, k) /= 8;
+ }
+
+ // e-y
+ if (i >= 1 && k >= 1) {
+ eyVolLiquid(i, j, k) = 0;
+ int base_i = 2 * i - 1;
+ int base_j = 2 * j;
+ int base_k = 2 * k - 1;
+ FOR_INT_IJK(2)
+ {
+ eyVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
+ }
+ eyVolLiquid(i, j, k) /= 8;
+ }
+
+ // e-z
+ if (i >= 1 && j >= 1) {
+ ezVolLiquid(i, j, k) = 0;
+ int base_i = 2 * i - 1;
+ int base_j = 2 * j - 1;
+ int base_k = 2 * k;
+ FOR_INT_IJK(2)
+ {
+ ezVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off);
+ }
+ ezVolLiquid(i, j, k) /= 8;
+ }
+ }
+ inline Grid<Real> &getArg0()
+ {
+ return cVolLiquid;
+ }
+ typedef Grid<Real> type0;
+ inline Grid<Real> &getArg1()
+ {
+ return uVolLiquid;
+ }
+ typedef Grid<Real> type1;
+ inline Grid<Real> &getArg2()
+ {
+ return vVolLiquid;
+ }
+ typedef Grid<Real> type2;
+ inline Grid<Real> &getArg3()
+ {
+ return wVolLiquid;
+ }
+ typedef Grid<Real> type3;
+ inline Grid<Real> &getArg4()
+ {
+ return exVolLiquid;
+ }
+ typedef Grid<Real> type4;
+ inline Grid<Real> &getArg5()
+ {
+ return eyVolLiquid;
+ }
+ typedef Grid<Real> type5;
+ inline Grid<Real> &getArg6()
+ {
+ return ezVolLiquid;
+ }
+ typedef Grid<Real> type6;
+ inline const Grid<Real> &getArg7()
+ {
+ return src;
+ }
+ typedef Grid<Real> type7;
+ void runMessage()
+ {
+ debMsg("Executing kernel KnUpdateVolumeGrid ", 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,
+ cVolLiquid,
+ uVolLiquid,
+ vVolLiquid,
+ wVolLiquid,
+ exVolLiquid,
+ eyVolLiquid,
+ ezVolLiquid,
+ src);
+ }
+ 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,
+ cVolLiquid,
+ uVolLiquid,
+ vVolLiquid,
+ wVolLiquid,
+ exVolLiquid,
+ eyVolLiquid,
+ ezVolLiquid,
+ src);
+ }
+ }
+ 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<Real> &cVolLiquid;
+ Grid<Real> &uVolLiquid;
+ Grid<Real> &vVolLiquid;
+ Grid<Real> &wVolLiquid;
+ Grid<Real> &exVolLiquid;
+ Grid<Real> &eyVolLiquid;
+ Grid<Real> &ezVolLiquid;
+ const Grid<Real> &src;
+};
+
+void computeWeights(const Grid<Real> &phi,
+ Grid<Real> &doubleSized,
+ Grid<Real> &cVolLiquid,
+ Grid<Real> &uVolLiquid,
+ Grid<Real> &vVolLiquid,
+ Grid<Real> &wVolLiquid,
+ Grid<Real> &exVolLiquid,
+ Grid<Real> &eyVolLiquid,
+ Grid<Real> &ezVolLiquid,
+ Real dx)
+{
+ KnEstimateVolumeFraction(doubleSized, phi, Vec3(0.25 * dx, 0.25 * dx, 0.25 * dx), 0.5 * dx);
+ KnUpdateVolumeGrid(cVolLiquid,
+ uVolLiquid,
+ vVolLiquid,
+ wVolLiquid,
+ exVolLiquid,
+ eyVolLiquid,
+ ezVolLiquid,
+ doubleSized);
+}
+
+struct KnUpdateFaceStates : public KernelBase {
+ KnUpdateFaceStates(const FlagGrid &flags,
+ Grid<int> &uState,
+ Grid<int> &vState,
+ Grid<int> &wState)
+ : KernelBase(&flags, 0), flags(flags), uState(uState), vState(vState), wState(wState)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(int i,
+ int j,
+ int k,
+ const FlagGrid &flags,
+ Grid<int> &uState,
+ Grid<int> &vState,
+ Grid<int> &wState) const
+ {
+ bool curObs = flags.isObstacle(i, j, k);
+ uState(i, j, k) = (i > 0 && !flags.isObstacle(i - 1, j, k) && !curObs) ?
+ FlagGrid::TypeFluid :
+ FlagGrid::TypeObstacle;
+ vState(i, j, k) = (j > 0 && !flags.isObstacle(i, j - 1, k) && !curObs) ?
+ FlagGrid::TypeFluid :
+ FlagGrid::TypeObstacle;
+ wState(i, j, k) = (k > 0 && !flags.isObstacle(i, j, k - 1) && !curObs) ?
+ FlagGrid::TypeFluid :
+ FlagGrid::TypeObstacle;
+ }
+ inline const FlagGrid &getArg0()
+ {
+ return flags;
+ }
+ typedef FlagGrid type0;
+ inline Grid<int> &getArg1()
+ {
+ return uState;
+ }
+ typedef Grid<int> type1;
+ inline Grid<int> &getArg2()
+ {
+ return vState;
+ }
+ typedef Grid<int> type2;
+ inline Grid<int> &getArg3()
+ {
+ return wState;
+ }
+ typedef Grid<int> type3;
+ void runMessage()
+ {
+ debMsg("Executing kernel KnUpdateFaceStates ", 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, flags, uState, vState, wState);
+ }
+ 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, flags, uState, vState, wState);
+ }
+ }
+ 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);
+ }
+ const FlagGrid &flags;
+ Grid<int> &uState;
+ Grid<int> &vState;
+ Grid<int> &wState;
+};
+
+struct KnApplyVelocities : public KernelBase {
+ KnApplyVelocities(MACGrid &dst,
+ const Grid<int> &uState,
+ const Grid<int> &vState,
+ const Grid<int> &wState,
+ Grid<Real> &srcU,
+ Grid<Real> &srcV,
+ Grid<Real> &srcW)
+ : KernelBase(&dst, 0),
+ dst(dst),
+ uState(uState),
+ vState(vState),
+ wState(wState),
+ srcU(srcU),
+ srcV(srcV),
+ srcW(srcW)
+ {
+ runMessage();
+ run();
+ }
+ inline void op(int i,
+ int j,
+ int k,
+ MACGrid &dst,
+ const Grid<int> &uState,
+ const Grid<int> &vState,
+ const Grid<int> &wState,
+ Grid<Real> &srcU,
+ Grid<Real> &srcV,
+ Grid<Real> &srcW) const
+ {
+ dst(i, j, k).x = (uState(i, j, k) == FlagGrid::TypeFluid) ? srcU(i, j, k) : 0;
+ dst(i, j, k).y = (vState(i, j, k) == FlagGrid::TypeFluid) ? srcV(i, j, k) : 0;
+ if (dst.is3D())
+ dst(i, j, k).z = (wState(i, j, k) == FlagGrid::TypeFluid) ? srcW(i, j, k) : 0;
+ }
+ inline MACGrid &getArg0()
+ {
+ return dst;
+ }
+ typedef MACGrid type0;
+ inline const Grid<int> &getArg1()
+ {
+ return uState;
+ }
+ typedef Grid<int> type1;
+ inline const Grid<int> &getArg2()
+ {
+ return vState;
+ }
+ typedef Grid<int> type2;
+ inline const Grid<int> &getArg3()
+ {
+ return wState;
+ }
+ typedef Grid<int> type3;
+ inline Grid<Real> &getArg4()
+ {
+ return srcU;
+ }
+ typedef Grid<Real> type4;
+ inline Grid<Real> &getArg5()
+ {
+ return srcV;
+ }
+ typedef Grid<Real> type5;
+ inline Grid<Real> &getArg6()
+ {
+ return srcW;
+ }
+ typedef Grid<Real> type6;
+ void runMessage()
+ {
+ debMsg("Executing kernel KnApplyVelocities ", 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, dst, uState, vState, wState, srcU, srcV, srcW);
+ }
+ 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, dst, uState, vState, wState, srcU, srcV, srcW);
+ }
+ }
+ 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);
+ }
+ MACGrid &dst;
+ const Grid<int> &uState;
+ const Grid<int> &vState;
+ const Grid<int> &wState;
+ Grid<Real> &srcU;
+ Grid<Real> &srcV;
+ Grid<Real> &srcW;
+};
+
+void solveViscosity(const FlagGrid &flags,
+ MACGrid &vel,
+ Grid<Real> &cVolLiquid,
+ Grid<Real> &uVolLiquid,
+ Grid<Real> &vVolLiquid,
+ Grid<Real> &wVolLiquid,
+ Grid<Real> &exVolLiquid,
+ Grid<Real> &eyVolLiquid,
+ Grid<Real> &ezVolLiquid,
+ Grid<Real> &viscosity,
+ const Real dt,
+ const Real dx,
+ const Real cgAccuracy,
+ const Real cgMaxIterFac)
+{
+ const Real factor = dt * square(1.0 / dx);
+ const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4);
+ GridCg<ApplyMatrixViscosityU> *uGcg;
+ GridCg<ApplyMatrixViscosityV> *vGcg;
+ GridCg<ApplyMatrixViscosityW> *wGcg;
+
+ // Tmp grids for CG solve in U, V, W dimensions
+ FluidSolver *parent = flags.getParent();
+ Grid<Real> uResidual(parent);
+ Grid<Real> vResidual(parent);
+ Grid<Real> wResidual(parent);
+ Grid<Real> uSearch(parent);
+ Grid<Real> vSearch(parent);
+ Grid<Real> wSearch(parent);
+ Grid<Real> uTmp(parent);
+ Grid<Real> vTmp(parent);
+ Grid<Real> wTmp(parent);
+ Grid<Real> uRhs(parent);
+ Grid<Real> vRhs(parent);
+ Grid<Real> wRhs(parent);
+
+ // A matrix U grids
+ Grid<Real> uA0(parent); // diagonal elements in A matrix
+ Grid<Real> uAplusi(parent); // neighbor at i+1
+ Grid<Real> uAplusj(parent); // neighbor at j+1
+ Grid<Real> uAplusk(parent); // neighbor at k+1
+ Grid<Real> uAminusi(parent); // neighbor at i-1
+ Grid<Real> uAminusj(parent); // neighbor at j-1
+ Grid<Real> uAminusk(parent); // neighbor at k-1
+ Grid<Real> uAhelper1(parent); // additional helper grids for off diagonal elements
+ Grid<Real> uAhelper2(parent);
+ Grid<Real> uAhelper3(parent);
+ Grid<Real> uAhelper4(parent);
+ Grid<Real> uAhelper5(parent);
+ Grid<Real> uAhelper6(parent);
+ Grid<Real> uAhelper7(parent);
+ Grid<Real> uAhelper8(parent);
+
+ // A matrix V grids
+ Grid<Real> vA0(parent);
+ Grid<Real> vAplusi(parent);
+ Grid<Real> vAplusj(parent);
+ Grid<Real> vAplusk(parent);
+ Grid<Real> vAminusi(parent);
+ Grid<Real> vAminusj(parent);
+ Grid<Real> vAminusk(parent);
+ Grid<Real> vAhelper1(parent);
+ Grid<Real> vAhelper2(parent);
+ Grid<Real> vAhelper3(parent);
+ Grid<Real> vAhelper4(parent);
+ Grid<Real> vAhelper5(parent);
+ Grid<Real> vAhelper6(parent);
+ Grid<Real> vAhelper7(parent);
+ Grid<Real> vAhelper8(parent);
+
+ // A matrix W grids
+ Grid<Real> wA0(parent);
+ Grid<Real> wAplusi(parent);
+ Grid<Real> wAplusj(parent);
+ Grid<Real> wAplusk(parent);
+ Grid<Real> wAminusi(parent);
+ Grid<Real> wAminusj(parent);
+ Grid<Real> wAminusk(parent);
+ Grid<Real> wAhelper1(parent);
+ Grid<Real> wAhelper2(parent);
+ Grid<Real> wAhelper3(parent);
+ Grid<Real> wAhelper4(parent);
+ Grid<Real> wAhelper5(parent);
+ Grid<Real> wAhelper6(parent);
+ Grid<Real> wAhelper7(parent);
+ Grid<Real> wAhelper8(parent);
+
+ // Solution grids for CG solvers
+ Grid<Real> uSolution(parent);
+ Grid<Real> vSolution(parent);
+ Grid<Real> wSolution(parent);
+
+ // Save state of voxel face (fluid or obstacle)
+ Grid<int> uState(parent);
+ Grid<int> vState(parent);
+ Grid<int> wState(parent);
+
+ // Save state of voxel face (fluid or obstacle)
+ KnUpdateFaceStates(flags, uState, vState, wState);
+
+ // Shorter names for flags, we will use them often
+ int isFluid = FlagGrid::TypeFluid;
+ int isObstacle = FlagGrid::TypeObstacle;
+
+ // Main viscosity loop: construct A matrices and rhs's in all dimensions
+ FOR_IJK_BND(flags, 1)
+ {
+
+ // U-terms: 2u_xx+ v_xy +uyy + u_zz + w_xz
+ if (uState(i, j, k) == isFluid) {
+
+ uRhs(i, j, k) = uVolLiquid(i, j, k) * vel(i, j, k).x;
+ uA0(i, j, k) = uVolLiquid(i, j, k);
+
+ Real viscRight = viscosity(i, j, k);
+ Real viscLeft = viscosity(i - 1, j, k);
+ Real volRight = cVolLiquid(i, j, k);
+ Real volLeft = cVolLiquid(i - 1, j, k);
+
+ Real viscTop = 0.25 * (viscosity(i - 1, j + 1, k) + viscosity(i - 1, j, k) +
+ viscosity(i, j + 1, k) + viscosity(i, j, k));
+ Real viscBottom = 0.25 * (viscosity(i - 1, j, k) + viscosity(i - 1, j - 1, k) +
+ viscosity(i, j, k) + viscosity(i, j - 1, k));
+ Real volTop = ezVolLiquid(i, j + 1, k);
+ Real volBottom = ezVolLiquid(i, j, k);
+
+ Real viscFront = 0.25 * (viscosity(i - 1, j, k + 1) + viscosity(i - 1, j, k) +
+ viscosity(i, j, k + 1) + viscosity(i, j, k));
+ Real viscBack = 0.25 * (viscosity(i - 1, j, k) + viscosity(i - 1, j, k - 1) +
+ viscosity(i, j, k) + viscosity(i, j, k - 1));
+ Real volFront = eyVolLiquid(i, j, k + 1);
+ Real volBack = eyVolLiquid(i, j, k);
+
+ Real factorRight = 2 * factor * viscRight * volRight;
+ Real factorLeft = 2 * factor * viscLeft * volLeft;
+ Real factorTop = factor * viscTop * volTop;
+ Real factorBottom = factor * viscBottom * volBottom;
+ Real factorFront = factor * viscFront * volFront;
+ Real factorBack = factor * viscBack * volBack;
+
+ // u_x_right
+ uA0(i, j, k) += factorRight;
+ if (uState(i + 1, j, k) == isFluid) {
+ uAplusi(i, j, k) += -factorRight;
+ }
+ else if (uState(i + 1, j, k) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight;
+ }
+
+ // u_x_left
+ uA0(i, j, k) += factorLeft;
+ if (uState(i - 1, j, k) == isFluid) {
+ uAminusi(i, j, k) += -factorLeft;
+ }
+ else if (uState(i - 1, j, k) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i - 1, j, k).x * factorLeft;
+ }
+
+ // u_y_top
+ uA0(i, j, k) += factorTop;
+ if (uState(i, j + 1, k) == isFluid) {
+ uAplusj(i, j, k) += -factorTop;
+ }
+ else if (uState(i, j + 1, k) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i, j + 1, k).x * factorTop;
+ }
+
+ // u_y_bottom
+ uA0(i, j, k) += factorBottom;
+ if (uState(i, j - 1, k) == isFluid) {
+ uAminusj(i, j, k) += -factorBottom;
+ }
+ else if (uState(i, j - 1, k) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i, j - 1, k).x * factorBottom;
+ }
+
+ // u_z_front
+ uA0(i, j, k) += factorFront;
+ if (uState(i, j, k + 1) == isFluid) {
+ uAplusk(i, j, k) += -factorFront;
+ }
+ else if (uState(i, j, k + 1) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i, j, k + 1).x * factorFront;
+ }
+
+ // u_z_back
+ uA0(i, j, k) += factorBack;
+ if (uState(i, j, k - 1) == isFluid) {
+ uAminusk(i, j, k) += -factorBack;
+ }
+ else if (uState(i, j, k - 1) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i, j, k - 1).x * factorBack;
+ }
+
+ // v_x_top
+ if (vState(i, j + 1, k) == isFluid) {
+ uAhelper1(i, j, k) += -factorTop;
+ }
+ else if (vState(i, j + 1, k) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop;
+ }
+
+ if (vState(i - 1, j + 1, k) == isFluid) {
+ uAhelper2(i, j, k) += factorTop;
+ }
+ else if (vState(i - 1, j + 1, k) == isObstacle) {
+ uRhs(i, j, k) -= vel(i - 1, j + 1, k).y * factorTop;
+ }
+
+ // v_x_bottom
+ if (vState(i, j, k) == isFluid) {
+ uAhelper3(i, j, k) += factorBottom;
+ }
+ else if (vState(i, j, k) == isObstacle) {
+ uRhs(i, j, k) -= vel(i, j, k).y * factorBottom;
+ }
+
+ if (vState(i - 1, j, k) == isFluid) {
+ uAhelper4(i, j, k) += -factorBottom;
+ }
+ else if (vState(i - 1, j, k) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i - 1, j, k).y * factorBottom;
+ }
+
+ // w_x_front
+ if (wState(i, j, k + 1) == isFluid) {
+ uAhelper5(i, j, k) += -factorFront;
+ }
+ else if (wState(i, j, k + 1) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront;
+ }
+
+ if (wState(i - 1, j, k + 1) == isFluid) {
+ uAhelper6(i, j, k) += factorFront;
+ }
+ else if (wState(i - 1, j, k + 1) == isObstacle) {
+ uRhs(i, j, k) -= vel(i - 1, j, k + 1).z * factorFront;
+ }
+
+ // w_x_back
+ if (wState(i, j, k) == isFluid) {
+ uAhelper7(i, j, k) += factorBack;
+ }
+ else if (wState(i, j, k) == isObstacle) {
+ uRhs(i, j, k) -= vel(i, j, k).z * factorBack;
+ }
+
+ if (wState(i - 1, j, k) == isFluid) {
+ uAhelper8(i, j, k) += -factorBack;
+ }
+ else if (wState(i - 1, j, k) == isObstacle) {
+ uRhs(i, j, k) -= -vel(i - 1, j, k).z * factorBack;
+ }
+ }
+
+ // V-terms: vxx + 2vyy + vzz + u_yx + w_yz
+ if (vState(i, j, k) == isFluid) {
+
+ vRhs(i, j, k) = vVolLiquid(i, j, k) * vel(i, j, k).y;
+ vA0(i, j, k) = vVolLiquid(i, j, k);
+
+ Real viscRight = 0.25 * (viscosity(i, j - 1, k) + viscosity(i + 1, j - 1, k) +
+ viscosity(i, j, k) + viscosity(i + 1, j, k));
+ Real viscLeft = 0.25 * (viscosity(i, j - 1, k) + viscosity(i - 1, j - 1, k) +
+ viscosity(i, j, k) + viscosity(i - 1, j, k));
+ Real volRight = ezVolLiquid(i + 1, j, k);
+ Real volLeft = ezVolLiquid(i, j, k);
+
+ Real viscTop = viscosity(i, j, k);
+ Real viscBottom = viscosity(i, j - 1, k);
+ Real volTop = cVolLiquid(i, j, k);
+ Real volBottom = cVolLiquid(i, j - 1, k);
+
+ Real viscFront = 0.25 * (viscosity(i, j - 1, k) + viscosity(i, j - 1, k + 1) +
+ viscosity(i, j, k) + viscosity(i, j, k + 1));
+ Real viscBack = 0.25 * (viscosity(i, j - 1, k) + viscosity(i, j - 1, k - 1) +
+ viscosity(i, j, k) + viscosity(i, j, k - 1));
+ Real volFront = exVolLiquid(i, j, k + 1);
+ Real volBack = exVolLiquid(i, j, k);
+
+ Real factorRight = factor * viscRight * volRight;
+ Real factorLeft = factor * viscLeft * volLeft;
+ Real factorTop = 2 * factor * viscTop * volTop;
+ Real factorBottom = 2 * factor * viscBottom * volBottom;
+ Real factorFront = factor * viscFront * volFront;
+ Real factorBack = factor * viscBack * volBack;
+
+ // v_x_right
+ vA0(i, j, k) += factorRight;
+ if (vState(i + 1, j, k) == isFluid) {
+ vAplusi(i, j, k) += -factorRight;
+ }
+ else if (vState(i + 1, j, k) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i + 1, j, k).y * factorRight;
+ }
+
+ // v_x_left
+ vA0(i, j, k) += factorLeft;
+ if (vState(i - 1, j, k) == isFluid) {
+ vAminusi(i, j, k) += -factorLeft;
+ }
+ else if (vState(i - 1, j, k) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i - 1, j, k).y * factorLeft;
+ }
+
+ // vy_top
+ vA0(i, j, k) += factorTop;
+ if (vState(i, j + 1, k) == isFluid) {
+ vAplusj(i, j, k) += -factorTop;
+ }
+ else if (vState(i, j + 1, k) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop;
+ }
+
+ // vy_bottom
+ vA0(i, j, k) += factorBottom;
+ if (vState(i, j - 1, k) == isFluid) {
+ vAminusj(i, j, k) += -factorBottom;
+ }
+ else if (vState(i, j - 1, k) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i, j - 1, k).y * factorBottom;
+ }
+
+ // v_z_front
+ vA0(i, j, k) += factorFront;
+ if (vState(i, j, k + 1) == isFluid) {
+ vAplusk(i, j, k) += -factorFront;
+ }
+ else if (vState(i, j, k + 1) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i, j, k + 1).y * factorFront;
+ }
+
+ // v_z_back
+ vA0(i, j, k) += factorBack;
+ if (vState(i, j, k - 1) == isFluid) {
+ vAminusk(i, j, k) += -factorBack;
+ }
+ else if (vState(i, j, k - 1) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i, j, k - 1).y * factorBack;
+ }
+
+ // u_y_right
+ if (uState(i + 1, j, k) == isFluid) {
+ vAhelper1(i, j, k) += -factorRight;
+ }
+ else if (uState(i + 1, j, k) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight;
+ }
+
+ if (uState(i + 1, j - 1, k) == isFluid) {
+ vAhelper2(i, j, k) += factorRight;
+ }
+ else if (uState(i + 1, j - 1, k) == isObstacle) {
+ vRhs(i, j, k) -= vel(i + 1, j - 1, k).x * factorRight;
+ }
+
+ // u_y_left
+ if (uState(i, j, k) == isFluid) {
+ vAhelper3(i, j, k) += factorLeft;
+ }
+ else if (uState(i, j, k) == isObstacle) {
+ vRhs(i, j, k) -= vel(i, j, k).x * factorLeft;
+ }
+
+ if (uState(i, j - 1, k) == isFluid) {
+ vAhelper4(i, j, k) += -factorLeft;
+ }
+ else if (uState(i, j - 1, k) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i, j - 1, k).x * factorLeft;
+ }
+
+ // w_y_front
+ if (wState(i, j, k + 1) == isFluid) {
+ vAhelper5(i, j, k) += -factorFront;
+ }
+ else if (wState(i, j, k + 1) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront;
+ }
+
+ if (wState(i, j - 1, k + 1) == isFluid) {
+ vAhelper6(i, j, k) += factorFront;
+ }
+ else if (wState(i, j - 1, k + 1) == isObstacle) {
+ vRhs(i, j, k) -= vel(i, j - 1, k + 1).z * factorFront;
+ }
+
+ // w_y_back
+ if (wState(i, j, k) == isFluid) {
+ vAhelper7(i, j, k) += factorBack;
+ }
+ else if (wState(i, j, k) == isObstacle) {
+ vRhs(i, j, k) -= vel(i, j, k).z * factorBack;
+ }
+
+ if (wState(i, j - 1, k) == isFluid) {
+ vAhelper8(i, j, k) += -factorBack;
+ }
+ else if (wState(i, j - 1, k) == isObstacle) {
+ vRhs(i, j, k) -= -vel(i, j - 1, k).z * factorBack;
+ }
+ }
+
+ // W-terms: wxx+ wyy+ 2wzz + u_zx + v_zy
+ if (wState(i, j, k) == isFluid) {
+
+ wRhs(i, j, k) = wVolLiquid(i, j, k) * vel(i, j, k).z;
+ wA0(i, j, k) = wVolLiquid(i, j, k);
+
+ Real viscRight = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) +
+ viscosity(i + 1, j, k) + viscosity(i + 1, j, k - 1));
+ Real viscLeft = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) +
+ viscosity(i - 1, j, k) + viscosity(i - 1, j, k - 1));
+ Real volRight = eyVolLiquid(i + 1, j, k);
+ Real volLeft = eyVolLiquid(i, j, k);
+
+ Real viscTop = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) + viscosity(i, j + 1, k) +
+ viscosity(i, j + 1, k - 1));
+ ;
+ Real viscBottom = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) +
+ viscosity(i, j - 1, k) + viscosity(i, j - 1, k - 1));
+ ;
+ Real volTop = exVolLiquid(i, j + 1, k);
+ Real volBottom = exVolLiquid(i, j, k);
+
+ Real viscFront = viscosity(i, j, k);
+ Real viscBack = viscosity(i, j, k - 1);
+ Real volFront = cVolLiquid(i, j, k);
+ Real volBack = cVolLiquid(i, j, k - 1);
+
+ Real factorRight = factor * viscRight * volRight;
+ Real factorLeft = factor * viscLeft * volLeft;
+ Real factorTop = factor * viscTop * volTop;
+ Real factorBottom = factor * viscBottom * volBottom;
+ Real factorFront = 2 * factor * viscFront * volFront;
+ Real factorBack = 2 * factor * viscBack * volBack;
+
+ // w_x_right
+ wA0(i, j, k) += factorRight;
+ if (wState(i + 1, j, k) == isFluid) {
+ wAplusi(i, j, k) += -factorRight;
+ }
+ else if (wState(i + 1, j, k) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i + 1, j, k).z * factorRight;
+ }
+
+ // w_x_left
+ wA0(i, j, k) += factorLeft;
+ if (wState(i - 1, j, k) == isFluid) {
+ wAminusi(i, j, k) += -factorLeft;
+ }
+ else if (wState(i - 1, j, k) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i - 1, j, k).z * factorLeft;
+ }
+
+ // w_y_top
+ wA0(i, j, k) += factorTop;
+ if (wState(i, j + 1, k) == isFluid) {
+ wAplusj(i, j, k) += -factorTop;
+ }
+ else if (wState(i, j + 1, k) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i, j + 1, k).z * factorTop;
+ }
+
+ // w_y_bottom
+ wA0(i, j, k) += factorBottom;
+ if (wState(i, j - 1, k) == isFluid) {
+ wAminusj(i, j, k) += -factorBottom;
+ }
+ else if (wState(i, j - 1, k) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i, j - 1, k).z * factorBottom;
+ }
+
+ // w_z_front
+ wA0(i, j, k) += factorFront;
+ if (wState(i, j, k + 1) == isFluid) {
+ wAplusk(i, j, k) += -factorFront;
+ }
+ else if (wState(i, j, k + 1) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront;
+ }
+
+ // w_z_back
+ wA0(i, j, k) += factorBack;
+ if (wState(i, j, k - 1) == isFluid) {
+ wAminusk(i, j, k) += -factorBack;
+ }
+ else if (wState(i, j, k - 1) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i, j, k - 1).z * factorBack;
+ }
+
+ // u_z_right
+ if (uState(i + 1, j, k) == isFluid) {
+ wAhelper1(i, j, k) += -factorRight;
+ }
+ else if (uState(i + 1, j, k) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight;
+ }
+
+ if (uState(i + 1, j, k - 1) == isFluid) {
+ wAhelper2(i, j, k) += factorRight;
+ }
+ else if (uState(i + 1, j, k - 1) == isObstacle) {
+ wRhs(i, j, k) -= vel(i + 1, j, k - 1).x * factorRight;
+ }
+
+ // u_z_left
+ if (uState(i, j, k) == isFluid) {
+ wAhelper3(i, j, k) += factorLeft;
+ }
+ else if (uState(i, j, k) == isObstacle) {
+ wRhs(i, j, k) -= vel(i, j, k).x * factorLeft;
+ }
+
+ if (uState(i, j, k - 1) == isFluid) {
+ wAhelper4(i, j, k) += -factorLeft;
+ }
+ else if (uState(i, j, k - 1) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i, j, k - 1).x * factorLeft;
+ }
+
+ // v_z_top
+ if (vState(i, j + 1, k) == isFluid) {
+ wAhelper5(i, j, k) += -factorTop;
+ }
+ else if (vState(i, j + 1, k) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop;
+ }
+
+ if (vState(i, j + 1, k - 1) == isFluid) {
+ wAhelper6(i, j, k) += factorTop;
+ }
+ else if (vState(i, j + 1, k - 1) == isObstacle) {
+ wRhs(i, j, k) -= vel(i, j + 1, k - 1).y * factorTop;
+ }
+
+ // v_z_bottom
+ if (vState(i, j, k) == isFluid) {
+ wAhelper7(i, j, k) += factorBottom;
+ }
+ else if (vState(i, j, k) == isObstacle) {
+ wRhs(i, j, k) -= vel(i, j, k).y * factorBottom;
+ }
+
+ if (vState(i, j, k - 1) == isFluid) {
+ wAhelper8(i, j, k) += -factorBottom;
+ }
+ else if (vState(i, j, k - 1) == isObstacle) {
+ wRhs(i, j, k) -= -vel(i, j, k - 1).y * factorBottom;
+ }
+ }
+ }
+
+ // CG solver for U
+ if (flags.is3D()) {
+ vector<Grid<Real> *> uMatA{&uA0,
+ &uAplusi,
+ &uAplusj,
+ &uAplusk,
+ &uAminusi,
+ &uAminusj,
+ &uAminusk,
+ &uAhelper1,
+ &uAhelper2,
+ &uAhelper3,
+ &uAhelper4,
+ &uAhelper5,
+ &uAhelper6,
+ &uAhelper7,
+ &uAhelper8};
+ vector<Grid<Real> *> uVecRhs{&vRhs, &wRhs};
+ uGcg = new GridCg<ApplyMatrixViscosityU>(
+ uSolution, uRhs, uResidual, uSearch, flags, uTmp, uMatA, uVecRhs);
+ }
+ else {
+ errMsg("2D Matrix application not yet supported in viscosity solver");
+ }
+
+ // CG solver for V
+ if (flags.is3D()) {
+ vector<Grid<Real> *> vMatA{&vA0,
+ &vAplusi,
+ &vAplusj,
+ &vAplusk,
+ &vAminusi,
+ &vAminusj,
+ &vAminusk,
+ &vAhelper1,
+ &vAhelper2,
+ &vAhelper3,
+ &vAhelper4,
+ &vAhelper5,
+ &vAhelper6,
+ &vAhelper7,
+ &vAhelper8};
+ vector<Grid<Real> *> vVecRhs{&uRhs, &wRhs};
+ vGcg = new GridCg<ApplyMatrixViscosityV>(
+ vSolution, vRhs, vResidual, vSearch, flags, vTmp, vMatA, vVecRhs);
+ }
+ else {
+ errMsg("2D Matrix application not yet supported in viscosity solver");
+ }
+
+ // CG solver for W
+ if (flags.is3D()) {
+ vector<Grid<Real> *> wMatA{&wA0,
+ &wAplusi,
+ &wAplusj,
+ &wAplusk,
+ &wAminusi,
+ &wAminusj,
+ &wAminusk,
+ &wAhelper1,
+ &wAhelper2,
+ &wAhelper3,
+ &wAhelper4,
+ &wAhelper5,
+ &wAhelper6,
+ &wAhelper7,
+ &wAhelper8};
+ vector<Grid<Real> *> wVecRhs{&uRhs, &vRhs};
+ wGcg = new GridCg<ApplyMatrixViscosityW>(
+ wSolution, wRhs, wResidual, wSearch, flags, wTmp, wMatA, wVecRhs);
+ }
+ else {
+ errMsg("2D Matrix application not yet supported in viscosity solver");
+ }
+
+ // Same accuracy for all dimensions
+ uGcg->setAccuracy(cgAccuracy);
+ vGcg->setAccuracy(cgAccuracy);
+ wGcg->setAccuracy(cgAccuracy);
+
+ // CG solve. Preconditioning not supported yet. Instead, U, V, W can optionally be solved in
+ // parallel.
+ for (int uIter = 0, vIter = 0, wIter = 0; uIter < maxIter || vIter < maxIter || wIter < maxIter;
+ uIter++, vIter++, wIter++) {
+#if ENABLE_PARALLEL == 1
+ parallel_block do_parallel
+#endif
+ if (uIter < maxIter && !uGcg->iterate()) uIter = maxIter;
+#if ENABLE_PARALLEL == 1
+ do_end do_parallel
+#endif
+ if (vIter < maxIter && !vGcg->iterate()) vIter = maxIter;
+#if ENABLE_PARALLEL == 1
+ do_end do_parallel
+#endif
+ if (wIter < maxIter && !wGcg->iterate()) wIter = maxIter;
+#if ENABLE_PARALLEL == 1
+ do_end block_end
+#endif
+
+ // Make sure that next CG iteration has updated rhs grids
+ uRhs.copyFrom(uSearch);
+ vRhs.copyFrom(vSearch);
+ wRhs.copyFrom(wSearch);
+ }
+ debMsg(
+ "Viscosity::solveViscosity done. "
+ "Iterations (u,v,w): ("
+ << uGcg->getIterations() << "," << vGcg->getIterations() << "," << wGcg->getIterations()
+ << "), "
+ "Residual (u,v,w): ("
+ << uGcg->getResNorm() << "," << vGcg->getResNorm() << "," << wGcg->getResNorm() << ")",
+ 2);
+
+ delete uGcg;
+ delete vGcg;
+ delete wGcg;
+
+ // Apply solutions to global velocity grid
+ KnApplyVelocities(vel, uState, vState, wState, uSolution, vSolution, wSolution);
+}
+
+//! To use the viscosity plugin, scenes must call this function before solving pressure.
+//! Note that the 'volumes' grid uses 2x the base resolution
+
+void applyViscosity(const FlagGrid &flags,
+ const Grid<Real> &phi,
+ MACGrid &vel,
+ Grid<Real> &volumes,
+ Grid<Real> &viscosity,
+ const Real cgAccuracy = 1e-9,
+ const Real cgMaxIterFac = 1.5)
+{
+ const Real dx = flags.getParent()->getDx();
+ const Real dt = flags.getParent()->getDt();
+
+ // Reserve temp grids for volume weight calculation
+ FluidSolver *parent = flags.getParent();
+ Grid<Real> cVolLiquid(parent);
+ Grid<Real> uVolLiquid(parent);
+ Grid<Real> vVolLiquid(parent);
+ Grid<Real> wVolLiquid(parent);
+ Grid<Real> exVolLiquid(parent);
+ Grid<Real> eyVolLiquid(parent);
+ Grid<Real> ezVolLiquid(parent);
+
+ // Ensure final weight grid gets cleared at every step
+ volumes.clear();
+
+ // Save viscous fluid volume in double-sized volumes grid
+ computeWeights(phi,
+ volumes,
+ cVolLiquid,
+ uVolLiquid,
+ vVolLiquid,
+ wVolLiquid,
+ exVolLiquid,
+ eyVolLiquid,
+ ezVolLiquid,
+ dx);
+
+ // Set up A matrix and rhs. Solve with CG. Update velocity grid.
+ solveViscosity(flags,
+ vel,
+ cVolLiquid,
+ uVolLiquid,
+ vVolLiquid,
+ wVolLiquid,
+ exVolLiquid,
+ eyVolLiquid,
+ ezVolLiquid,
+ viscosity,
+ dt,
+ dx,
+ cgAccuracy,
+ cgMaxIterFac);
+}
+static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
+{
+ try {
+ PbArgs _args(_linargs, _kwds);
+ FluidSolver *parent = _args.obtainParent();
+ bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
+ pbPreparePlugin(parent, "applyViscosity", !noTiming);
+ PyObject *_retval = nullptr;
+ {
+ ArgLocker _lock;
+ const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
+ const Grid<Real> &phi = *_args.getPtr<Grid<Real>>("phi", 1, &_lock);
+ MACGrid &vel = *_args.getPtr<MACGrid>("vel", 2, &_lock);
+ Grid<Real> &volumes = *_args.getPtr<Grid<Real>>("volumes", 3, &_lock);
+ Grid<Real> &viscosity = *_args.getPtr<Grid<Real>>("viscosity", 4, &_lock);
+ const Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 5, 1e-9, &_lock);
+ const Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 6, 1.5, &_lock);
+ _retval = getPyNone();
+ applyViscosity(flags, phi, vel, volumes, viscosity, cgAccuracy, cgMaxIterFac);
+ _args.check();
+ }
+ pbFinalizePlugin(parent, "applyViscosity", !noTiming);
+ return _retval;
+ }
+ catch (std::exception &e) {
+ pbSetError("applyViscosity", e.what());
+ return 0;
+ }
+}
+static const Pb::Register _RP_applyViscosity("", "applyViscosity", _W_0);
+extern "C" {
+void PbRegister_applyViscosity()
+{
+ KEEP_UNUSED(_RP_applyViscosity);
+}
+}
+
+} // namespace Manta
+
+#if ENABLE_PARALLEL == 1
+
+# undef parallel_block
+# undef do_parallel
+# undef do_end
+# undef block_end
+# undef parallel_for
+# undef parallel_end
+
+#endif
diff --git a/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp b/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp
index d5d7d597a7c..18222c4ccda 100644
--- a/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp
+++ b/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp
@@ -576,8 +576,10 @@ void VICintegration(VortexSheetMesh &mesh,
// prepare CG solver
const int maxIter = (int)(cgMaxIterFac * vel.getSize().max());
+ vector<Grid<Real> *> matA{&A0, &Ai, &Aj, &Ak};
+
GridCgInterface *gcg = new GridCg<ApplyMatrix>(
- solution, rhs, residual, search, flags, temp1, &A0, &Ai, &Aj, &Ak);
+ solution, rhs, residual, search, flags, temp1, matA);
gcg->setAccuracy(cgAccuracy);
gcg->setUseL2Norm(true);
gcg->setICPreconditioner(
diff --git a/extern/mantaflow/preprocessed/plugin/waves.cpp b/extern/mantaflow/preprocessed/plugin/waves.cpp
index 6fb3b16c742..53c56b8c506 100644
--- a/extern/mantaflow/preprocessed/plugin/waves.cpp
+++ b/extern/mantaflow/preprocessed/plugin/waves.cpp
@@ -423,10 +423,15 @@ void cgSolveWE(const FlagGrid &flags,
const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4);
GridCgInterface *gcg;
- if (flags.is3D())
- gcg = new GridCg<ApplyMatrix>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
- else
- gcg = new GridCg<ApplyMatrix2D>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+ vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
+
+ if (flags.is3D()) {
+ matA.push_back(&Ak);
+ gcg = new GridCg<ApplyMatrix>(out, rhs, residual, search, flags, tmp, matA);
+ }
+ else {
+ gcg = new GridCg<ApplyMatrix2D>(out, rhs, residual, search, flags, tmp, matA);
+ }
gcg->setAccuracy(cgAccuracy);
diff --git a/extern/mantaflow/preprocessed/registration.cpp b/extern/mantaflow/preprocessed/registration.cpp
index d5dae479f0e..fd32463b95f 100644
--- a/extern/mantaflow/preprocessed/registration.cpp
+++ b/extern/mantaflow/preprocessed/registration.cpp
@@ -145,6 +145,7 @@ extern void PbRegister_flipComputeSurfaceNormals();
extern void PbRegister_flipUpdateNeighborRatio();
extern void PbRegister_particleSurfaceTurbulence();
extern void PbRegister_debugCheckParts();
+extern void PbRegister_applyViscosity();
extern void PbRegister_markAsFixed();
extern void PbRegister_texcoordInflow();
extern void PbRegister_meshSmokeInflow();
@@ -342,6 +343,7 @@ void MantaEnsureRegistration()
PbRegister_flipUpdateNeighborRatio();
PbRegister_particleSurfaceTurbulence();
PbRegister_debugCheckParts();
+ PbRegister_applyViscosity();
PbRegister_markAsFixed();
PbRegister_texcoordInflow();
PbRegister_meshSmokeInflow();
diff --git a/intern/mantaflow/intern/MANTA_main.cpp b/intern/mantaflow/intern/MANTA_main.cpp
index d49915e2452..d5d41e71f22 100644
--- a/intern/mantaflow/intern/MANTA_main.cpp
+++ b/intern/mantaflow/intern/MANTA_main.cpp
@@ -72,6 +72,7 @@ MANTA::MANTA(int *res, FluidModifierData *fmd)
mUsingFractions = (fds->flags & FLUID_DOMAIN_USE_FRACTIONS) && mUsingLiquid;
mUsingMesh = (fds->flags & FLUID_DOMAIN_USE_MESH) && mUsingLiquid;
mUsingDiffusion = (fds->flags & FLUID_DOMAIN_USE_DIFFUSION) && mUsingLiquid;
+ mUsingViscosity = (fds->flags & FLUID_DOMAIN_USE_VISCOSITY) && mUsingLiquid;
mUsingMVel = (fds->flags & FLUID_DOMAIN_USE_SPEED_VECTORS) && mUsingLiquid;
mUsingGuiding = (fds->flags & FLUID_DOMAIN_USE_GUIDE);
mUsingDrops = (fds->particle_type & FLUID_DOMAIN_PARTICLE_SPRAY) && mUsingLiquid;
@@ -221,6 +222,10 @@ MANTA::MANTA(int *res, FluidModifierData *fmd)
initSuccess &= initLiquidMesh();
}
+ if (mUsingViscosity) {
+ initSuccess &= initLiquidViscosity();
+ }
+
if (mUsingDiffusion) {
initSuccess &= initCurvature();
}
@@ -440,6 +445,17 @@ bool MANTA::initLiquidMesh(FluidModifierData *fmd)
return runPythonString(pythonCommands);
}
+bool MANTA::initLiquidViscosity(FluidModifierData *fmd)
+{
+ vector<string> pythonCommands;
+ string tmpString = fluid_variables_viscosity + fluid_solver_viscosity + liquid_alloc_viscosity;
+ string finalString = parseScript(tmpString, fmd);
+ pythonCommands.push_back(finalString);
+
+ mUsingViscosity = true;
+ return runPythonString(pythonCommands);
+}
+
bool MANTA::initCurvature(FluidModifierData *fmd)
{
std::vector<std::string> pythonCommands;
@@ -871,8 +887,10 @@ void MANTA::initializeRNAMap(FluidModifierData *fmd)
mRNAMap["CACHE_DIR"] = cacheDirectory;
mRNAMap["COMPRESSION_OPENVDB"] = vdbCompressionMethod;
mRNAMap["PRECISION_OPENVDB"] = vdbPrecisionHalf;
- mRNAMap["CLIP_OPENVDB"] = to_string(fds->clipping);
+ mRNAMap["CLIP_OPENVDB"] = to_string(fds->clipping);
mRNAMap["PP_PARTICLE_MAXIMUM"] = to_string(fds->sys_particle_maximum);
+ mRNAMap["USING_VISCOSITY"] = getBooleanString(fds->flags & FLUID_DOMAIN_USE_VISCOSITY);
+ mRNAMap["VISCOSITY_VALUE"] = to_string(fds->viscosity_value);
/* Fluid object names. */
mRNAMap["NAME_FLAGS"] = FLUID_NAME_FLAGS;
@@ -1728,6 +1746,7 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
bool guiding = fds->active_fields & FLUID_DOMAIN_ACTIVE_GUIDE;
bool invel = fds->active_fields & FLUID_DOMAIN_ACTIVE_INVEL;
bool outflow = fds->active_fields & FLUID_DOMAIN_ACTIVE_OUTFLOW;
+ bool viscosity = fds->flags & FLUID_DOMAIN_USE_VISCOSITY;
string manta_script;
@@ -1742,6 +1761,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
manta_script += fluid_variables_particles + liquid_variables_particles;
if (guiding)
manta_script += fluid_variables_guiding;
+ if (viscosity)
+ manta_script += fluid_variables_viscosity;
/* Solvers. */
manta_script += header_solvers + fluid_solver;
@@ -1751,6 +1772,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
manta_script += fluid_solver_particles;
if (guiding)
manta_script += fluid_solver_guiding;
+ if (viscosity)
+ manta_script += fluid_solver_viscosity;
/* Grids. */
manta_script += header_grids + fluid_alloc + liquid_alloc;
@@ -1768,6 +1791,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
manta_script += fluid_alloc_invel;
if (outflow)
manta_script += fluid_alloc_outflow;
+ if (viscosity)
+ manta_script += liquid_alloc_viscosity;
/* Domain init. */
manta_script += header_gridinit + liquid_init_phi;
diff --git a/intern/mantaflow/intern/MANTA_main.h b/intern/mantaflow/intern/MANTA_main.h
index 163b168e43d..7129a545243 100644
--- a/intern/mantaflow/intern/MANTA_main.h
+++ b/intern/mantaflow/intern/MANTA_main.h
@@ -67,6 +67,7 @@ struct MANTA {
bool initColorsHigh(struct FluidModifierData *fmd = nullptr);
bool initLiquid(FluidModifierData *fmd = nullptr);
bool initLiquidMesh(FluidModifierData *fmd = nullptr);
+ bool initLiquidViscosity(FluidModifierData *fmd = nullptr);
bool initObstacle(FluidModifierData *fmd = nullptr);
bool initCurvature(FluidModifierData *fmd = nullptr);
bool initGuiding(FluidModifierData *fmd = nullptr);
@@ -757,6 +758,7 @@ struct MANTA {
bool mUsingNoise;
bool mUsingMesh;
bool mUsingDiffusion;
+ bool mUsingViscosity;
bool mUsingMVel;
bool mUsingLiquid;
bool mUsingSmoke;
diff --git a/intern/mantaflow/intern/strings/fluid_script.h b/intern/mantaflow/intern/strings/fluid_script.h
index 8c9025dd435..73e8552d260 100644
--- a/intern/mantaflow/intern/strings/fluid_script.h
+++ b/intern/mantaflow/intern/strings/fluid_script.h
@@ -79,6 +79,11 @@ const std::string fluid_solver_guiding =
mantaMsg('Solver guiding')\n\
sg$ID$ = Solver(name='solver_guiding$ID$', gridSize=gs_sg$ID$)\n";
+const std::string fluid_solver_viscosity =
+ "\n\
+mantaMsg('Solver viscosity')\n\
+sv$ID$ = Solver(name='solver_viscosity$ID$', gridSize=gs_sv$ID$, dim=dim_s$ID$)\n";
+
//////////////////////////////////////////////////////////////////////
// VARIABLES
//////////////////////////////////////////////////////////////////////
@@ -133,7 +138,7 @@ end_frame_s$ID$ = $END_FRAME$\n\
\n\
# Fluid diffusion / viscosity\n\
domainSize_s$ID$ = $FLUID_DOMAIN_SIZE$ # longest domain side in meters\n\
-viscosity_s$ID$ = $FLUID_VISCOSITY$ / (domainSize_s$ID$*domainSize_s$ID$) # kinematic viscosity in m^2/s\n\
+kinViscosity_s$ID$ = $FLUID_VISCOSITY$ / (domainSize_s$ID$*domainSize_s$ID$) # kinematic viscosity in m^2/s\n\
\n\
# Factors to convert Blender units to Manta units\n\
ratioMetersToRes_s$ID$ = float(domainSize_s$ID$) / float(res_s$ID$) # [meters / cells]\n\
@@ -199,6 +204,10 @@ tau_sg$ID$ = 1.0\n\
sigma_sg$ID$ = 0.99/tau_sg$ID$\n\
theta_sg$ID$ = 1.0\n";
+const std::string fluid_variables_viscosity =
+ "\n\
+gs_sv$ID$ = vec3($RESX$*2, $RESY$*2, $RESZ$*2)\n";
+
const std::string fluid_with_obstacle =
"\n\
using_obstacle_s$ID$ = True\n";
diff --git a/intern/mantaflow/intern/strings/liquid_script.h b/intern/mantaflow/intern/strings/liquid_script.h
index 8e0d113d680..c44727bd47e 100644
--- a/intern/mantaflow/intern/strings/liquid_script.h
+++ b/intern/mantaflow/intern/strings/liquid_script.h
@@ -40,6 +40,8 @@ radiusFactor_s$ID$ = $PARTICLE_RADIUS$\n\
using_mesh_s$ID$ = $USING_MESH$\n\
using_final_mesh_s$ID$ = $USING_IMPROVED_MESH$\n\
using_fractions_s$ID$ = $USING_FRACTIONS$\n\
+using_apic_s$ID$ = $USING_APIC$\n\
+using_viscosity_s$ID$ = $USING_VISCOSITY$\n\
fracThreshold_s$ID$ = $FRACTIONS_THRESHOLD$\n\
fracDistance_s$ID$ = $FRACTIONS_DISTANCE$\n\
flipRatio_s$ID$ = $FLIP_RATIO$\n\
@@ -51,7 +53,7 @@ smoothenNeg_s$ID$ = $MESH_SMOOTHEN_NEG$\n\
randomness_s$ID$ = $PARTICLE_RANDOMNESS$\n\
surfaceTension_s$ID$ = $LIQUID_SURFACE_TENSION$\n\
maxSysParticles_s$ID$ = $PP_PARTICLE_MAXIMUM$\n\
-using_apic_s$ID$ = $USING_APIC$\n";
+viscosityValue_s$ID$ = $VISCOSITY_VALUE$\n";
const std::string liquid_variables_particles =
"\n\
@@ -135,6 +137,13 @@ liquid_mesh_dict_s$ID$ = { 'lMesh' : mesh_sm$ID$ }\n\
if using_speedvectors_s$ID$:\n\
liquid_meshvel_dict_s$ID$ = { 'lVelMesh' : mVel_mesh$ID$ }\n";
+const std::string liquid_alloc_viscosity =
+ "\n\
+# Viscosity grids\n\
+volumes_s$ID$ = sv$ID$.create(RealGrid)\n\
+viscosity_s$ID$ = s$ID$.create(RealGrid)\n\
+viscosity_s$ID$.setConst(viscosityValue_s$ID$)\n";
+
const std::string liquid_alloc_curvature =
"\n\
mantaMsg('Liquid alloc curvature')\n\
@@ -306,7 +315,7 @@ def liquid_step_$ID$():\n\
if using_diffusion_s$ID$:\n\
mantaMsg('Viscosity')\n\
# diffusion param for solve = const * dt / dx^2\n\
- alphaV = viscosity_s$ID$ * s$ID$.timestep * float(res_s$ID$*res_s$ID$)\n\
+ alphaV = kinViscosity_s$ID$ * s$ID$.timestep * float(res_s$ID$*res_s$ID$)\n\
setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\
cgSolveDiffusion(flags_s$ID$, vel_s$ID$, alphaV)\n\
\n\
@@ -315,7 +324,11 @@ def liquid_step_$ID$():\n\
curvature_s$ID$.clamp(-1.0, 1.0)\n\
\n\
setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\
+ if using_viscosity_s$ID$:\n\
+ viscosity_s$ID$.setConst(viscosityValue_s$ID$)\n\
+ applyViscosity(flags=flags_s$ID$, phi=phi_s$ID$, vel=vel_s$ID$, volumes=volumes_s$ID$, viscosity=viscosity_s$ID$)\n\
\n\
+ setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\
if using_guiding_s$ID$:\n\
mantaMsg('Guiding and pressure')\n\
PD_fluid_guiding(vel=vel_s$ID$, velT=velT_s$ID$, flags=flags_s$ID$, phi=phi_s$ID$, curv=curvature_s$ID$, surfTens=surfaceTension_s$ID$, fractions=fractions_s$ID$, weight=weightGuide_s$ID$, blurRadius=beta_sg$ID$, pressure=pressure_s$ID$, tau=tau_sg$ID$, sigma=sigma_sg$ID$, theta=theta_sg$ID$, zeroPressureFixing=domainClosed_s$ID$)\n\
diff --git a/release/scripts/startup/bl_ui/properties_physics_fluid.py b/release/scripts/startup/bl_ui/properties_physics_fluid.py
index 3afe7a47028..5c3975ffb76 100644
--- a/release/scripts/startup/bl_ui/properties_physics_fluid.py
+++ b/release/scripts/startup/bl_ui/properties_physics_fluid.py
@@ -1006,6 +1006,46 @@ class PHYSICS_PT_particles(PhysicButtonsPanel, Panel):
split.operator("fluid.free_particles", text="Free Particles")
+class PHYSICS_PT_viscosity(PhysicButtonsPanel, Panel):
+ bl_label = "Viscosity"
+ bl_parent_id = 'PHYSICS_PT_liquid'
+ bl_options = {'DEFAULT_CLOSED'}
+ COMPAT_ENGINES = {'BLENDER_RENDER', 'BLENDER_EEVEE', 'BLENDER_WORKBENCH'}
+
+ @classmethod
+ def poll(cls, context):
+ # Fluid viscosity only enabled for liquids
+ if not PhysicButtonsPanel.poll_liquid_domain(context):
+ return False
+
+ return (context.engine in cls.COMPAT_ENGINES)
+
+ def draw_header(self, context):
+ md = context.fluid.domain_settings
+ domain = context.fluid.domain_settings
+ is_baking_any = domain.is_cache_baking_any
+ has_baked_any = domain.has_cache_baked_any
+ self.layout.enabled = not is_baking_any and not has_baked_any
+ self.layout.prop(md, "use_viscosity", text="")
+
+ def draw(self, context):
+ layout = self.layout
+ layout.use_property_split = True
+
+ domain = context.fluid.domain_settings
+ layout.active = domain.use_viscosity
+
+ is_baking_any = domain.is_cache_baking_any
+ has_baked_any = domain.has_cache_baked_any
+ has_baked_data = domain.has_cache_baked_data
+
+ flow = layout.grid_flow(row_major=True, columns=0, even_columns=True, even_rows=False, align=False)
+ flow.enabled = not is_baking_any and not has_baked_any and not has_baked_data
+
+ col = flow.column(align=True)
+ col.prop(domain, "viscosity_value", text="Strength")
+
+
class PHYSICS_PT_diffusion(PhysicButtonsPanel, Panel):
bl_label = "Diffusion"
bl_parent_id = 'PHYSICS_PT_liquid'
@@ -1470,6 +1510,7 @@ classes = (
PHYSICS_PT_noise,
PHYSICS_PT_fire,
PHYSICS_PT_liquid,
+ PHYSICS_PT_viscosity,
PHYSICS_PT_diffusion,
PHYSICS_PT_particles,
PHYSICS_PT_mesh,
diff --git a/source/blender/blenkernel/BKE_blender_version.h b/source/blender/blenkernel/BKE_blender_version.h
index 1ed4d1183a1..06a11fc8d1d 100644
--- a/source/blender/blenkernel/BKE_blender_version.h
+++ b/source/blender/blenkernel/BKE_blender_version.h
@@ -39,7 +39,7 @@ extern "C" {
/* Blender file format version. */
#define BLENDER_FILE_VERSION BLENDER_VERSION
-#define BLENDER_FILE_SUBVERSION 8
+#define BLENDER_FILE_SUBVERSION 9
/* Minimum Blender version that supports reading file written with the current
* version. Older Blender versions will test this and show a warning if the file
diff --git a/source/blender/blenkernel/intern/fluid.c b/source/blender/blenkernel/intern/fluid.c
index 1bdc2a5dfd2..b93e6148b86 100644
--- a/source/blender/blenkernel/intern/fluid.c
+++ b/source/blender/blenkernel/intern/fluid.c
@@ -5011,6 +5011,9 @@ void BKE_fluid_modifier_copy(const struct FluidModifierData *fmd,
tfds->sys_particle_maximum = fds->sys_particle_maximum;
tfds->simulation_method = fds->simulation_method;
+ /* viscosity options */
+ tfds->viscosity_value = fds->viscosity_value;
+
/* diffusion options*/
tfds->surface_tension = fds->surface_tension;
tfds->viscosity_base = fds->viscosity_base;
diff --git a/source/blender/blenloader/intern/versioning_290.c b/source/blender/blenloader/intern/versioning_290.c
index 8a85278a931..63ef436d8e2 100644
--- a/source/blender/blenloader/intern/versioning_290.c
+++ b/source/blender/blenloader/intern/versioning_290.c
@@ -1430,18 +1430,7 @@ void blo_do_versions_290(FileData *fd, Library *UNUSED(lib), Main *bmain)
}
}
- /**
- * Versioning code until next subversion bump goes here.
- *
- * \note Be sure to check when bumping the version:
- * - "versioning_userdef.c", #blo_do_versions_userdef
- * - "versioning_userdef.c", #do_versions_theme
- *
- * \note Keep this message at the bottom of the function.
- */
- {
- /* Keep this block, even when empty. */
-
+ if (!MAIN_VERSION_ATLEAST(bmain, 292, 9)) {
FOREACH_NODETREE_BEGIN (bmain, ntree, id) {
if (ntree->type == NTREE_GEOMETRY) {
LISTBASE_FOREACH (bNode *, node, &ntree->nodes) {
@@ -1474,5 +1463,32 @@ void blo_do_versions_290(FileData *fd, Library *UNUSED(lib), Main *bmain)
}
}
}
+
+ /* Ensure that new viscosity strength field is initialized correctly. */
+ if (!DNA_struct_elem_find(fd->filesdna, "FluidModifierData", "float", "viscosity_value")) {
+ LISTBASE_FOREACH (Object *, ob, &bmain->objects) {
+ LISTBASE_FOREACH (ModifierData *, md, &ob->modifiers) {
+ if (md->type == eModifierType_Fluid) {
+ FluidModifierData *fmd = (FluidModifierData *)md;
+ if (fmd->domain != NULL) {
+ fmd->domain->viscosity_value = 0.05;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /**
+ * Versioning code until next subversion bump goes here.
+ *
+ * \note Be sure to check when bumping the version:
+ * - "versioning_userdef.c", #blo_do_versions_userdef
+ * - "versioning_userdef.c", #do_versions_theme
+ *
+ * \note Keep this message at the bottom of the function.
+ */
+ {
+ /* Keep this block, even when empty. */
}
}
diff --git a/source/blender/blenloader/intern/versioning_userdef.c b/source/blender/blenloader/intern/versioning_userdef.c
index f1572c563bf..b9f09e0326d 100644
--- a/source/blender/blenloader/intern/versioning_userdef.c
+++ b/source/blender/blenloader/intern/versioning_userdef.c
@@ -820,6 +820,12 @@ void blo_do_versions_userdef(UserDef *userdef)
userdef->uiflag &= ~USER_UIFLAG_UNUSED_3;
}
+ if (!USER_VERSION_ATLEAST(292, 9)) {
+ if (BLI_listbase_is_empty(&userdef->asset_libraries)) {
+ BKE_preferences_asset_library_default_add(userdef);
+ }
+ }
+
/**
* Versioning code until next subversion bump goes here.
*
@@ -831,9 +837,6 @@ void blo_do_versions_userdef(UserDef *userdef)
*/
{
/* Keep this block, even when empty. */
- if (BLI_listbase_is_empty(&userdef->asset_libraries)) {
- BKE_preferences_asset_library_default_add(userdef);
- }
}
LISTBASE_FOREACH (bTheme *, btheme, &userdef->themes) {
diff --git a/source/blender/makesdna/DNA_fluid_defaults.h b/source/blender/makesdna/DNA_fluid_defaults.h
index 2ee83cf3387..9454342654c 100644
--- a/source/blender/makesdna/DNA_fluid_defaults.h
+++ b/source/blender/makesdna/DNA_fluid_defaults.h
@@ -113,6 +113,7 @@
.flip_ratio = 0.97f, \
.sys_particle_maximum = 0, \
.simulation_method = FLUID_DOMAIN_METHOD_FLIP, \
+ .viscosity_value = 0.05f, \
.surface_tension = 0.0f, \
.viscosity_base = 1.0f, \
.viscosity_exponent = 6.0f, \
diff --git a/source/blender/makesdna/DNA_fluid_types.h b/source/blender/makesdna/DNA_fluid_types.h
index 98007a2b0b1..2ba0d15fbb3 100644
--- a/source/blender/makesdna/DNA_fluid_types.h
+++ b/source/blender/makesdna/DNA_fluid_types.h
@@ -52,6 +52,7 @@ enum {
FLUID_DOMAIN_DELETE_IN_OBSTACLE = (1 << 14), /* Delete fluid inside obstacles. */
FLUID_DOMAIN_USE_DIFFUSION = (1 << 15), /* Use diffusion (e.g. viscosity, surface tension). */
FLUID_DOMAIN_USE_RESUMABLE_CACHE = (1 << 16), /* Determine if cache should be resumable. */
+ FLUID_DOMAIN_USE_VISCOSITY = (1 << 17), /* Use viscosity. */
};
/**
@@ -289,7 +290,7 @@ enum {
#define FLUID_NAME_GUIDING "fluid_guiding"
/* Fluid object names.*/
-#define FLUID_NAME_FLAGS "flags" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_FLAGS "flags" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_VELOCITY "velocity" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_VEL "vel"
#define FLUID_NAME_VELOCITYTMP "velocity_previous" /* == OpenVDB grid attribute name. */
@@ -300,7 +301,7 @@ enum {
#define FLUID_NAME_PHIOBS "phi_obstacle" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_PHISIN "phiSIn"
#define FLUID_NAME_PHIIN "phi_inflow" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_PHIOUT "phi_out" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_PHIOUT "phi_out" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_FORCES "forces"
#define FLUID_NAME_FORCE_X "x_force"
#define FLUID_NAME_FORCE_Y "y_force"
@@ -322,37 +323,37 @@ enum {
#define FLUID_NAME_PHIOUTIN "phi_out_inflow"
/* Smoke object names. */
-#define FLUID_NAME_SHADOW "shadow" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_SHADOW "shadow" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_EMISSION "emission" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_EMISSIONIN "emissionIn"
-#define FLUID_NAME_DENSITY "density" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_DENSITY "density" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_DENSITYIN "density_inflow" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_HEAT "heat"
#define FLUID_NAME_HEATIN "heatIn"
-#define FLUID_NAME_TEMPERATURE "temperature" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_TEMPERATURE "temperature" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_TEMPERATUREIN "temperature_inflow" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORR "color_r" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORG "color_g" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORB "color_b" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORRIN "color_r_inflow" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORGIN "color_g_inflow" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORBIN "color_b_inflow" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_FLAME "flame" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_FUEL "fuel" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_REACT "react" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_FUELIN "fuel_inflow" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_REACTIN "react_inflow" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_COLORR "color_r" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_COLORG "color_g" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_COLORB "color_b" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_COLORRIN "color_r_inflow" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_COLORGIN "color_g_inflow" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_COLORBIN "color_b_inflow" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_FLAME "flame" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_FUEL "fuel" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_REACT "react" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_FUELIN "fuel_inflow" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_REACTIN "react_inflow" /* == OpenVDB grid attribute name. */
/* Liquid object names. */
#define FLUID_NAME_PHIPARTS "phi_particles" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_PHI "phi" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_PHITMP "phi_previous" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_PHI "phi" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_PHITMP "phi_previous" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_VELOCITYOLD "velOld"
#define FLUID_NAME_VELOCITYPARTS "velParts"
#define FLUID_NAME_MAPWEIGHTS "mapWeights"
#define FLUID_NAME_PP "pp"
#define FLUID_NAME_PVEL "pVel"
-#define FLUID_NAME_PARTS "particles" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_PARTS "particles" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_PARTSVELOCITY "particles_velocity" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_PINDEX "pindex"
#define FLUID_NAME_GPI "gpi"
@@ -375,8 +376,8 @@ enum {
#define FLUID_NAME_TEXTURE_U2 "textureU2"
#define FLUID_NAME_TEXTURE_V2 "textureV2"
#define FLUID_NAME_TEXTURE_W2 "textureW2"
-#define FLUID_NAME_UV0 "uv_grid_0" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_UV1 "uv_grid_1" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_UV0 "uv_grid_0" /* == OpenVDB grid attribute name. */
+#define FLUID_NAME_UV1 "uv_grid_1" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_COLORR_NOISE "color_r_noise" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_COLORG_NOISE "color_g_noise" /* == OpenVDB grid attribute name. */
#define FLUID_NAME_COLORB_NOISE "color_b_noise" /* == OpenVDB grid attribute name. */
@@ -596,6 +597,10 @@ typedef struct FluidDomainSettings {
short simulation_method;
char _pad4[6];
+ /* Viscosity options. */
+ float viscosity_value;
+ char _pad5[4];
+
/* Diffusion options. */
float surface_tension;
float viscosity_base;
@@ -610,7 +615,7 @@ typedef struct FluidDomainSettings {
int mesh_scale;
int totvert;
short mesh_generator;
- char _pad5[6]; /* Unused. */
+ char _pad6[6]; /* Unused. */
/* Secondary particle options. */
int particle_type;
@@ -631,7 +636,7 @@ typedef struct FluidDomainSettings {
int sndparticle_update_radius;
char sndparticle_boundary;
char sndparticle_combined_export;
- char _pad6[6]; /* Unused. */
+ char _pad7[6]; /* Unused. */
/* Fluid guiding options. */
float guide_alpha; /* Guiding weight scalar (determines strength). */
@@ -639,7 +644,7 @@ typedef struct FluidDomainSettings {
float guide_vel_factor; /* Multiply guiding velocity by this factor. */
int guide_res[3]; /* Res for velocity guide grids - independent from base res. */
short guide_source;
- char _pad7[2]; /* Unused. */
+ char _pad8[2]; /* Unused. */
/* Cache options. */
int cache_frame_start;
@@ -659,7 +664,7 @@ typedef struct FluidDomainSettings {
char error[64]; /* Bake error description. */
short cache_type;
char cache_id[4]; /* Run-time only */
- char _pad8[2];
+ char _pad9[2]; /* Unused. */
/* Time options. */
float dt;
@@ -694,19 +699,19 @@ typedef struct FluidDomainSettings {
char interp_method;
char gridlines_color_field; /* Simulation field used to color map onto gridlines. */
char gridlines_cell_filter;
- char _pad9[7];
+ char _pad10[7]; /* Unused. */
/* OpenVDB cache options. */
int openvdb_compression;
float clipping;
char openvdb_data_depth;
- char _pad10[7]; /* Unused. */
+ char _pad11[7]; /* Unused. */
/* -- Deprecated / unsed options (below). -- */
/* View options. */
int viewsettings;
- char _pad11[4]; /* Unused. */
+ char _pad12[4]; /* Unused. */
/* Pointcache options. */
/* Smoke uses only one cache from now on (index [0]), but keeping the array for now for reading
@@ -716,7 +721,7 @@ typedef struct FluidDomainSettings {
int cache_comp;
int cache_high_comp;
char cache_file_format;
- char _pad12[7]; /* Unused. */
+ char _pad13[7]; /* Unused. */
} FluidDomainSettings;
diff --git a/source/blender/makesrna/intern/rna_fluid.c b/source/blender/makesrna/intern/rna_fluid.c
index bb8280ede91..afe564dff0a 100644
--- a/source/blender/makesrna/intern/rna_fluid.c
+++ b/source/blender/makesrna/intern/rna_fluid.c
@@ -1929,6 +1929,22 @@ static void rna_def_fluid_domain_settings(BlenderRNA *brna)
"Maximum number of fluid particles that are allowed in this simulation");
RNA_def_property_update(prop, NC_OBJECT | ND_MODIFIER, "rna_Fluid_datacache_reset");
+ /* viscosity options */
+
+ prop = RNA_def_property(srna, "use_viscosity", PROP_BOOLEAN, PROP_NONE);
+ RNA_def_property_boolean_sdna(prop, NULL, "flags", FLUID_DOMAIN_USE_VISCOSITY);
+ RNA_def_property_ui_text(prop, "Use Viscosity", "Enable fluid viscosity settings");
+ RNA_def_property_update(prop, NC_OBJECT | ND_MODIFIER, "rna_Fluid_datacache_reset");
+
+ prop = RNA_def_property(srna, "viscosity_value", PROP_FLOAT, PROP_NONE);
+ RNA_def_property_range(prop, 0.0, 10.0);
+ RNA_def_property_ui_range(prop, 0.0, 5.0, 0.01, 3);
+ RNA_def_property_ui_text(prop,
+ "Strength",
+ "Viscosity of liquid (higher values result in more viscous fluids, a "
+ "value of 0 will still apply some viscosity)");
+ RNA_def_property_update(prop, NC_OBJECT | ND_MODIFIER, "rna_Fluid_datacache_reset");
+
/* diffusion options */
prop = RNA_def_property(srna, "use_diffusion", PROP_BOOLEAN, PROP_NONE);
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 03:08:46 +0300