Page MenuHome
Differential D3850 Diff 19899 extern/mantaflow/preprocessed/plugin/extforces.cpp

Changeset View

Standalone View

View Options

extern/mantaflow/preprocessed/plugin/extforces.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
* GNU General Public License (GPL)
* http://www.gnu.org/licenses
*
* Set boundary conditions, gravity
*
******************************************************************************/
#include "vectorbase.h"
#include "grid.h"
#include "commonkernels.h"
#include "particle.h"
using namespace std;
namespace Manta {
//! add constant force between fl/fl and fl/em cells
struct KnApplyForceField : public KernelBase {
KnApplyForceField(const FlagGrid &flags,
MACGrid &vel,
const Grid<Vec3> &force,
const Grid<Real> *include,
bool additive,
bool isMAC)
: KernelBase(&flags, 1),
flags(flags),
vel(vel),
force(force),
include(include),
additive(additive),
isMAC(isMAC)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
const FlagGrid &flags,
MACGrid &vel,
const Grid<Vec3> &force,
const Grid<Real> *include,
bool additive,
bool isMAC) const
{
bool curFluid = flags.isFluid(i, j, k);
bool curEmpty = flags.isEmpty(i, j, k);
if (!curFluid && !curEmpty)
return;
if (include && ((*include)(i, j, k) > 0.))
return;
Real forceX = (isMAC) ? force(i, j, k).x : 0.5 * (force(i - 1, j, k).x + force(i, j, k).x);
Real forceY = (isMAC) ? force(i, j, k).y : 0.5 * (force(i, j - 1, k).y + force(i, j, k).y);
Real forceZ = 0.;
if (vel.is3D())
forceZ = (isMAC) ? force(i, j, k).z : 0.5 * (force(i, j, k - 1).z + force(i, j, k).z);
if (flags.isFluid(i - 1, j, k) || (curFluid && flags.isEmpty(i - 1, j, k)))
vel(i, j, k).x = (additive) ? vel(i, j, k).x + forceX : forceX;
if (flags.isFluid(i, j - 1, k) || (curFluid && flags.isEmpty(i, j - 1, k)))
vel(i, j, k).y = (additive) ? vel(i, j, k).y + forceY : forceY;
if (vel.is3D() && (flags.isFluid(i, j, k - 1) || (curFluid && flags.isEmpty(i, j, k - 1))))
vel(i, j, k).z = (additive) ? vel(i, j, k).z + forceZ : forceZ;
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline MACGrid &getArg1()
{
return vel;
}
typedef MACGrid type1;
inline const Grid<Vec3> &getArg2()
{
return force;
}
typedef Grid<Vec3> type2;
inline const Grid<Real> *getArg3()
{
return include;
}
typedef Grid<Real> type3;
inline bool &getArg4()
{
return additive;
}
typedef bool type4;
inline bool &getArg5()
{
return isMAC;
}
typedef bool type5;
void runMessage()
{
debMsg("Executing kernel KnApplyForceField ", 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, vel, force, include, additive, isMAC);
}
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, vel, force, include, additive, isMAC);
}
}
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;
MACGrid &vel;
const Grid<Vec3> &force;
const Grid<Real> *include;
bool additive;
bool isMAC;
};
//! add constant force between fl/fl and fl/em cells
struct KnApplyForce : public KernelBase {
KnApplyForce(
const FlagGrid &flags, MACGrid &vel, Vec3 force, const Grid<Real> *exclude, bool additive)
: KernelBase(&flags, 1),
flags(flags),
vel(vel),
force(force),
exclude(exclude),
additive(additive)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
const FlagGrid &flags,
MACGrid &vel,
Vec3 force,
const Grid<Real> *exclude,
bool additive) const
{
bool curFluid = flags.isFluid(i, j, k);
bool curEmpty = flags.isEmpty(i, j, k);
if (!curFluid && !curEmpty)
return;
if (exclude && ((*exclude)(i, j, k) < 0.))
return;
if (flags.isFluid(i - 1, j, k) || (curFluid && flags.isEmpty(i - 1, j, k)))
vel(i, j, k).x = (additive) ? vel(i, j, k).x + force.x : force.x;
if (flags.isFluid(i, j - 1, k) || (curFluid && flags.isEmpty(i, j - 1, k)))
vel(i, j, k).y = (additive) ? vel(i, j, k).y + force.y : force.y;
if (vel.is3D() && (flags.isFluid(i, j, k - 1) || (curFluid && flags.isEmpty(i, j, k - 1))))
vel(i, j, k).z = (additive) ? vel(i, j, k).z + force.z : force.z;
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline MACGrid &getArg1()
{
return vel;
}
typedef MACGrid type1;
inline Vec3 &getArg2()
{
return force;
}
typedef Vec3 type2;
inline const Grid<Real> *getArg3()
{
return exclude;
}
typedef Grid<Real> type3;
inline bool &getArg4()
{
return additive;
}
typedef bool type4;
void runMessage()
{
debMsg("Executing kernel KnApplyForce ", 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, vel, force, exclude, additive);
}
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, vel, force, exclude, additive);
}
}
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;
MACGrid &vel;
Vec3 force;
const Grid<Real> *exclude;
bool additive;
};
//! add gravity forces to all fluid cells, automatically adapts to different grid sizes
void addGravity(const FlagGrid &flags,
MACGrid &vel,
Vec3 gravity,
const Grid<Real> *exclude = NULL)
{
Vec3 f = gravity * flags.getParent()->getDt() / flags.getDx();
KnApplyForce(flags, vel, f, exclude, true);
}
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, "addGravity", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
Vec3 gravity = _args.get<Vec3>("gravity", 2, &_lock);
const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 3, NULL, &_lock);
_retval = getPyNone();
addGravity(flags, vel, gravity, exclude);
_args.check();
}
pbFinalizePlugin(parent, "addGravity", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("addGravity", e.what());
return 0;
}
}
static const Pb::Register _RP_addGravity("", "addGravity", _W_0);
extern "C" {
void PbRegister_addGravity()
{
KEEP_UNUSED(_RP_addGravity);
}
}
//! add gravity forces to all fluid cells , but dont account for changing cell size
void addGravityNoScale(const FlagGrid &flags,
MACGrid &vel,
const Vec3 &gravity,
const Grid<Real> *exclude = NULL)
{
const Vec3 f = gravity * flags.getParent()->getDt();
KnApplyForce(flags, vel, f, exclude, true);
}
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, "addGravityNoScale", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
const Vec3 &gravity = _args.get<Vec3>("gravity", 2, &_lock);
const Grid<Real> *exclude = _args.getPtrOpt<Grid<Real>>("exclude", 3, NULL, &_lock);
_retval = getPyNone();
addGravityNoScale(flags, vel, gravity, exclude);
_args.check();
}
pbFinalizePlugin(parent, "addGravityNoScale", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("addGravityNoScale", e.what());
return 0;
}
}
static const Pb::Register _RP_addGravityNoScale("", "addGravityNoScale", _W_1);
extern "C" {
void PbRegister_addGravityNoScale()
{
KEEP_UNUSED(_RP_addGravityNoScale);
}
}
//! kernel to add Buoyancy force
struct KnAddBuoyancy : public KernelBase {
KnAddBuoyancy(const FlagGrid &flags, const Grid<Real> &factor, MACGrid &vel, Vec3 strength)
: KernelBase(&flags, 1), flags(flags), factor(factor), vel(vel), strength(strength)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
const FlagGrid &flags,
const Grid<Real> &factor,
MACGrid &vel,
Vec3 strength) const
{
if (!flags.isFluid(i, j, k))
return;
if (flags.isFluid(i - 1, j, k))
vel(i, j, k).x += (0.5 * strength.x) * (factor(i, j, k) + factor(i - 1, j, k));
if (flags.isFluid(i, j - 1, k))
vel(i, j, k).y += (0.5 * strength.y) * (factor(i, j, k) + factor(i, j - 1, k));
if (vel.is3D() && flags.isFluid(i, j, k - 1))
vel(i, j, k).z += (0.5 * strength.z) * (factor(i, j, k) + factor(i, j, k - 1));
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline const Grid<Real> &getArg1()
{
return factor;
}
typedef Grid<Real> type1;
inline MACGrid &getArg2()
{
return vel;
}
typedef MACGrid type2;
inline Vec3 &getArg3()
{
return strength;
}
typedef Vec3 type3;
void runMessage()
{
debMsg("Executing kernel KnAddBuoyancy ", 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, factor, vel, strength);
}
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, factor, vel, strength);
}
}
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;
const Grid<Real> &factor;
MACGrid &vel;
Vec3 strength;
};
//! add Buoyancy force based on fctor (e.g. smoke density)
void addBuoyancy(const FlagGrid &flags,
const Grid<Real> &density,
MACGrid &vel,
Vec3 gravity,
Real coefficient = 1.)
{
Vec3 f = -gravity * flags.getParent()->getDt() / flags.getParent()->getDx() * coefficient;
KnAddBuoyancy(flags, density, vel, f);
}
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, "addBuoyancy", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
const Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 1, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 2, &_lock);
Vec3 gravity = _args.get<Vec3>("gravity", 3, &_lock);
Real coefficient = _args.getOpt<Real>("coefficient", 4, 1., &_lock);
_retval = getPyNone();
addBuoyancy(flags, density, vel, gravity, coefficient);
_args.check();
}
pbFinalizePlugin(parent, "addBuoyancy", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("addBuoyancy", e.what());
return 0;
}
}
static const Pb::Register _RP_addBuoyancy("", "addBuoyancy", _W_2);
extern "C" {
void PbRegister_addBuoyancy()
{
KEEP_UNUSED(_RP_addBuoyancy);
}
}
// inflow / outflow boundaries
//! helper to parse openbounds string [xXyYzZ] , convert to vec3
inline void convertDescToVec(const string &desc, Vector3D<bool> &lo, Vector3D<bool> &up)
{
for (size_t i = 0; i < desc.size(); i++) {
if (desc[i] == 'x')
lo.x = true;
else if (desc[i] == 'y')
lo.y = true;
else if (desc[i] == 'z')
lo.z = true;
else if (desc[i] == 'X')
up.x = true;
else if (desc[i] == 'Y')
up.y = true;
else if (desc[i] == 'Z')
up.z = true;
else
errMsg("invalid character in boundary description string. Only [xyzXYZ] allowed.");
}
}
//! add empty and outflow flag to cells of open boundaries
void setOpenBound(FlagGrid &flags,
int bWidth,
string openBound = "",
int type = FlagGrid::TypeOutflow | FlagGrid::TypeEmpty)
{
if (openBound == "")
return;
Vector3D<bool> lo, up;
convertDescToVec(openBound, lo, up);
FOR_IJK(flags)
{
bool loX = lo.x && i <= bWidth; // a cell which belongs to the lower x open bound
bool loY = lo.y && j <= bWidth;
bool upX = up.x && i >= flags.getSizeX() - bWidth -
1; // a cell which belongs to the upper x open bound
bool upY = up.y && j >= flags.getSizeY() - bWidth - 1;
bool innerI = i > bWidth &&
i < flags.getSizeX() - bWidth -
1; // a cell which does not belong to the lower or upper x bound
bool innerJ = j > bWidth && j < flags.getSizeY() - bWidth - 1;
// when setting boundaries to open: don't set shared part of wall to empty if neighboring wall
// is not open
if ((!flags.is3D()) && (loX || upX || loY || upY)) {
if ((loX || upX || innerI) && (loY || upY || innerJ) && flags.isObstacle(i, j, k))
flags(i, j, k) = type;
}
else {
bool loZ = lo.z && k <= bWidth; // a cell which belongs to the lower z open bound
bool upZ = up.z && k >= flags.getSizeZ() - bWidth -
1; // a cell which belongs to the upper z open bound
bool innerK = k > bWidth &&
k < flags.getSizeZ() - bWidth -
1; // a cell which does not belong to the lower or upper z bound
if (loX || upX || loY || upY || loZ || upZ) {
if ((loX || upX || innerI) && (loY || upY || innerJ) && (loZ || upZ || innerK) &&
flags.isObstacle(i, j, k))
flags(i, j, k) = type;
}
}
}
}
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, "setOpenBound", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
int bWidth = _args.get<int>("bWidth", 1, &_lock);
string openBound = _args.getOpt<string>("openBound", 2, "", &_lock);
int type = _args.getOpt<int>("type", 3, FlagGrid::TypeOutflow | FlagGrid::TypeEmpty, &_lock);
_retval = getPyNone();
setOpenBound(flags, bWidth, openBound, type);
_args.check();
}
pbFinalizePlugin(parent, "setOpenBound", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("setOpenBound", e.what());
return 0;
}
}
static const Pb::Register _RP_setOpenBound("", "setOpenBound", _W_3);
extern "C" {
void PbRegister_setOpenBound()
{
KEEP_UNUSED(_RP_setOpenBound);
}
}
//! delete fluid and ensure empty flag in outflow cells, delete particles and density and set phi
//! to 0.5
void resetOutflow(FlagGrid &flags,
Grid<Real> *phi = 0,
BasicParticleSystem *parts = 0,
Grid<Real> *real = 0,
Grid<int> *index = 0,
ParticleIndexSystem *indexSys = 0)
{
// check if phi and parts -> pindex and gpi already created -> access particles from cell index,
// avoid extra looping over particles
if (parts && (!index || !indexSys)) {
if (phi)
debMsg(
"resetOpenBound for phi and particles, but missing index and indexSys for enhanced "
"particle access!",
1);
for (int idx = 0; idx < (int)parts->size(); idx++)
if (parts->isActive(idx) && flags.isInBounds(parts->getPos(idx)) &&
flags.isOutflow(parts->getPos(idx)))
parts->kill(idx);
}
FOR_IJK(flags)
{
if (flags.isOutflow(i, j, k)) {
flags(i, j, k) = (flags(i, j, k) | FlagGrid::TypeEmpty) &
~FlagGrid::TypeFluid; // make sure there is not fluid flag set and to reset
// the empty flag
// the particles in a cell i,j,k are particles[index(i,j,k)] to particles[index(i+1,j,k)-1]
if (parts && index && indexSys) {
int isysIdxS = index->index(i, j, k);
int pStart = (*index)(isysIdxS), pEnd = 0;
if (flags.isInBounds(isysIdxS + 1))
pEnd = (*index)(isysIdxS + 1);
else
pEnd = indexSys->size();
// now loop over particles in cell
for (int p = pStart; p < pEnd; ++p) {
int psrc = (*indexSys)[p].sourceIndex;
if (parts->isActive(psrc) && flags.isInBounds(parts->getPos(psrc)))
parts->kill(psrc);
}
}
if (phi)
(*phi)(i, j, k) = 0.5;
if (real)
(*real)(i, j, k) = 0;
}
}
if (parts)
parts->doCompress();
}
static PyObject *_W_4(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, "resetOutflow", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 1, 0, &_lock);
BasicParticleSystem *parts = _args.getPtrOpt<BasicParticleSystem>("parts", 2, 0, &_lock);
Grid<Real> *real = _args.getPtrOpt<Grid<Real>>("real", 3, 0, &_lock);
Grid<int> *index = _args.getPtrOpt<Grid<int>>("index", 4, 0, &_lock);
ParticleIndexSystem *indexSys = _args.getPtrOpt<ParticleIndexSystem>(
"indexSys", 5, 0, &_lock);
_retval = getPyNone();
resetOutflow(flags, phi, parts, real, index, indexSys);
_args.check();
}
pbFinalizePlugin(parent, "resetOutflow", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("resetOutflow", e.what());
return 0;
}
}
static const Pb::Register _RP_resetOutflow("", "resetOutflow", _W_4);
extern "C" {
void PbRegister_resetOutflow()
{
KEEP_UNUSED(_RP_resetOutflow);
}
}
//! enforce a constant inflow/outflow at the grid boundaries
struct KnSetInflow : public KernelBase {
KnSetInflow(MACGrid &vel, int dim, int p0, const Vec3 &val)
: KernelBase(&vel, 0), vel(vel), dim(dim), p0(p0), val(val)
{
runMessage();
run();
}
inline void op(int i, int j, int k, MACGrid &vel, int dim, int p0, const Vec3 &val) const
{
Vec3i p(i, j, k);
if (p[dim] == p0 || p[dim] == p0 + 1)
vel(i, j, k) = val;
}
inline MACGrid &getArg0()
{
return vel;
}
typedef MACGrid type0;
inline int &getArg1()
{
return dim;
}
typedef int type1;
inline int &getArg2()
{
return p0;
}
typedef int type2;
inline const Vec3 &getArg3()
{
return val;
}
typedef Vec3 type3;
void runMessage()
{
debMsg("Executing kernel KnSetInflow ", 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, vel, dim, p0, val);
}
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, vel, dim, p0, val);
}
}
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 &vel;
int dim;
int p0;
const Vec3 &val;
};
//! enforce a constant inflow/outflow at the grid boundaries
void setInflowBcs(MACGrid &vel, string dir, Vec3 value)
{
for (size_t i = 0; i < dir.size(); i++) {
if (dir[i] >= 'x' && dir[i] <= 'z') {
int dim = dir[i] - 'x';
KnSetInflow(vel, dim, 0, value);
}
else if (dir[i] >= 'X' && dir[i] <= 'Z') {
int dim = dir[i] - 'X';
KnSetInflow(vel, dim, vel.getSize()[dim] - 1, value);
}
else
errMsg("invalid character in direction string. Only [xyzXYZ] allowed.");
}
}
static PyObject *_W_5(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, "setInflowBcs", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
string dir = _args.get<string>("dir", 1, &_lock);
Vec3 value = _args.get<Vec3>("value", 2, &_lock);
_retval = getPyNone();
setInflowBcs(vel, dir, value);
_args.check();
}
pbFinalizePlugin(parent, "setInflowBcs", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("setInflowBcs", e.what());
return 0;
}
}
static const Pb::Register _RP_setInflowBcs("", "setInflowBcs", _W_5);
extern "C" {
void PbRegister_setInflowBcs()
{
KEEP_UNUSED(_RP_setInflowBcs);
}
}
// set obstacle boundary conditions
//! set no-stick wall boundary condition between ob/fl and ob/ob cells
struct KnSetWallBcs : public KernelBase {
KnSetWallBcs(const FlagGrid &flags, MACGrid &vel, const MACGrid *obvel)
: KernelBase(&flags, 0), flags(flags), vel(vel), obvel(obvel)
{
runMessage();
run();
}
inline void op(
int i, int j, int k, const FlagGrid &flags, MACGrid &vel, const MACGrid *obvel) const
{
bool curFluid = flags.isFluid(i, j, k);
bool curObs = flags.isObstacle(i, j, k);
Vec3 bcsVel(0., 0., 0.);
if (!curFluid && !curObs)
return;
if (obvel) {
bcsVel.x = (*obvel)(i, j, k).x;
bcsVel.y = (*obvel)(i, j, k).y;
if ((*obvel).is3D())
bcsVel.z = (*obvel)(i, j, k).z;
}
// we use i>0 instead of bnd=1 to check outer wall
if (i > 0 && flags.isObstacle(i - 1, j, k))
vel(i, j, k).x = bcsVel.x;
if (i > 0 && curObs && flags.isFluid(i - 1, j, k))
vel(i, j, k).x = bcsVel.x;
if (j > 0 && flags.isObstacle(i, j - 1, k))
vel(i, j, k).y = bcsVel.y;
if (j > 0 && curObs && flags.isFluid(i, j - 1, k))
vel(i, j, k).y = bcsVel.y;
if (!vel.is3D()) {
vel(i, j, k).z = 0;
}
else {
if (k > 0 && flags.isObstacle(i, j, k - 1))
vel(i, j, k).z = bcsVel.z;
if (k > 0 && curObs && flags.isFluid(i, j, k - 1))
vel(i, j, k).z = bcsVel.z;
}
if (curFluid) {
if ((i > 0 && flags.isStick(i - 1, j, k)) ||
(i < flags.getSizeX() - 1 && flags.isStick(i + 1, j, k)))
vel(i, j, k).y = vel(i, j, k).z = 0;
if ((j > 0 && flags.isStick(i, j - 1, k)) ||
(j < flags.getSizeY() - 1 && flags.isStick(i, j + 1, k)))
vel(i, j, k).x = vel(i, j, k).z = 0;
if (vel.is3D() && ((k > 0 && flags.isStick(i, j, k - 1)) ||
(k < flags.getSizeZ() - 1 && flags.isStick(i, j, k + 1))))
vel(i, j, k).x = vel(i, j, k).y = 0;
}
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline MACGrid &getArg1()
{
return vel;
}
typedef MACGrid type1;
inline const MACGrid *getArg2()
{
return obvel;
}
typedef MACGrid type2;
void runMessage()
{
debMsg("Executing kernel KnSetWallBcs ", 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, vel, obvel);
}
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, vel, obvel);
}
}
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;
MACGrid &vel;
const MACGrid *obvel;
};
//! set wall BCs for fill fraction mode, note - only needs obstacle SDF
struct KnSetWallBcsFrac : public KernelBase {
KnSetWallBcsFrac(const FlagGrid &flags,
const MACGrid &vel,
MACGrid &velTarget,
const MACGrid *obvel,
const Grid<Real> *phiObs,
const int &boundaryWidth = 0)
: KernelBase(&flags, 0),
flags(flags),
vel(vel),
velTarget(velTarget),
obvel(obvel),
phiObs(phiObs),
boundaryWidth(boundaryWidth)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
const FlagGrid &flags,
const MACGrid &vel,
MACGrid &velTarget,
const MACGrid *obvel,
const Grid<Real> *phiObs,
const int &boundaryWidth = 0) const
{
bool curFluid = flags.isFluid(i, j, k);
bool curObs = flags.isObstacle(i, j, k);
velTarget(i, j, k) = vel(i, j, k);
if (!curFluid && !curObs)
return;
// zero normal component in all obstacle regions
if (flags.isInBounds(Vec3i(i, j, k), 1)) {
if (curObs | flags.isObstacle(i - 1, j, k)) {
Vec3 dphi(0., 0., 0.);
const Real tmp1 = (phiObs->get(i, j, k) + phiObs->get(i - 1, j, k)) * .5;
Real tmp2 = (phiObs->get(i, j + 1, k) + phiObs->get(i - 1, j + 1, k)) * .5;
Real phi1 = (tmp1 + tmp2) * .5;
tmp2 = (phiObs->get(i, j - 1, k) + phiObs->get(i - 1, j - 1, k)) * .5;
Real phi2 = (tmp1 + tmp2) * .5;
dphi.x = phiObs->get(i, j, k) - phiObs->get(i - 1, j, k);
dphi.y = phi1 - phi2;
if (phiObs->is3D()) {
tmp2 = (phiObs->get(i, j, k + 1) + phiObs->get(i - 1, j, k + 1)) * .5;
phi1 = (tmp1 + tmp2) * .5;
tmp2 = (phiObs->get(i, j, k - 1) + phiObs->get(i - 1, j, k - 1)) * .5;
phi2 = (tmp1 + tmp2) * .5;
dphi.z = phi1 - phi2;
}
normalize(dphi);
Vec3 velMAC = vel.getAtMACX(i, j, k);
velTarget(i, j, k).x = velMAC.x - dot(dphi, velMAC) * dphi.x;
}
if (curObs | flags.isObstacle(i, j - 1, k)) {
Vec3 dphi(0., 0., 0.);
const Real tmp1 = (phiObs->get(i, j, k) + phiObs->get(i, j - 1, k)) * .5;
Real tmp2 = (phiObs->get(i + 1, j, k) + phiObs->get(i + 1, j - 1, k)) * .5;
Real phi1 = (tmp1 + tmp2) * .5;
tmp2 = (phiObs->get(i - 1, j, k) + phiObs->get(i - 1, j - 1, k)) * .5;
Real phi2 = (tmp1 + tmp2) * .5;
dphi.x = phi1 - phi2;
dphi.y = phiObs->get(i, j, k) - phiObs->get(i, j - 1, k);
if (phiObs->is3D()) {
tmp2 = (phiObs->get(i, j, k + 1) + phiObs->get(i, j - 1, k + 1)) * .5;
phi1 = (tmp1 + tmp2) * .5;
tmp2 = (phiObs->get(i, j, k - 1) + phiObs->get(i, j - 1, k - 1)) * .5;
phi2 = (tmp1 + tmp2) * .5;
dphi.z = phi1 - phi2;
}
normalize(dphi);
Vec3 velMAC = vel.getAtMACY(i, j, k);
velTarget(i, j, k).y = velMAC.y - dot(dphi, velMAC) * dphi.y;
}
if (phiObs->is3D() && (curObs | flags.isObstacle(i, j, k - 1))) {
Vec3 dphi(0., 0., 0.);
const Real tmp1 = (phiObs->get(i, j, k) + phiObs->get(i, j, k - 1)) * .5;
Real tmp2;
tmp2 = (phiObs->get(i + 1, j, k) + phiObs->get(i + 1, j, k - 1)) * .5;
Real phi1 = (tmp1 + tmp2) * .5;
tmp2 = (phiObs->get(i - 1, j, k) + phiObs->get(i - 1, j, k - 1)) * .5;
Real phi2 = (tmp1 + tmp2) * .5;
dphi.x = phi1 - phi2;
tmp2 = (phiObs->get(i, j + 1, k) + phiObs->get(i, j + 1, k - 1)) * .5;
phi1 = (tmp1 + tmp2) * .5;
tmp2 = (phiObs->get(i, j - 1, k) + phiObs->get(i, j - 1, k - 1)) * .5;
phi2 = (tmp1 + tmp2) * .5;
dphi.y = phi1 - phi2;
dphi.z = phiObs->get(i, j, k) - phiObs->get(i, j, k - 1);
normalize(dphi);
Vec3 velMAC = vel.getAtMACZ(i, j, k);
velTarget(i, j, k).z = velMAC.z - dot(dphi, velMAC) * dphi.z;
}
} // not at boundary
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline const MACGrid &getArg1()
{
return vel;
}
typedef MACGrid type1;
inline MACGrid &getArg2()
{
return velTarget;
}
typedef MACGrid type2;
inline const MACGrid *getArg3()
{
return obvel;
}
typedef MACGrid type3;
inline const Grid<Real> *getArg4()
{
return phiObs;
}
typedef Grid<Real> type4;
inline const int &getArg5()
{
return boundaryWidth;
}
typedef int type5;
void runMessage()
{
debMsg("Executing kernel KnSetWallBcsFrac ", 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, vel, velTarget, obvel, phiObs, boundaryWidth);
}
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, vel, velTarget, obvel, phiObs, boundaryWidth);
}
}
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;
const MACGrid &vel;
MACGrid &velTarget;
const MACGrid *obvel;
const Grid<Real> *phiObs;
const int &boundaryWidth;
};
//! set zero normal velocity boundary condition on walls
// (optionally with second order accuracy using the obstacle SDF , fractions grid currentlyl not
// needed)
void setWallBcs(const FlagGrid &flags,
MACGrid &vel,
const MACGrid *obvel = 0,
const MACGrid *fractions = 0,
const Grid<Real> *phiObs = 0,
int boundaryWidth = 0)
{
if (!phiObs || !fractions) {
KnSetWallBcs(flags, vel, obvel);
}
else {
MACGrid tmpvel(vel.getParent());
KnSetWallBcsFrac(flags, vel, tmpvel, obvel, phiObs, boundaryWidth);
vel.swap(tmpvel);
}
}
static PyObject *_W_6(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, "setWallBcs", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 2, 0, &_lock);
const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 3, 0, &_lock);
const Grid<Real> *phiObs = _args.getPtrOpt<Grid<Real>>("phiObs", 4, 0, &_lock);
int boundaryWidth = _args.getOpt<int>("boundaryWidth", 5, 0, &_lock);
_retval = getPyNone();
setWallBcs(flags, vel, obvel, fractions, phiObs, boundaryWidth);
_args.check();
}
pbFinalizePlugin(parent, "setWallBcs", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("setWallBcs", e.what());
return 0;
}
}
static const Pb::Register _RP_setWallBcs("", "setWallBcs", _W_6);
extern "C" {
void PbRegister_setWallBcs()
{
KEEP_UNUSED(_RP_setWallBcs);
}
}
//! add Forces between fl/fl and fl/em cells (interpolate cell centered forces to MAC grid)
struct KnAddForceIfLower : public KernelBase {
KnAddForceIfLower(const FlagGrid &flags, MACGrid &vel, const Grid<Vec3> &force)
: KernelBase(&flags, 1), flags(flags), vel(vel), force(force)
{
runMessage();
run();
}
inline void op(
int i, int j, int k, const FlagGrid &flags, MACGrid &vel, const Grid<Vec3> &force) const
{
bool curFluid = flags.isFluid(i, j, k);
bool curEmpty = flags.isEmpty(i, j, k);
if (!curFluid && !curEmpty)
return;
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);
}
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);
}
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);
}
}
inline const FlagGrid &getArg0()
{
return flags;
}
typedef FlagGrid type0;
inline MACGrid &getArg1()
{
return vel;
}
typedef MACGrid type1;
inline const Grid<Vec3> &getArg2()
{
return force;
}
typedef Grid<Vec3> type2;
void runMessage()
{
debMsg("Executing kernel KnAddForceIfLower ", 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, vel, force);
}
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, vel, force);
}
}
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;
MACGrid &vel;
const Grid<Vec3> &force;
};
// Initial velocity for smoke
void setInitialVelocity(const FlagGrid &flags, MACGrid &vel, const Grid<Vec3> &invel)
{
KnAddForceIfLower(flags, vel, invel);
}
static PyObject *_W_7(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, "setInitialVelocity", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
const Grid<Vec3> &invel = *_args.getPtr<Grid<Vec3>>("invel", 2, &_lock);
_retval = getPyNone();
setInitialVelocity(flags, vel, invel);
_args.check();
}
pbFinalizePlugin(parent, "setInitialVelocity", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("setInitialVelocity", e.what());
return 0;
}
}
static const Pb::Register _RP_setInitialVelocity("", "setInitialVelocity", _W_7);
extern "C" {
void PbRegister_setInitialVelocity()
{
KEEP_UNUSED(_RP_setInitialVelocity);
}
}
//! Kernel: gradient norm operator
struct KnConfForce : public KernelBase {
KnConfForce(Grid<Vec3> &force,
const Grid<Real> &grid,
const Grid<Vec3> &curl,
Real str,
const Grid<Real> *strGrid)
: KernelBase(&force, 1), force(force), grid(grid), curl(curl), str(str), strGrid(strGrid)
{
runMessage();
run();
}
inline void op(int i,
int j,
int k,
Grid<Vec3> &force,
const Grid<Real> &grid,
const Grid<Vec3> &curl,
Real str,
const Grid<Real> *strGrid) const
{
Vec3 grad = 0.5 * Vec3(grid(i + 1, j, k) - grid(i - 1, j, k),
grid(i, j + 1, k) - grid(i, j - 1, k),
0.);
if (grid.is3D())
grad[2] = 0.5 * (grid(i, j, k + 1) - grid(i, j, k - 1));
normalize(grad);
if (strGrid)
str += (*strGrid)(i, j, k);
force(i, j, k) = str * cross(grad, curl(i, j, k));
}
inline Grid<Vec3> &getArg0()
{
return force;
}
typedef Grid<Vec3> type0;
inline const Grid<Real> &getArg1()
{
return grid;
}
typedef Grid<Real> type1;
inline const Grid<Vec3> &getArg2()
{
return curl;
}
typedef Grid<Vec3> type2;
inline Real &getArg3()
{
return str;
}
typedef Real type3;
inline const Grid<Real> *getArg4()
{
return strGrid;
}
typedef Grid<Real> type4;
void runMessage()
{
debMsg("Executing kernel KnConfForce ", 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, force, grid, curl, str, strGrid);
}
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, force, grid, curl, str, strGrid);
}
}
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);
}
Grid<Vec3> &force;
const Grid<Real> &grid;
const Grid<Vec3> &curl;
Real str;
const Grid<Real> *strGrid;
};
void vorticityConfinement(MACGrid &vel,
const FlagGrid &flags,
Real strengthGlobal = 0,
const Grid<Real> *strengthCell = NULL)
{
Grid<Vec3> velCenter(flags.getParent()), curl(flags.getParent()), force(flags.getParent());
Grid<Real> norm(flags.getParent());
GetCentered(velCenter, vel);
CurlOp(velCenter, curl);
GridNorm(norm, curl);
KnConfForce(force, norm, curl, strengthGlobal, strengthCell);
KnApplyForceField(flags, vel, force, NULL, true, false);
}
static PyObject *_W_8(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, "vorticityConfinement", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock);
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 1, &_lock);
Real strengthGlobal = _args.getOpt<Real>("strengthGlobal", 2, 0, &_lock);
const Grid<Real> *strengthCell = _args.getPtrOpt<Grid<Real>>(
"strengthCell", 3, NULL, &_lock);
_retval = getPyNone();
vorticityConfinement(vel, flags, strengthGlobal, strengthCell);
_args.check();
}
pbFinalizePlugin(parent, "vorticityConfinement", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("vorticityConfinement", e.what());
return 0;
}
}
static const Pb::Register _RP_vorticityConfinement("", "vorticityConfinement", _W_8);
extern "C" {
void PbRegister_vorticityConfinement()
{
KEEP_UNUSED(_RP_vorticityConfinement);
}
}
void addForceField(const FlagGrid &flags,
MACGrid &vel,
const Grid<Vec3> &force,
const Grid<Real> *region = NULL,
bool isMAC = false)
{
KnApplyForceField(flags, vel, force, region, true, isMAC);
}
static PyObject *_W_9(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, "addForceField", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
const Grid<Vec3> &force = *_args.getPtr<Grid<Vec3>>("force", 2, &_lock);
const Grid<Real> *region = _args.getPtrOpt<Grid<Real>>("region", 3, NULL, &_lock);
bool isMAC = _args.getOpt<bool>("isMAC", 4, false, &_lock);
_retval = getPyNone();
addForceField(flags, vel, force, region, isMAC);
_args.check();
}
pbFinalizePlugin(parent, "addForceField", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("addForceField", e.what());
return 0;
}
}
static const Pb::Register _RP_addForceField("", "addForceField", _W_9);
extern "C" {
void PbRegister_addForceField()
{
KEEP_UNUSED(_RP_addForceField);
}
}
void setForceField(const FlagGrid &flags,
MACGrid &vel,
const Grid<Vec3> &force,
const Grid<Real> *region = NULL,
bool isMAC = false)
{
KnApplyForceField(flags, vel, force, region, false, isMAC);
}
static PyObject *_W_10(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, "setForceField", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
const Grid<Vec3> &force = *_args.getPtr<Grid<Vec3>>("force", 2, &_lock);
const Grid<Real> *region = _args.getPtrOpt<Grid<Real>>("region", 3, NULL, &_lock);
bool isMAC = _args.getOpt<bool>("isMAC", 4, false, &_lock);
_retval = getPyNone();
setForceField(flags, vel, force, region, isMAC);
_args.check();
}
pbFinalizePlugin(parent, "setForceField", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("setForceField", e.what());
return 0;
}
}
static const Pb::Register _RP_setForceField("", "setForceField", _W_10);
extern "C" {
void PbRegister_setForceField()
{
KEEP_UNUSED(_RP_setForceField);
}
}
void dissolveSmoke(const FlagGrid &flags,
Grid<Real> &density,
Grid<Real> *heat = NULL,
Grid<Real> *red = NULL,
Grid<Real> *green = NULL,
Grid<Real> *blue = NULL,
int speed = 5,
bool logFalloff = true)
{
float dydx = 1.0f / (float)speed; // max density/speed = dydx
float fac = 1.0f - dydx;
FOR_IJK_BND(density, 0)
{
bool curFluid = flags.isFluid(i, j, k);
if (!curFluid)
continue;
if (logFalloff) {
density(i, j, k) *= fac;
if (heat) {
(*heat)(i, j, k) *= fac;
}
if (red) {
(*red)(i, j, k) *= fac;
(*green)(i, j, k) *= fac;
(*blue)(i, j, k) *= fac;
}
}
else { // linear falloff
float d = density(i, j, k);
density(i, j, k) -= dydx;
if (density(i, j, k) < 0.0f)
density(i, j, k) = 0.0f;
if (heat) {
if (fabs((*heat)(i, j, k)) < dydx)
(*heat)(i, j, k) = 0.0f;
else if ((*heat)(i, j, k) > 0.0f)
(*heat)(i, j, k) -= dydx;
else if ((*heat)(i, j, k) < 0.0f)
(*heat)(i, j, k) += dydx;
}
if (red && notZero(d)) {
(*red)(i, j, k) *= (density(i, j, k) / d);
(*green)(i, j, k) *= (density(i, j, k) / d);
(*blue)(i, j, k) *= (density(i, j, k) / d);
}
}
}
}
static PyObject *_W_11(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, "dissolveSmoke", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 1, &_lock);
Grid<Real> *heat = _args.getPtrOpt<Grid<Real>>("heat", 2, NULL, &_lock);
Grid<Real> *red = _args.getPtrOpt<Grid<Real>>("red", 3, NULL, &_lock);
Grid<Real> *green = _args.getPtrOpt<Grid<Real>>("green", 4, NULL, &_lock);
Grid<Real> *blue = _args.getPtrOpt<Grid<Real>>("blue", 5, NULL, &_lock);
int speed = _args.getOpt<int>("speed", 6, 5, &_lock);
bool logFalloff = _args.getOpt<bool>("logFalloff", 7, true, &_lock);
_retval = getPyNone();
dissolveSmoke(flags, density, heat, red, green, blue, speed, logFalloff);
_args.check();
}
pbFinalizePlugin(parent, "dissolveSmoke", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("dissolveSmoke", e.what());
return 0;
}
}
static const Pb::Register _RP_dissolveSmoke("", "dissolveSmoke", _W_11);
extern "C" {
void PbRegister_dissolveSmoke()
{
KEEP_UNUSED(_RP_dissolveSmoke);
}
}
} // namespace Manta