Changeset View
Changeset View
Standalone View
Standalone View
extern/mantaflow/preprocessed/plugin/fluidguiding.cpp
- This file was added.
| // DO NOT EDIT ! | |||||
| // This file is generated using the MantaFlow preprocessor (prep generate). | |||||
| /****************************************************************************** | |||||
| * | |||||
| * MantaFlow fluid solver framework | |||||
| * Copyright 2011 Tobias Pfaff, Nils Thuerey | |||||
| * | |||||
| * This program is free software, distributed under the terms of the | |||||
| * Apache License, Version 2.0 | |||||
| * http://www.apache.org/licenses/LICENSE-2.0 | |||||
| * | |||||
| * Plugins for pressure correction: solve_pressure, and ghost fluid helpers | |||||
| * | |||||
| ******************************************************************************/ | |||||
| #include "vectorbase.h" | |||||
| #include "grid.h" | |||||
| #include "kernel.h" | |||||
| #include "conjugategrad.h" | |||||
| #include "rcmatrix.h" | |||||
| using namespace std; | |||||
| namespace Manta { | |||||
| // only supports a single blur size for now, globals stored here | |||||
| bool gBlurPrecomputed = false; | |||||
| int gBlurKernelRadius = -1; | |||||
| Matrix gBlurKernel; | |||||
| // ***************************************************************************** | |||||
| // Helper functions for fluid guiding | |||||
| //! creates a 1D (horizontal) Gaussian blur kernel of size n and standard deviation sigma | |||||
| Matrix get1DGaussianBlurKernel(const int n, const int sigma) | |||||
| { | |||||
| Matrix x(n), y(n); | |||||
| for (int j = 0; j < n; j++) { | |||||
| x.add_to_element(0, j, -(n - 1) * 0.5); | |||||
| y.add_to_element(0, j, j - (n - 1) * 0.5); | |||||
| } | |||||
| Matrix G(n); | |||||
| Real sumG = 0; | |||||
| for (int j = 0; j < n; j++) { | |||||
| G.add_to_element(0, | |||||
| j, | |||||
| 1 / (2 * M_PI * sigma * sigma) * | |||||
| exp(-(x(0, j) * x(0, j) + y(0, j) * y(0, j)) / (2 * sigma * sigma))); | |||||
| sumG += G(0, j); | |||||
| } | |||||
| G = G * (1.0 / sumG); | |||||
| return G; | |||||
| } | |||||
| //! convolves in with 1D kernel (centred at the kernel's midpoint) in the x-direction | |||||
| //! (out must be a grid of zeros) | |||||
| struct apply1DKernelDirX : public KernelBase { | |||||
| apply1DKernelDirX(const MACGrid &in, MACGrid &out, const Matrix &kernel) | |||||
| : KernelBase(&in, 0), in(in), out(out), kernel(kernel) | |||||
| { | |||||
| runMessage(); | |||||
| run(); | |||||
| } | |||||
| inline void op(int i, int j, int k, const MACGrid &in, MACGrid &out, const Matrix &kernel) const | |||||
| { | |||||
| int nx = in.getSizeX(); | |||||
| int kn = kernel.n; | |||||
| int kCentre = kn / 2; | |||||
| for (int m = 0, ind = kn - 1, ii = i - kCentre; m < kn; m++, ind--, ii++) { | |||||
| if (ii < 0) | |||||
| continue; | |||||
| else if (ii >= nx) | |||||
| break; | |||||
| else | |||||
| out(i, j, k) += in(ii, j, k) * kernel(0, ind); | |||||
| } | |||||
| } | |||||
| inline const MACGrid &getArg0() | |||||
| { | |||||
| return in; | |||||
| } | |||||
| typedef MACGrid type0; | |||||
| inline MACGrid &getArg1() | |||||
| { | |||||
| return out; | |||||
| } | |||||
| typedef MACGrid type1; | |||||
| inline const Matrix &getArg2() | |||||
| { | |||||
| return kernel; | |||||
| } | |||||
| typedef Matrix type2; | |||||
| void runMessage() | |||||
| { | |||||
| debMsg("Executing kernel apply1DKernelDirX ", 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, in, out, kernel); | |||||
| } | |||||
| 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, in, out, kernel); | |||||
| } | |||||
| } | |||||
| 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 MACGrid ∈ | |||||
| MACGrid &out; | |||||
| const Matrix &kernel; | |||||
| }; | |||||
| //! convolves in with 1D kernel (centred at the kernel's midpoint) in the y-direction | |||||
| //! (out must be a grid of zeros) | |||||
| struct apply1DKernelDirY : public KernelBase { | |||||
| apply1DKernelDirY(const MACGrid &in, MACGrid &out, const Matrix &kernel) | |||||
| : KernelBase(&in, 0), in(in), out(out), kernel(kernel) | |||||
| { | |||||
| runMessage(); | |||||
| run(); | |||||
| } | |||||
| inline void op(int i, int j, int k, const MACGrid &in, MACGrid &out, const Matrix &kernel) const | |||||
| { | |||||
| int ny = in.getSizeY(); | |||||
| int kn = kernel.n; | |||||
| int kCentre = kn / 2; | |||||
| for (int m = 0, ind = kn - 1, jj = j - kCentre; m < kn; m++, ind--, jj++) { | |||||
| if (jj < 0) | |||||
| continue; | |||||
| else if (jj >= ny) | |||||
| break; | |||||
| else | |||||
| out(i, j, k) += in(i, jj, k) * kernel(0, ind); | |||||
| } | |||||
| } | |||||
| inline const MACGrid &getArg0() | |||||
| { | |||||
| return in; | |||||
| } | |||||
| typedef MACGrid type0; | |||||
| inline MACGrid &getArg1() | |||||
| { | |||||
| return out; | |||||
| } | |||||
| typedef MACGrid type1; | |||||
| inline const Matrix &getArg2() | |||||
| { | |||||
| return kernel; | |||||
| } | |||||
| typedef Matrix type2; | |||||
| void runMessage() | |||||
| { | |||||
| debMsg("Executing kernel apply1DKernelDirY ", 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, in, out, kernel); | |||||
| } | |||||
| 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, in, out, kernel); | |||||
| } | |||||
| } | |||||
| 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 MACGrid ∈ | |||||
| MACGrid &out; | |||||
| const Matrix &kernel; | |||||
| }; | |||||
| //! convolves in with 1D kernel (centred at the kernel's midpoint) in the z-direction | |||||
| //! (out must be a grid of zeros) | |||||
| struct apply1DKernelDirZ : public KernelBase { | |||||
| apply1DKernelDirZ(const MACGrid &in, MACGrid &out, const Matrix &kernel) | |||||
| : KernelBase(&in, 0), in(in), out(out), kernel(kernel) | |||||
| { | |||||
| runMessage(); | |||||
| run(); | |||||
| } | |||||
| inline void op(int i, int j, int k, const MACGrid &in, MACGrid &out, const Matrix &kernel) const | |||||
| { | |||||
| int nz = in.getSizeZ(); | |||||
| int kn = kernel.n; | |||||
| int kCentre = kn / 2; | |||||
| for (int m = 0, ind = kn - 1, kk = k - kCentre; m < kn; m++, ind--, kk++) { | |||||
| if (kk < 0) | |||||
| continue; | |||||
| else if (kk >= nz) | |||||
| break; | |||||
| else | |||||
| out(i, j, k) += in(i, j, kk) * kernel(0, ind); | |||||
| } | |||||
| } | |||||
| inline const MACGrid &getArg0() | |||||
| { | |||||
| return in; | |||||
| } | |||||
| typedef MACGrid type0; | |||||
| inline MACGrid &getArg1() | |||||
| { | |||||
| return out; | |||||
| } | |||||
| typedef MACGrid type1; | |||||
| inline const Matrix &getArg2() | |||||
| { | |||||
| return kernel; | |||||
| } | |||||
| typedef Matrix type2; | |||||
| void runMessage() | |||||
| { | |||||
| debMsg("Executing kernel apply1DKernelDirZ ", 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, in, out, kernel); | |||||
| } | |||||
| 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, in, out, kernel); | |||||
| } | |||||
| } | |||||
| 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 MACGrid ∈ | |||||
| MACGrid &out; | |||||
| const Matrix &kernel; | |||||
| }; | |||||
| //! Apply separable Gaussian blur in 2D | |||||
| void applySeparableKernel2D(MACGrid &grid, const FlagGrid &flags, const Matrix &kernel) | |||||
| { | |||||
| // int nx = grid.getSizeX(), ny = grid.getSizeY(); | |||||
| // int kn = kernel.n; | |||||
| // int kCentre = kn / 2; | |||||
| FluidSolver *parent = grid.getParent(); | |||||
| MACGrid orig = MACGrid(parent); | |||||
| orig.copyFrom(grid); | |||||
| MACGrid gridX = MACGrid(parent); | |||||
| apply1DKernelDirX(grid, gridX, kernel); | |||||
| MACGrid gridXY = MACGrid(parent); | |||||
| apply1DKernelDirY(gridX, gridXY, kernel); | |||||
| grid.copyFrom(gridXY); | |||||
| FOR_IJK(grid) | |||||
| { | |||||
| if ((i > 0 && flags.isObstacle(i - 1, j, k)) || (j > 0 && flags.isObstacle(i, j - 1, k)) || | |||||
| flags.isObstacle(i, j, k)) { | |||||
| grid(i, j, k).x = orig(i, j, k).x; | |||||
| grid(i, j, k).y = orig(i, j, k).y; | |||||
| grid(i, j, k).z = orig(i, j, k).z; | |||||
| } | |||||
| } | |||||
| } | |||||
| //! Apply separable Gaussian blur in 3D | |||||
| void applySeparableKernel3D(MACGrid &grid, const FlagGrid &flags, const Matrix &kernel) | |||||
| { | |||||
| // int nx = grid.getSizeX(), ny = grid.getSizeY(), nz = grid.getSizeZ(); | |||||
| // int kn = kernel.n; | |||||
| // int kCentre = kn / 2; | |||||
| FluidSolver *parent = grid.getParent(); | |||||
| MACGrid orig = MACGrid(parent); | |||||
| orig.copyFrom(grid); | |||||
| MACGrid gridX = MACGrid(parent); | |||||
| apply1DKernelDirX(grid, gridX, kernel); | |||||
| MACGrid gridXY = MACGrid(parent); | |||||
| apply1DKernelDirY(gridX, gridXY, kernel); | |||||
| MACGrid gridXYZ = MACGrid(parent); | |||||
| apply1DKernelDirZ(gridXY, gridXYZ, kernel); | |||||
| grid.copyFrom(gridXYZ); | |||||
| FOR_IJK(grid) | |||||
| { | |||||
| if ((i > 0 && flags.isObstacle(i - 1, j, k)) || (j > 0 && flags.isObstacle(i, j - 1, k)) || | |||||
| (k > 0 && flags.isObstacle(i, j, k - 1)) || flags.isObstacle(i, j, k)) { | |||||
| grid(i, j, k).x = orig(i, j, k).x; | |||||
| grid(i, j, k).y = orig(i, j, k).y; | |||||
| grid(i, j, k).z = orig(i, j, k).z; | |||||
| } | |||||
| } | |||||
| } | |||||
| //! Apply separable Gaussian blur in 2D or 3D depending on input dimensions | |||||
| void applySeparableKernel(MACGrid &grid, const FlagGrid &flags, const Matrix &kernel) | |||||
| { | |||||
| if (!grid.is3D()) | |||||
| applySeparableKernel2D(grid, flags, kernel); | |||||
| else | |||||
| applySeparableKernel3D(grid, flags, kernel); | |||||
| } | |||||
| //! Compute r-norm for the stopping criterion | |||||
| Real getRNorm(const MACGrid &x, const MACGrid &z) | |||||
| { | |||||
| MACGrid r = MACGrid(x.getParent()); | |||||
| r.copyFrom(x); | |||||
| r.sub(z); | |||||
| return r.getMaxAbs(); | |||||
| } | |||||
| //! Compute s-norm for the stopping criterion | |||||
| Real getSNorm(const Real rho, const MACGrid &z, const MACGrid &z_prev) | |||||
| { | |||||
| MACGrid s = MACGrid(z_prev.getParent()); | |||||
| s.copyFrom(z_prev); | |||||
| s.sub(z); | |||||
| s.multConst(rho); | |||||
| return s.getMaxAbs(); | |||||
| } | |||||
| //! Compute primal eps for the stopping criterion | |||||
| Real getEpsPri(const Real eps_abs, const Real eps_rel, const MACGrid &x, const MACGrid &z) | |||||
| { | |||||
| Real max_norm = max(x.getMaxAbs(), z.getMaxAbs()); | |||||
| Real eps_pri = sqrt(x.is3D() ? 3.0 : 2.0) * eps_abs + eps_rel * max_norm; | |||||
| return eps_pri; | |||||
| } | |||||
| //! Compute dual eps for the stopping criterion | |||||
| Real getEpsDual(const Real eps_abs, const Real eps_rel, const MACGrid &y) | |||||
| { | |||||
| Real eps_dual = sqrt(y.is3D() ? 3.0 : 2.0) * eps_abs + eps_rel * y.getMaxAbs(); | |||||
| return eps_dual; | |||||
| } | |||||
| //! Create a spiral velocity field in 2D as a test scene (optionally in 3D) | |||||
| void getSpiralVelocity(const FlagGrid &flags, | |||||
| MACGrid &vel, | |||||
| Real strength = 1.0, | |||||
| bool with3D = false) | |||||
| { | |||||
| int nx = flags.getSizeX(), ny = flags.getSizeY(), nz = 1; | |||||
| if (with3D) | |||||
| nz = flags.getSizeZ(); | |||||
| Real midX = 0.5 * (Real)(nx - 1); | |||||
| Real midY = 0.5 * (Real)(ny - 1); | |||||
| Real midZ = 0.5 * (Real)(nz - 1); | |||||
| for (int i = 0; i < nx; i++) { | |||||
| for (int j = 0; j < ny; j++) { | |||||
| for (int k = 0; k < nz; k++) { | |||||
| int idx = flags.index(i, j, k); | |||||
| Real diffX = midX - i; | |||||
| Real diffY = midY - j; | |||||
| Real hypotenuse = sqrt(diffX * diffX + diffY * diffY); | |||||
| if (hypotenuse > 0) { | |||||
| vel[idx].x = diffY / hypotenuse; | |||||
| vel[idx].y = -diffX / hypotenuse; | |||||
| } | |||||
| } | |||||
| } | |||||
| } | |||||
| vel.multConst(strength); | |||||
| } | |||||
| 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, "getSpiralVelocity", !noTiming); | |||||
| PyObject *_retval = 0; | |||||
| { | |||||
| ArgLocker _lock; | |||||
| const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock); | |||||
| MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock); | |||||
| Real strength = _args.getOpt<Real>("strength", 2, 1.0, &_lock); | |||||
| bool with3D = _args.getOpt<bool>("with3D", 3, false, &_lock); | |||||
| _retval = getPyNone(); | |||||
| getSpiralVelocity(flags, vel, strength, with3D); | |||||
| _args.check(); | |||||
| } | |||||
| pbFinalizePlugin(parent, "getSpiralVelocity", !noTiming); | |||||
| return _retval; | |||||
| } | |||||
| catch (std::exception &e) { | |||||
| pbSetError("getSpiralVelocity", e.what()); | |||||
| return 0; | |||||
| } | |||||
| } | |||||
| static const Pb::Register _RP_getSpiralVelocity("", "getSpiralVelocity", _W_0); | |||||
| extern "C" { | |||||
| void PbRegister_getSpiralVelocity() | |||||
| { | |||||
| KEEP_UNUSED(_RP_getSpiralVelocity); | |||||
| } | |||||
| } | |||||
| //! Set the guiding weight W as a gradient in the y-direction | |||||
| void setGradientYWeight( | |||||
| Grid<Real> &W, const int minY, const int maxY, const Real valAtMin, const Real valAtMax) | |||||
| { | |||||
| FOR_IJK(W) | |||||
| { | |||||
| if (minY <= j && j <= maxY) { | |||||
| Real val = valAtMin; | |||||
| if (valAtMax != valAtMin) { | |||||
| Real ratio = (Real)(j - minY) / (Real)(maxY - minY); | |||||
| val = ratio * valAtMax + (1.0 - ratio) * valAtMin; | |||||
| } | |||||
| W(i, j, k) = val; | |||||
| } | |||||
| } | |||||
| } | |||||
| static PyObject *_W_1(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, "setGradientYWeight", !noTiming); | |||||
| PyObject *_retval = 0; | |||||
| { | |||||
| ArgLocker _lock; | |||||
| Grid<Real> &W = *_args.getPtr<Grid<Real>>("W", 0, &_lock); | |||||
| const int minY = _args.get<int>("minY", 1, &_lock); | |||||
| const int maxY = _args.get<int>("maxY", 2, &_lock); | |||||
| const Real valAtMin = _args.get<Real>("valAtMin", 3, &_lock); | |||||
| const Real valAtMax = _args.get<Real>("valAtMax", 4, &_lock); | |||||
| _retval = getPyNone(); | |||||
| setGradientYWeight(W, minY, maxY, valAtMin, valAtMax); | |||||
| _args.check(); | |||||
| } | |||||
| pbFinalizePlugin(parent, "setGradientYWeight", !noTiming); | |||||
| return _retval; | |||||
| } | |||||
| catch (std::exception &e) { | |||||
| pbSetError("setGradientYWeight", e.what()); | |||||
| return 0; | |||||
| } | |||||
| } | |||||
| static const Pb::Register _RP_setGradientYWeight("", "setGradientYWeight", _W_1); | |||||
| extern "C" { | |||||
| void PbRegister_setGradientYWeight() | |||||
| { | |||||
| KEEP_UNUSED(_RP_setGradientYWeight); | |||||
| } | |||||
| } | |||||
| // ***************************************************************************** | |||||
| // More helper functions for fluid guiding | |||||
| //! Apply Gaussian blur (either 2D or 3D) in a separable way | |||||
| void applySeparableGaussianBlur(MACGrid &grid, const FlagGrid &flags, const Matrix &kernel1D) | |||||
| { | |||||
| assertMsg(gBlurPrecomputed, "Error - blue kernel not precomputed"); | |||||
| applySeparableKernel(grid, flags, kernel1D); | |||||
| } | |||||
| //! Precomputation performed before the first PD iteration | |||||
| void ADMM_precompute_Separable(int blurRadius) | |||||
| { | |||||
| if (gBlurPrecomputed) { | |||||
| assertMsg(gBlurKernelRadius == blurRadius, | |||||
| "More than a single blur radius not supported at the moment."); | |||||
| return; | |||||
| } | |||||
| int kernelSize = 2 * blurRadius + 1; | |||||
| gBlurKernel = get1DGaussianBlurKernel(kernelSize, kernelSize); | |||||
| gBlurPrecomputed = true; | |||||
| gBlurKernelRadius = blurRadius; | |||||
| } | |||||
| //! Apply approximate multiplication of inverse(M) | |||||
| void applyApproxInvM(MACGrid &v, const FlagGrid &flags, const MACGrid &invA) | |||||
| { | |||||
| MACGrid v_new = MACGrid(v.getParent()); | |||||
| v_new.copyFrom(v); | |||||
| v_new.mult(invA); | |||||
| applySeparableGaussianBlur(v_new, flags, gBlurKernel); | |||||
| applySeparableGaussianBlur(v_new, flags, gBlurKernel); | |||||
| v_new.multConst(2.0); | |||||
| v_new.mult(invA); | |||||
| v.mult(invA); | |||||
| v.sub(v_new); | |||||
| } | |||||
| //! Precompute Q, a reused quantity in the PD iterations | |||||
| //! Q = 2*G*G*(velT-velC)-sigma*velC | |||||
| void precomputeQ(MACGrid &Q, | |||||
| const FlagGrid &flags, | |||||
| const MACGrid &velT_region, | |||||
| const MACGrid &velC, | |||||
| const Matrix &gBlurKernel, | |||||
| const Real sigma) | |||||
| { | |||||
| Q.copyFrom(velT_region); | |||||
| Q.sub(velC); | |||||
| applySeparableGaussianBlur(Q, flags, gBlurKernel); | |||||
| applySeparableGaussianBlur(Q, flags, gBlurKernel); | |||||
| Q.multConst(2.0); | |||||
| Q.addScaled(velC, -sigma); | |||||
| } | |||||
| //! Precompute inverse(A), a reused quantity in the PD iterations | |||||
| //! A = 2*S^2 + p*I, invA = elementwise 1/A | |||||
| void precomputeInvA(MACGrid &invA, const Grid<Real> &weight, const Real sigma) | |||||
| { | |||||
| FOR_IJK(invA) | |||||
| { | |||||
| Real val = 2 * weight(i, j, k) * weight(i, j, k) + sigma; | |||||
| if (val < 0.01) | |||||
| val = 0.01; | |||||
| Real invVal = 1.0 / val; | |||||
| invA(i, j, k).x = invVal; | |||||
| invA(i, j, k).y = invVal; | |||||
| invA(i, j, k).z = invVal; | |||||
| } | |||||
| } | |||||
| //! proximal operator of f , guiding | |||||
| void prox_f(MACGrid &v, | |||||
| const FlagGrid &flags, | |||||
| const MACGrid &Q, | |||||
| const MACGrid &velC, | |||||
| const Real sigma, | |||||
| const MACGrid &invA) | |||||
| { | |||||
| v.multConst(sigma); | |||||
| v.add(Q); | |||||
| applyApproxInvM(v, flags, invA); | |||||
| v.add(velC); | |||||
| } | |||||
| // ***************************************************************************** | |||||
| // re-uses main pressure solve from pressure.cpp | |||||
| void solvePressure(MACGrid &vel, | |||||
| Grid<Real> &pressure, | |||||
| const FlagGrid &flags, | |||||
| Real cgAccuracy = 1e-3, | |||||
| const Grid<Real> *phi = 0, | |||||
| const Grid<Real> *perCellCorr = 0, | |||||
| const MACGrid *fractions = 0, | |||||
| const MACGrid *obvel = 0, | |||||
| Real gfClamp = 1e-04, | |||||
| Real cgMaxIterFac = 1.5, | |||||
| bool precondition = true, | |||||
| int preconditioner = 1, | |||||
| bool enforceCompatibility = false, | |||||
| bool useL2Norm = false, | |||||
| bool zeroPressureFixing = false, | |||||
| const Grid<Real> *curv = NULL, | |||||
| const Real surfTens = 0.0, | |||||
| Grid<Real> *retRhs = NULL); | |||||
| //! Main function for fluid guiding , includes "regular" pressure solve | |||||
| void PD_fluid_guiding(MACGrid &vel, | |||||
| MACGrid &velT, | |||||
| Grid<Real> &pressure, | |||||
| FlagGrid &flags, | |||||
| Grid<Real> &weight, | |||||
| int blurRadius = 5, | |||||
| Real theta = 1.0, | |||||
| Real tau = 1.0, | |||||
| Real sigma = 1.0, | |||||
| Real epsRel = 1e-3, | |||||
| Real epsAbs = 1e-3, | |||||
| int maxIters = 200, | |||||
| Grid<Real> *phi = 0, | |||||
| Grid<Real> *perCellCorr = 0, | |||||
| MACGrid *fractions = 0, | |||||
| MACGrid *obvel = 0, | |||||
| Real gfClamp = 1e-04, | |||||
| Real cgMaxIterFac = 1.5, | |||||
| Real cgAccuracy = 1e-3, | |||||
| int preconditioner = 1, | |||||
| bool zeroPressureFixing = false, | |||||
| const Grid<Real> *curv = NULL, | |||||
| const Real surfTens = 0.) | |||||
| { | |||||
| FluidSolver *parent = vel.getParent(); | |||||
| // initialize dual/slack variables | |||||
| MACGrid velC = MACGrid(parent); | |||||
| velC.copyFrom(vel); | |||||
| MACGrid x = MACGrid(parent); | |||||
| MACGrid y = MACGrid(parent); | |||||
| MACGrid z = MACGrid(parent); | |||||
| MACGrid x0 = MACGrid(parent); | |||||
| MACGrid z0 = MACGrid(parent); | |||||
| // precomputation | |||||
| ADMM_precompute_Separable(blurRadius); | |||||
| MACGrid Q = MACGrid(parent); | |||||
| precomputeQ(Q, flags, velT, velC, gBlurKernel, sigma); | |||||
| MACGrid invA = MACGrid(parent); | |||||
| precomputeInvA(invA, weight, sigma); | |||||
| // loop | |||||
| int iter = 0; | |||||
| for (iter = 0; iter < maxIters; iter++) { | |||||
| // x-update | |||||
| x0.copyFrom(x); | |||||
| x.multConst(1.0 / sigma); | |||||
| x.add(y); | |||||
| prox_f(x, flags, Q, velC, sigma, invA); | |||||
| x.multConst(-sigma); | |||||
| x.addScaled(y, sigma); | |||||
| x.add(x0); | |||||
| // z-update | |||||
| z0.copyFrom(z); | |||||
| z.addScaled(x, -tau); | |||||
| Real cgAccuracyAdaptive = cgAccuracy; | |||||
| solvePressure(z, | |||||
| pressure, | |||||
| flags, | |||||
| cgAccuracyAdaptive, | |||||
| phi, | |||||
| perCellCorr, | |||||
| fractions, | |||||
| obvel, | |||||
| gfClamp, | |||||
| cgMaxIterFac, | |||||
| true, | |||||
| preconditioner, | |||||
| false, | |||||
| false, | |||||
| zeroPressureFixing, | |||||
| curv, | |||||
| surfTens); | |||||
| // y-update | |||||
| y.copyFrom(z); | |||||
| y.sub(z0); | |||||
| y.multConst(theta); | |||||
| y.add(z); | |||||
| // stopping criterion | |||||
| bool stop = (iter > 0 && getRNorm(z, z0) < getEpsDual(epsAbs, epsRel, z)); | |||||
| if (stop || (iter == maxIters - 1)) | |||||
| break; | |||||
| } | |||||
| // vel_new = z | |||||
| vel.copyFrom(z); | |||||
| debMsg("PD_fluid_guiding iterations:" << iter, 1); | |||||
| } | |||||
| static PyObject *_W_2(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, "PD_fluid_guiding", !noTiming); | |||||
| PyObject *_retval = 0; | |||||
| { | |||||
| ArgLocker _lock; | |||||
| MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock); | |||||
| MACGrid &velT = *_args.getPtr<MACGrid>("velT", 1, &_lock); | |||||
| Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 2, &_lock); | |||||
| FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock); | |||||
| Grid<Real> &weight = *_args.getPtr<Grid<Real>>("weight", 4, &_lock); | |||||
| int blurRadius = _args.getOpt<int>("blurRadius", 5, 5, &_lock); | |||||
| Real theta = _args.getOpt<Real>("theta", 6, 1.0, &_lock); | |||||
| Real tau = _args.getOpt<Real>("tau", 7, 1.0, &_lock); | |||||
| Real sigma = _args.getOpt<Real>("sigma", 8, 1.0, &_lock); | |||||
| Real epsRel = _args.getOpt<Real>("epsRel", 9, 1e-3, &_lock); | |||||
| Real epsAbs = _args.getOpt<Real>("epsAbs", 10, 1e-3, &_lock); | |||||
| int maxIters = _args.getOpt<int>("maxIters", 11, 200, &_lock); | |||||
| Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 12, 0, &_lock); | |||||
| Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 13, 0, &_lock); | |||||
| MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 14, 0, &_lock); | |||||
| MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 15, 0, &_lock); | |||||
| Real gfClamp = _args.getOpt<Real>("gfClamp", 16, 1e-04, &_lock); | |||||
| Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 17, 1.5, &_lock); | |||||
| Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 18, 1e-3, &_lock); | |||||
| int preconditioner = _args.getOpt<int>("preconditioner", 19, 1, &_lock); | |||||
| bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 20, false, &_lock); | |||||
| const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 21, NULL, &_lock); | |||||
| const Real surfTens = _args.getOpt<Real>("surfTens", 22, 0., &_lock); | |||||
| _retval = getPyNone(); | |||||
| PD_fluid_guiding(vel, | |||||
| velT, | |||||
| pressure, | |||||
| flags, | |||||
| weight, | |||||
| blurRadius, | |||||
| theta, | |||||
| tau, | |||||
| sigma, | |||||
| epsRel, | |||||
| epsAbs, | |||||
| maxIters, | |||||
| phi, | |||||
| perCellCorr, | |||||
| fractions, | |||||
| obvel, | |||||
| gfClamp, | |||||
| cgMaxIterFac, | |||||
| cgAccuracy, | |||||
| preconditioner, | |||||
| zeroPressureFixing, | |||||
| curv, | |||||
| surfTens); | |||||
| _args.check(); | |||||
| } | |||||
| pbFinalizePlugin(parent, "PD_fluid_guiding", !noTiming); | |||||
| return _retval; | |||||
| } | |||||
| catch (std::exception &e) { | |||||
| pbSetError("PD_fluid_guiding", e.what()); | |||||
| return 0; | |||||
| } | |||||
| } | |||||
| static const Pb::Register _RP_PD_fluid_guiding("", "PD_fluid_guiding", _W_2); | |||||
| extern "C" { | |||||
| void PbRegister_PD_fluid_guiding() | |||||
| { | |||||
| KEEP_UNUSED(_RP_PD_fluid_guiding); | |||||
| } | |||||
| } | |||||
| //! reset precomputation | |||||
| void releaseBlurPrecomp() | |||||
| { | |||||
| gBlurPrecomputed = false; | |||||
| gBlurKernelRadius = -1; | |||||
| gBlurKernel = 0.f; | |||||
| } | |||||
| static PyObject *_W_3(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, "releaseBlurPrecomp", !noTiming); | |||||
| PyObject *_retval = 0; | |||||
| { | |||||
| ArgLocker _lock; | |||||
| _retval = getPyNone(); | |||||
| releaseBlurPrecomp(); | |||||
| _args.check(); | |||||
| } | |||||
| pbFinalizePlugin(parent, "releaseBlurPrecomp", !noTiming); | |||||
| return _retval; | |||||
| } | |||||
| catch (std::exception &e) { | |||||
| pbSetError("releaseBlurPrecomp", e.what()); | |||||
| return 0; | |||||
| } | |||||
| } | |||||
| static const Pb::Register _RP_releaseBlurPrecomp("", "releaseBlurPrecomp", _W_3); | |||||
| extern "C" { | |||||
| void PbRegister_releaseBlurPrecomp() | |||||
| { | |||||
| KEEP_UNUSED(_RP_releaseBlurPrecomp); | |||||
| } | |||||
| } | |||||
| } // namespace Manta | |||||