Page MenuHome
Differential D3850 Diff 19901 extern/mantaflow/preprocessed/vortexpart.cpp

Changeset View

Standalone View

View Options

extern/mantaflow/preprocessed/vortexpart.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
*
* 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