diff options
Diffstat (limited to 'extern/mantaflow/preprocessed/plugin/apic.cpp')
| -rw-r--r-- | extern/mantaflow/preprocessed/plugin/apic.cpp | 496 |
1 files changed, 496 insertions, 0 deletions
diff --git a/extern/mantaflow/preprocessed/plugin/apic.cpp b/extern/mantaflow/preprocessed/plugin/apic.cpp new file mode 100644 index 00000000000..6ff893014c9 --- /dev/null +++ b/extern/mantaflow/preprocessed/plugin/apic.cpp @@ -0,0 +1,496 @@ + + +// 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 |