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

Changeset View

Standalone View

View Options

extern/mantaflow/preprocessed/plugin/apic.cpp

  • This file was added.
// DO NOT EDIT !
// This file is generated using the MantaFlow preprocessor (prep generate).
// ----------------------------------------------------------------------------
//
// MantaFlow fluid solver framework
// Copyright 2016-2017 Kiwon Um, 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
//
// Affine Particle-In-Cell
//
// ----------------------------------------------------------------------------
#include "particle.h"
#include "grid.h"
namespace Manta {
struct knApicMapLinearVec3ToMACGrid : public KernelBase {
knApicMapLinearVec3ToMACGrid(const BasicParticleSystem &p,
MACGrid &mg,
MACGrid &vg,
const ParticleDataImpl<Vec3> &vp,
const ParticleDataImpl<Vec3> &cpx,
const ParticleDataImpl<Vec3> &cpy,
const ParticleDataImpl<Vec3> &cpz,
const ParticleDataImpl<int> *ptype,
const int exclude)
: KernelBase(p.size()),
p(p),
mg(mg),
vg(vg),
vp(vp),
cpx(cpx),
cpy(cpy),
cpz(cpz),
ptype(ptype),
exclude(exclude)
{
runMessage();
run();
}
inline void op(IndexInt idx,
const BasicParticleSystem &p,
MACGrid &mg,
MACGrid &vg,
const ParticleDataImpl<Vec3> &vp,
const ParticleDataImpl<Vec3> &cpx,
const ParticleDataImpl<Vec3> &cpy,
const ParticleDataImpl<Vec3> &cpz,
const ParticleDataImpl<int> *ptype,
const int exclude)
{
if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
return;
const IndexInt dX[2] = {0, vg.getStrideX()};
const IndexInt dY[2] = {0, vg.getStrideY()};
const IndexInt dZ[2] = {0, vg.getStrideZ()};
const Vec3 &pos = p[idx].pos, &vel = vp[idx];
const IndexInt fi = static_cast<IndexInt>(pos.x), fj = static_cast<IndexInt>(pos.y),
fk = static_cast<IndexInt>(pos.z);
const IndexInt ci = static_cast<IndexInt>(pos.x - 0.5),
cj = static_cast<IndexInt>(pos.y - 0.5),
ck = static_cast<IndexInt>(pos.z - 0.5);
const Real wfi = clamp(pos.x - fi, Real(0), Real(1)),
wfj = clamp(pos.y - fj, Real(0), Real(1)),
wfk = clamp(pos.z - fk, Real(0), Real(1));
const Real wci = clamp(Real(pos.x - ci - 0.5), Real(0), Real(1)),
wcj = clamp(Real(pos.y - cj - 0.5), Real(0), Real(1)),
wck = clamp(Real(pos.z - ck - 0.5), Real(0), Real(1));
// TODO: check index for safety
{ // u-face
const IndexInt gidx = fi * dX[1] + cj * dY[1] + ck * dZ[1];
const Vec3 gpos(fi, cj + 0.5, ck + 0.5);
const Real wi[2] = {Real(1) - wfi, wfi};
const Real wj[2] = {Real(1) - wcj, wcj};
const Real wk[2] = {Real(1) - wck, wck};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k) {
const Real w = wi[i] * wj[j] * wk[k];
mg[gidx + dX[i] + dY[j] + dZ[k]].x += w;
vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * vel.x;
vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * dot(cpx[idx], gpos + Vec3(i, j, k) - pos);
}
}
{ // v-face
const IndexInt gidx = ci * dX[1] + fj * dY[1] + ck * dZ[1];
const Vec3 gpos(ci + 0.5, fj, ck + 0.5);
const Real wi[2] = {Real(1) - wci, wci};
const Real wj[2] = {Real(1) - wfj, wfj};
const Real wk[2] = {Real(1) - wck, wck};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k) {
const Real w = wi[i] * wj[j] * wk[k];
mg[gidx + dX[i] + dY[j] + dZ[k]].y += w;
vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * vel.y;
vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * dot(cpy[idx], gpos + Vec3(i, j, k) - pos);
}
}
if (!vg.is3D())
return;
{ // w-face
const IndexInt gidx = ci * dX[1] + cj * dY[1] + fk * dZ[1];
const Vec3 gpos(ci + 0.5, cj + 0.5, fk);
const Real wi[2] = {Real(1) - wci, wci};
const Real wj[2] = {Real(1) - wcj, wcj};
const Real wk[2] = {Real(1) - wfk, wfk};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k) {
const Real w = wi[i] * wj[j] * wk[k];
mg[gidx + dX[i] + dY[j] + dZ[k]].z += w;
vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * vel.z;
vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * dot(cpz[idx], gpos + Vec3(i, j, k) - pos);
}
}
}
inline const BasicParticleSystem &getArg0()
{
return p;
}
typedef BasicParticleSystem type0;
inline MACGrid &getArg1()
{
return mg;
}
typedef MACGrid type1;
inline MACGrid &getArg2()
{
return vg;
}
typedef MACGrid type2;
inline const ParticleDataImpl<Vec3> &getArg3()
{
return vp;
}
typedef ParticleDataImpl<Vec3> type3;
inline const ParticleDataImpl<Vec3> &getArg4()
{
return cpx;
}
typedef ParticleDataImpl<Vec3> type4;
inline const ParticleDataImpl<Vec3> &getArg5()
{
return cpy;
}
typedef ParticleDataImpl<Vec3> type5;
inline const ParticleDataImpl<Vec3> &getArg6()
{
return cpz;
}
typedef ParticleDataImpl<Vec3> type6;
inline const ParticleDataImpl<int> *getArg7()
{
return ptype;
}
typedef ParticleDataImpl<int> type7;
inline const int &getArg8()
{
return exclude;
}
typedef int type8;
void runMessage()
{
debMsg("Executing kernel knApicMapLinearVec3ToMACGrid ", 3);
debMsg("Kernel range"
<< " size " << size << " ",
4);
};
void run()
{
const IndexInt _sz = size;
for (IndexInt i = 0; i < _sz; i++)
op(i, p, mg, vg, vp, cpx, cpy, cpz, ptype, exclude);
}
const BasicParticleSystem &p;
MACGrid &mg;
MACGrid &vg;
const ParticleDataImpl<Vec3> &vp;
const ParticleDataImpl<Vec3> &cpx;
const ParticleDataImpl<Vec3> &cpy;
const ParticleDataImpl<Vec3> &cpz;
const ParticleDataImpl<int> *ptype;
const int exclude;
};
void apicMapPartsToMAC(const FlagGrid &flags,
MACGrid &vel,
const BasicParticleSystem &parts,
const ParticleDataImpl<Vec3> &partVel,
const ParticleDataImpl<Vec3> &cpx,
const ParticleDataImpl<Vec3> &cpy,
const ParticleDataImpl<Vec3> &cpz,
MACGrid *mass = NULL,
const ParticleDataImpl<int> *ptype = NULL,
const int exclude = 0)
{
// affine map
// let's assume that the particle mass is constant, 1.0
const bool freeMass = !mass;
if (!mass)
mass = new MACGrid(flags.getParent());
else
mass->clear();
vel.clear();
knApicMapLinearVec3ToMACGrid(parts, *mass, vel, partVel, cpx, cpy, cpz, ptype, exclude);
mass->stomp(VECTOR_EPSILON);
vel.safeDivide(*mass);
if (freeMass)
delete mass;
}
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, "apicMapPartsToMAC", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock);
const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 2, &_lock);
const ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
"partVel", 3, &_lock);
const ParticleDataImpl<Vec3> &cpx = *_args.getPtr<ParticleDataImpl<Vec3>>("cpx", 4, &_lock);
const ParticleDataImpl<Vec3> &cpy = *_args.getPtr<ParticleDataImpl<Vec3>>("cpy", 5, &_lock);
const ParticleDataImpl<Vec3> &cpz = *_args.getPtr<ParticleDataImpl<Vec3>>("cpz", 6, &_lock);
MACGrid *mass = _args.getPtrOpt<MACGrid>("mass", 7, NULL, &_lock);
const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
"ptype", 8, NULL, &_lock);
const int exclude = _args.getOpt<int>("exclude", 9, 0, &_lock);
_retval = getPyNone();
apicMapPartsToMAC(flags, vel, parts, partVel, cpx, cpy, cpz, mass, ptype, exclude);
_args.check();
}
pbFinalizePlugin(parent, "apicMapPartsToMAC", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("apicMapPartsToMAC", e.what());
return 0;
}
}
static const Pb::Register _RP_apicMapPartsToMAC("", "apicMapPartsToMAC", _W_0);
extern "C" {
void PbRegister_apicMapPartsToMAC()
{
KEEP_UNUSED(_RP_apicMapPartsToMAC);
}
}
struct knApicMapLinearMACGridToVec3 : public KernelBase {
knApicMapLinearMACGridToVec3(ParticleDataImpl<Vec3> &vp,
ParticleDataImpl<Vec3> &cpx,
ParticleDataImpl<Vec3> &cpy,
ParticleDataImpl<Vec3> &cpz,
const BasicParticleSystem &p,
const MACGrid &vg,
const FlagGrid &flags,
const ParticleDataImpl<int> *ptype,
const int exclude)
: KernelBase(vp.size()),
vp(vp),
cpx(cpx),
cpy(cpy),
cpz(cpz),
p(p),
vg(vg),
flags(flags),
ptype(ptype),
exclude(exclude)
{
runMessage();
run();
}
inline void op(IndexInt idx,
ParticleDataImpl<Vec3> &vp,
ParticleDataImpl<Vec3> &cpx,
ParticleDataImpl<Vec3> &cpy,
ParticleDataImpl<Vec3> &cpz,
const BasicParticleSystem &p,
const MACGrid &vg,
const FlagGrid &flags,
const ParticleDataImpl<int> *ptype,
const int exclude) const
{
if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
return;
vp[idx] = cpx[idx] = cpy[idx] = cpz[idx] = Vec3(Real(0));
const IndexInt dX[2] = {0, vg.getStrideX()}, dY[2] = {0, vg.getStrideY()},
dZ[2] = {0, vg.getStrideZ()};
const Real gw[2] = {-Real(1), Real(1)};
const Vec3 &pos = p[idx].pos;
const IndexInt fi = static_cast<IndexInt>(pos.x), fj = static_cast<IndexInt>(pos.y),
fk = static_cast<IndexInt>(pos.z);
const IndexInt ci = static_cast<IndexInt>(pos.x - 0.5),
cj = static_cast<IndexInt>(pos.y - 0.5),
ck = static_cast<IndexInt>(pos.z - 0.5);
const Real wfi = clamp(pos.x - fi, Real(0), Real(1)),
wfj = clamp(pos.y - fj, Real(0), Real(1)),
wfk = clamp(pos.z - fk, Real(0), Real(1));
const Real wci = clamp(Real(pos.x - ci - 0.5), Real(0), Real(1)),
wcj = clamp(Real(pos.y - cj - 0.5), Real(0), Real(1)),
wck = clamp(Real(pos.z - ck - 0.5), Real(0), Real(1));
// TODO: check index for safety
{ // u
const IndexInt gidx = fi * dX[1] + cj * dY[1] + ck * dZ[1];
const Real wx[2] = {Real(1) - wfi, wfi};
const Real wy[2] = {Real(1) - wcj, wcj};
const Real wz[2] = {Real(1) - wck, wck};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k) {
const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
Real vgx = vg[vidx].x;
vp[idx].x += wx[i] * wy[j] * wz[k] * vgx;
cpx[idx].x += gw[i] * wy[j] * wz[k] * vgx;
cpx[idx].y += wx[i] * gw[j] * wz[k] * vgx;
cpx[idx].z += wx[i] * wy[j] * gw[k] * vgx;
}
}
{ // v
const IndexInt gidx = ci * dX[1] + fj * dY[1] + ck * dZ[1];
const Real wx[2] = {Real(1) - wci, wci};
const Real wy[2] = {Real(1) - wfj, wfj};
const Real wz[2] = {Real(1) - wck, wck};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k) {
const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
Real vgy = vg[vidx].y;
vp[idx].y += wx[i] * wy[j] * wz[k] * vgy;
cpy[idx].x += gw[i] * wy[j] * wz[k] * vgy;
cpy[idx].y += wx[i] * gw[j] * wz[k] * vgy;
cpy[idx].z += wx[i] * wy[j] * gw[k] * vgy;
}
}
if (!vg.is3D())
return;
{ // w
const IndexInt gidx = ci * dX[1] + cj * dY[1] + fk * dZ[1];
const Real wx[2] = {Real(1) - wci, wci};
const Real wy[2] = {Real(1) - wcj, wcj};
const Real wz[2] = {Real(1) - wfk, wfk};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k) {
const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
Real vgz = vg[vidx].z;
vp[idx].z += wx[i] * wy[j] * wz[k] * vgz;
cpz[idx].x += gw[i] * wy[j] * wz[k] * vgz;
cpz[idx].y += wx[i] * gw[j] * wz[k] * vgz;
cpz[idx].z += wx[i] * wy[j] * gw[k] * vgz;
}
}
}
inline ParticleDataImpl<Vec3> &getArg0()
{
return vp;
}
typedef ParticleDataImpl<Vec3> type0;
inline ParticleDataImpl<Vec3> &getArg1()
{
return cpx;
}
typedef ParticleDataImpl<Vec3> type1;
inline ParticleDataImpl<Vec3> &getArg2()
{
return cpy;
}
typedef ParticleDataImpl<Vec3> type2;
inline ParticleDataImpl<Vec3> &getArg3()
{
return cpz;
}
typedef ParticleDataImpl<Vec3> type3;
inline const BasicParticleSystem &getArg4()
{
return p;
}
typedef BasicParticleSystem type4;
inline const MACGrid &getArg5()
{
return vg;
}
typedef MACGrid type5;
inline const FlagGrid &getArg6()
{
return flags;
}
typedef FlagGrid type6;
inline const ParticleDataImpl<int> *getArg7()
{
return ptype;
}
typedef ParticleDataImpl<int> type7;
inline const int &getArg8()
{
return exclude;
}
typedef int type8;
void runMessage()
{
debMsg("Executing kernel knApicMapLinearMACGridToVec3 ", 3);
debMsg("Kernel range"
<< " size " << size << " ",
4);
};
void operator()(const tbb::blocked_range<IndexInt> &__r) const
{
for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
op(idx, vp, cpx, cpy, cpz, p, vg, flags, ptype, exclude);
}
void run()
{
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
}
ParticleDataImpl<Vec3> &vp;
ParticleDataImpl<Vec3> &cpx;
ParticleDataImpl<Vec3> &cpy;
ParticleDataImpl<Vec3> &cpz;
const BasicParticleSystem &p;
const MACGrid &vg;
const FlagGrid &flags;
const ParticleDataImpl<int> *ptype;
const int exclude;
};
void apicMapMACGridToParts(ParticleDataImpl<Vec3> &partVel,
ParticleDataImpl<Vec3> &cpx,
ParticleDataImpl<Vec3> &cpy,
ParticleDataImpl<Vec3> &cpz,
const BasicParticleSystem &parts,
const MACGrid &vel,
const FlagGrid &flags,
const ParticleDataImpl<int> *ptype = NULL,
const int exclude = 0)
{
knApicMapLinearMACGridToVec3(partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude);
}
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, "apicMapMACGridToParts", !noTiming);
PyObject *_retval = 0;
{
ArgLocker _lock;
ParticleDataImpl<Vec3> &partVel = *_args.getPtr<ParticleDataImpl<Vec3>>(
"partVel", 0, &_lock);
ParticleDataImpl<Vec3> &cpx = *_args.getPtr<ParticleDataImpl<Vec3>>("cpx", 1, &_lock);
ParticleDataImpl<Vec3> &cpy = *_args.getPtr<ParticleDataImpl<Vec3>>("cpy", 2, &_lock);
ParticleDataImpl<Vec3> &cpz = *_args.getPtr<ParticleDataImpl<Vec3>>("cpz", 3, &_lock);
const BasicParticleSystem &parts = *_args.getPtr<BasicParticleSystem>("parts", 4, &_lock);
const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 5, &_lock);
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 6, &_lock);
const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
"ptype", 7, NULL, &_lock);
const int exclude = _args.getOpt<int>("exclude", 8, 0, &_lock);
_retval = getPyNone();
apicMapMACGridToParts(partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude);
_args.check();
}
pbFinalizePlugin(parent, "apicMapMACGridToParts", !noTiming);
return _retval;
}
catch (std::exception &e) {
pbSetError("apicMapMACGridToParts", e.what());
return 0;
}
}
static const Pb::Register _RP_apicMapMACGridToParts("", "apicMapMACGridToParts", _W_1);
extern "C" {
void PbRegister_apicMapMACGridToParts()
{
KEEP_UNUSED(_RP_apicMapMACGridToParts);
}
}
} // namespace Manta