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

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/extern
diff options
context:
space:
mode:
Diffstat (limited to 'extern')
-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
12 files changed, 1965 insertions, 123 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();
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-10 18:56:03 +0300