diff options
Diffstat (limited to 'extern/mantaflow/preprocessed/vortexpart.cpp')
| -rw-r--r-- | extern/mantaflow/preprocessed/vortexpart.cpp | 251 |
1 files changed, 251 insertions, 0 deletions
diff --git a/extern/mantaflow/preprocessed/vortexpart.cpp b/extern/mantaflow/preprocessed/vortexpart.cpp new file mode 100644 index 00000000000..0eba2743ee8 --- /dev/null +++ b/extern/mantaflow/preprocessed/vortexpart.cpp @@ -0,0 +1,251 @@ + + +// 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 + * + * Vortex particles + * (warning, the vortex methods are currently experimental, and not fully supported!) + * + ******************************************************************************/ + +#include "vortexpart.h" +#include "integrator.h" +#include "mesh.h" + +using namespace std; +namespace Manta { + +// vortex particle effect: (cyl coord around wp) +// u = -|wp|*rho*exp( (-rho^2-z^2)/(2sigma^2) ) e_phi +inline Vec3 VortexKernel(const Vec3 &p, const vector<VortexParticleData> &vp, Real scale) +{ + Vec3 u(0.0); + for (size_t i = 0; i < vp.size(); i++) { + if (vp[i].flag & ParticleBase::PDELETE) + continue; + + // cutoff radius + const Vec3 r = p - vp[i].pos; + const Real rlen2 = normSquare(r); + const Real sigma2 = square(vp[i].sigma); + if (rlen2 > 6.0 * sigma2 || rlen2 < 1e-8) + continue; + + // split vortex strength + Vec3 vortNorm = vp[i].vorticity; + Real strength = normalize(vortNorm) * scale; + + // transform in cylinder coordinate system + const Real rlen = sqrt(rlen2); + const Real z = dot(r, vortNorm); + const Vec3 ePhi = cross(r, vortNorm) / rlen; + const Real rho2 = rlen2 - z * z; + + Real vortex = 0; + if (rho2 > 1e-10) { + // evaluate Kernel + vortex = strength * sqrt(rho2) * exp(rlen2 * -0.5 / sigma2); + } + u += vortex * ePhi; + } + return u; +} + +struct _KnVpAdvectMesh : public KernelBase { + _KnVpAdvectMesh(const KernelBase &base, + vector<Node> &nodes, + const vector<VortexParticleData> &vp, + Real scale, + vector<Vec3> &u) + : KernelBase(base), nodes(nodes), vp(vp), scale(scale), u(u) + { + } + inline void op(IndexInt idx, + vector<Node> &nodes, + const vector<VortexParticleData> &vp, + Real scale, + vector<Vec3> &u) const + { + if (nodes[idx].flags & Mesh::NfFixed) + u[idx] = 0.0; + else + u[idx] = VortexKernel(nodes[idx].pos, vp, scale); + } + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++) + op(idx, nodes, vp, scale, u); + } + void run() + { + tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this); + } + vector<Node> &nodes; + const vector<VortexParticleData> &vp; + Real scale; + vector<Vec3> &u; +}; +struct KnVpAdvectMesh : public KernelBase { + KnVpAdvectMesh(vector<Node> &nodes, const vector<VortexParticleData> &vp, Real scale) + : KernelBase(nodes.size()), + _inner(KernelBase(nodes.size()), nodes, vp, scale, u), + nodes(nodes), + vp(vp), + scale(scale), + u((size)) + { + runMessage(); + run(); + } + void run() + { + _inner.run(); + } + inline operator vector<Vec3>() + { + return u; + } + inline vector<Vec3> &getRet() + { + return u; + } + inline vector<Node> &getArg0() + { + return nodes; + } + typedef vector<Node> type0; + inline const vector<VortexParticleData> &getArg1() + { + return vp; + } + typedef vector<VortexParticleData> type1; + inline Real &getArg2() + { + return scale; + } + typedef Real type2; + void runMessage() + { + debMsg("Executing kernel KnVpAdvectMesh ", 3); + debMsg("Kernel range" + << " size " << size << " ", + 4); + }; + _KnVpAdvectMesh _inner; + vector<Node> &nodes; + const vector<VortexParticleData> &vp; + Real scale; + vector<Vec3> u; +}; + +struct _KnVpAdvectSelf : public KernelBase { + _KnVpAdvectSelf(const KernelBase &base, + vector<VortexParticleData> &vp, + Real scale, + vector<Vec3> &u) + : KernelBase(base), vp(vp), scale(scale), u(u) + { + } + inline void op(IndexInt idx, vector<VortexParticleData> &vp, Real scale, vector<Vec3> &u) const + { + if (vp[idx].flag & ParticleBase::PDELETE) + u[idx] = 0.0; + else + u[idx] = VortexKernel(vp[idx].pos, vp, scale); + } + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++) + op(idx, vp, scale, u); + } + void run() + { + tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this); + } + vector<VortexParticleData> &vp; + Real scale; + vector<Vec3> &u; +}; +struct KnVpAdvectSelf : public KernelBase { + KnVpAdvectSelf(vector<VortexParticleData> &vp, Real scale) + : KernelBase(vp.size()), + _inner(KernelBase(vp.size()), vp, scale, u), + vp(vp), + scale(scale), + u((size)) + { + runMessage(); + run(); + } + void run() + { + _inner.run(); + } + inline operator vector<Vec3>() + { + return u; + } + inline vector<Vec3> &getRet() + { + return u; + } + inline vector<VortexParticleData> &getArg0() + { + return vp; + } + typedef vector<VortexParticleData> type0; + inline Real &getArg1() + { + return scale; + } + typedef Real type1; + void runMessage() + { + debMsg("Executing kernel KnVpAdvectSelf ", 3); + debMsg("Kernel range" + << " size " << size << " ", + 4); + }; + _KnVpAdvectSelf _inner; + vector<VortexParticleData> &vp; + Real scale; + vector<Vec3> u; +}; + +VortexParticleSystem::VortexParticleSystem(FluidSolver *parent) + : ParticleSystem<VortexParticleData>(parent) +{ +} + +void VortexParticleSystem::advectSelf(Real scale, int integrationMode) +{ + KnVpAdvectSelf kernel(mData, scale * getParent()->getDt()); + integratePointSet(kernel, integrationMode); +} + +void VortexParticleSystem::applyToMesh(Mesh &mesh, Real scale, int integrationMode) +{ + KnVpAdvectMesh kernel(mesh.getNodeData(), mData, scale * getParent()->getDt()); + integratePointSet(kernel, integrationMode); +} + +ParticleBase *VortexParticleSystem::clone() +{ + VortexParticleSystem *nm = new VortexParticleSystem(getParent()); + compress(); + + nm->mData = mData; + nm->setName(getName()); + return nm; +} + +} // namespace Manta |