Sebastián Barschkis
635694c0ff
Mainly updated the Mantaflow version. It includes the new viscosity solver plugin based on the method from 'Accurate Viscous Free Surfaces for Buckling, Coiling, and Rotating Liquids' (Batty & Bridson). In the UI, this update adds a new 'Viscosity' section to the fluid modifier UI (liquid domains only). For now, there is a single 'strength' value to control the viscosity of liquids.
698 lines
21 KiB
C++
698 lines
21 KiB
C++
|
|
|
|
// 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
|
|
*
|
|
* Plugins for using vortex sheet meshes
|
|
*
|
|
******************************************************************************/
|
|
|
|
#include <iostream>
|
|
#include "vortexsheet.h"
|
|
#include "vortexpart.h"
|
|
#include "shapes.h"
|
|
#include "commonkernels.h"
|
|
#include "conjugategrad.h"
|
|
#include "randomstream.h"
|
|
#include "levelset.h"
|
|
|
|
using namespace std;
|
|
|
|
namespace Manta {
|
|
|
|
//! Mark area of mesh inside shape as fixed nodes.
|
|
//! Remove all other fixed nodes if 'exclusive' is set
|
|
|
|
void markAsFixed(Mesh &mesh, const Shape *shape, bool exclusive = true)
|
|
{
|
|
for (int i = 0; i < mesh.numNodes(); i++) {
|
|
if (shape->isInside(mesh.nodes(i).pos))
|
|
mesh.nodes(i).flags |= Mesh::NfFixed;
|
|
else if (exclusive)
|
|
mesh.nodes(i).flags &= ~Mesh::NfFixed;
|
|
}
|
|
}
|
|
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, "markAsFixed", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
Mesh &mesh = *_args.getPtr<Mesh>("mesh", 0, &_lock);
|
|
const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
|
|
bool exclusive = _args.getOpt<bool>("exclusive", 2, true, &_lock);
|
|
_retval = getPyNone();
|
|
markAsFixed(mesh, shape, exclusive);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "markAsFixed", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("markAsFixed", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_markAsFixed("", "markAsFixed", _W_0);
|
|
extern "C" {
|
|
void PbRegister_markAsFixed()
|
|
{
|
|
KEEP_UNUSED(_RP_markAsFixed);
|
|
}
|
|
}
|
|
|
|
//! Adapt texture coordinates of mesh inside shape
|
|
//! to obtain an effective inflow effect
|
|
|
|
void texcoordInflow(VortexSheetMesh &mesh, const Shape *shape, const MACGrid &vel)
|
|
{
|
|
static Vec3 t0 = Vec3::Zero;
|
|
|
|
// get mean velocity
|
|
int cnt = 0;
|
|
Vec3 meanV(0.0);
|
|
FOR_IJK(vel)
|
|
{
|
|
if (shape->isInsideGrid(i, j, k)) {
|
|
cnt++;
|
|
meanV += vel.getCentered(i, j, k);
|
|
}
|
|
}
|
|
meanV /= (Real)cnt;
|
|
t0 -= mesh.getParent()->getDt() * meanV;
|
|
mesh.setReferenceTexOffset(t0);
|
|
|
|
// apply mean velocity
|
|
for (int i = 0; i < mesh.numNodes(); i++) {
|
|
if (shape->isInside(mesh.nodes(i).pos)) {
|
|
Vec3 tc = mesh.nodes(i).pos + t0;
|
|
mesh.tex1(i) = tc;
|
|
mesh.tex2(i) = tc;
|
|
}
|
|
}
|
|
}
|
|
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, "texcoordInflow", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
|
|
const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
|
|
const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 2, &_lock);
|
|
_retval = getPyNone();
|
|
texcoordInflow(mesh, shape, vel);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "texcoordInflow", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("texcoordInflow", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_texcoordInflow("", "texcoordInflow", _W_1);
|
|
extern "C" {
|
|
void PbRegister_texcoordInflow()
|
|
{
|
|
KEEP_UNUSED(_RP_texcoordInflow);
|
|
}
|
|
}
|
|
|
|
;
|
|
|
|
//! Init smoke density values of the mesh surface inside source shape
|
|
|
|
void meshSmokeInflow(VortexSheetMesh &mesh, const Shape *shape, Real amount)
|
|
{
|
|
for (int t = 0; t < mesh.numTris(); t++) {
|
|
if (shape->isInside(mesh.getFaceCenter(t)))
|
|
mesh.sheet(t).smokeAmount = amount;
|
|
}
|
|
}
|
|
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, "meshSmokeInflow", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
|
|
const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
|
|
Real amount = _args.get<Real>("amount", 2, &_lock);
|
|
_retval = getPyNone();
|
|
meshSmokeInflow(mesh, shape, amount);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "meshSmokeInflow", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("meshSmokeInflow", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_meshSmokeInflow("", "meshSmokeInflow", _W_2);
|
|
extern "C" {
|
|
void PbRegister_meshSmokeInflow()
|
|
{
|
|
KEEP_UNUSED(_RP_meshSmokeInflow);
|
|
}
|
|
}
|
|
|
|
struct KnAcceleration : public KernelBase {
|
|
KnAcceleration(MACGrid &a, const MACGrid &v1, const MACGrid &v0, const Real idt)
|
|
: KernelBase(&a, 0), a(a), v1(v1), v0(v0), idt(idt)
|
|
{
|
|
runMessage();
|
|
run();
|
|
}
|
|
inline void op(
|
|
IndexInt idx, MACGrid &a, const MACGrid &v1, const MACGrid &v0, const Real idt) const
|
|
{
|
|
a[idx] = (v1[idx] - v0[idx]) * idt;
|
|
}
|
|
inline MACGrid &getArg0()
|
|
{
|
|
return a;
|
|
}
|
|
typedef MACGrid type0;
|
|
inline const MACGrid &getArg1()
|
|
{
|
|
return v1;
|
|
}
|
|
typedef MACGrid type1;
|
|
inline const MACGrid &getArg2()
|
|
{
|
|
return v0;
|
|
}
|
|
typedef MACGrid type2;
|
|
inline const Real &getArg3()
|
|
{
|
|
return idt;
|
|
}
|
|
typedef Real type3;
|
|
void runMessage()
|
|
{
|
|
debMsg("Executing kernel KnAcceleration ", 3);
|
|
debMsg("Kernel range"
|
|
<< " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
|
|
4);
|
|
};
|
|
void operator()(const tbb::blocked_range<IndexInt> &__r) const
|
|
{
|
|
for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
|
|
op(idx, a, v1, v0, idt);
|
|
}
|
|
void run()
|
|
{
|
|
tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
|
|
}
|
|
MACGrid &a;
|
|
const MACGrid &v1;
|
|
const MACGrid &v0;
|
|
const Real idt;
|
|
};
|
|
|
|
//! Add vorticity to vortex sheets based on buoyancy
|
|
|
|
void vorticitySource(VortexSheetMesh &mesh,
|
|
Vec3 gravity,
|
|
const MACGrid *vel = nullptr,
|
|
const MACGrid *velOld = nullptr,
|
|
Real scale = 0.1,
|
|
Real maxAmount = 0,
|
|
Real mult = 1.0)
|
|
{
|
|
Real dt = mesh.getParent()->getDt();
|
|
Real dx = mesh.getParent()->getDx();
|
|
MACGrid acceleration(mesh.getParent());
|
|
if (vel)
|
|
KnAcceleration(acceleration, *vel, *velOld, 1.0 / dt);
|
|
const Real A = -1.0;
|
|
Real maxV = 0, meanV = 0;
|
|
|
|
for (int t = 0; t < mesh.numTris(); t++) {
|
|
Vec3 fn = mesh.getFaceNormal(t);
|
|
Vec3 source;
|
|
if (vel) {
|
|
Vec3 a = acceleration.getInterpolated(mesh.getFaceCenter(t));
|
|
source = A * cross(fn, a - gravity) * scale;
|
|
}
|
|
else {
|
|
source = A * cross(fn, -gravity) * scale;
|
|
}
|
|
|
|
if (mesh.isTriangleFixed(t))
|
|
source = 0;
|
|
|
|
mesh.sheet(t).vorticity *= mult;
|
|
mesh.sheet(t).vorticity += dt * source / dx;
|
|
// upper limit
|
|
Real v = norm(mesh.sheet(t).vorticity);
|
|
if (maxAmount > 0 && v > maxAmount)
|
|
mesh.sheet(t).vorticity *= maxAmount / v;
|
|
|
|
// stats
|
|
if (v > maxV)
|
|
maxV = v;
|
|
meanV += v;
|
|
}
|
|
|
|
cout << "vorticity: max " << maxV << " / mean " << meanV / mesh.numTris() << endl;
|
|
}
|
|
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, "vorticitySource", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
|
|
Vec3 gravity = _args.get<Vec3>("gravity", 1, &_lock);
|
|
const MACGrid *vel = _args.getPtrOpt<MACGrid>("vel", 2, nullptr, &_lock);
|
|
const MACGrid *velOld = _args.getPtrOpt<MACGrid>("velOld", 3, nullptr, &_lock);
|
|
Real scale = _args.getOpt<Real>("scale", 4, 0.1, &_lock);
|
|
Real maxAmount = _args.getOpt<Real>("maxAmount", 5, 0, &_lock);
|
|
Real mult = _args.getOpt<Real>("mult", 6, 1.0, &_lock);
|
|
_retval = getPyNone();
|
|
vorticitySource(mesh, gravity, vel, velOld, scale, maxAmount, mult);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "vorticitySource", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("vorticitySource", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_vorticitySource("", "vorticitySource", _W_3);
|
|
extern "C" {
|
|
void PbRegister_vorticitySource()
|
|
{
|
|
KEEP_UNUSED(_RP_vorticitySource);
|
|
}
|
|
}
|
|
|
|
void smoothVorticity(VortexSheetMesh &mesh, int iter = 1, Real sigma = 0.2, Real alpha = 0.8)
|
|
{
|
|
const Real mult = -0.5 / sigma / sigma;
|
|
|
|
// pre-calculate positions and weights
|
|
vector<Vec3> vort(mesh.numTris()), pos(mesh.numTris());
|
|
vector<Real> weights(3 * mesh.numTris());
|
|
vector<int> index(3 * mesh.numTris());
|
|
for (int i = 0; i < mesh.numTris(); i++) {
|
|
pos[i] = mesh.getFaceCenter(i);
|
|
mesh.sheet(i).vorticitySmoothed = mesh.sheet(i).vorticity;
|
|
}
|
|
for (int i = 0; i < mesh.numTris(); i++) {
|
|
for (int c = 0; c < 3; c++) {
|
|
int oc = mesh.corners(i, c).opposite;
|
|
if (oc >= 0) {
|
|
int t = mesh.corners(oc).tri;
|
|
weights[3 * i + c] = exp(normSquare(pos[t] - pos[i]) * mult);
|
|
index[3 * i + c] = t;
|
|
}
|
|
else {
|
|
weights[3 * i + c] = 0;
|
|
index[3 * i + c] = 0;
|
|
}
|
|
}
|
|
}
|
|
|
|
for (int it = 0; it < iter; ++it) {
|
|
// first, preload
|
|
for (int i = 0; i < mesh.numTris(); i++)
|
|
vort[i] = mesh.sheet(i).vorticitySmoothed;
|
|
|
|
for (int i = 0, idx = 0; i < mesh.numTris(); i++) {
|
|
// loop over adjacent tris
|
|
Real sum = 1.0f;
|
|
Vec3 v = vort[i];
|
|
for (int c = 0; c < 3; c++, idx++) {
|
|
Real w = weights[index[idx]];
|
|
v += w * vort[index[idx]];
|
|
sum += w;
|
|
}
|
|
mesh.sheet(i).vorticitySmoothed = v / sum;
|
|
}
|
|
}
|
|
for (int i = 0; i < mesh.numTris(); i++)
|
|
mesh.sheet(i).vorticitySmoothed *= alpha;
|
|
}
|
|
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, "smoothVorticity", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
|
|
int iter = _args.getOpt<int>("iter", 1, 1, &_lock);
|
|
Real sigma = _args.getOpt<Real>("sigma", 2, 0.2, &_lock);
|
|
Real alpha = _args.getOpt<Real>("alpha", 3, 0.8, &_lock);
|
|
_retval = getPyNone();
|
|
smoothVorticity(mesh, iter, sigma, alpha);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "smoothVorticity", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("smoothVorticity", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_smoothVorticity("", "smoothVorticity", _W_4);
|
|
extern "C" {
|
|
void PbRegister_smoothVorticity()
|
|
{
|
|
KEEP_UNUSED(_RP_smoothVorticity);
|
|
}
|
|
}
|
|
|
|
//! Seed Vortex Particles inside shape with K41 characteristics
|
|
void VPseedK41(VortexParticleSystem &system,
|
|
const Shape *shape,
|
|
Real strength = 0,
|
|
Real sigma0 = 0.2,
|
|
Real sigma1 = 1.0,
|
|
Real probability = 1.0,
|
|
Real N = 3.0)
|
|
{
|
|
Grid<Real> temp(system.getParent());
|
|
const Real dt = system.getParent()->getDt();
|
|
static RandomStream rand(3489572);
|
|
Real s0 = pow((Real)sigma0, (Real)(-N + 1.0));
|
|
Real s1 = pow((Real)sigma1, (Real)(-N + 1.0));
|
|
|
|
FOR_IJK(temp)
|
|
{
|
|
if (shape->isInsideGrid(i, j, k)) {
|
|
if (rand.getReal() < probability * dt) {
|
|
Real p = rand.getReal();
|
|
Real sigma = pow((1.0 - p) * s0 + p * s1, 1. / (-N + 1.0));
|
|
Vec3 randDir(rand.getReal(), rand.getReal(), rand.getReal());
|
|
Vec3 posUpd(i + rand.getReal(), j + rand.getReal(), k + rand.getReal());
|
|
normalize(randDir);
|
|
Vec3 vorticity = randDir * strength * pow((Real)sigma, (Real)(-10. / 6. + N / 2.0));
|
|
system.add(VortexParticleData(posUpd, vorticity, sigma));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
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, "VPseedK41", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
VortexParticleSystem &system = *_args.getPtr<VortexParticleSystem>("system", 0, &_lock);
|
|
const Shape *shape = _args.getPtr<Shape>("shape", 1, &_lock);
|
|
Real strength = _args.getOpt<Real>("strength", 2, 0, &_lock);
|
|
Real sigma0 = _args.getOpt<Real>("sigma0", 3, 0.2, &_lock);
|
|
Real sigma1 = _args.getOpt<Real>("sigma1", 4, 1.0, &_lock);
|
|
Real probability = _args.getOpt<Real>("probability", 5, 1.0, &_lock);
|
|
Real N = _args.getOpt<Real>("N", 6, 3.0, &_lock);
|
|
_retval = getPyNone();
|
|
VPseedK41(system, shape, strength, sigma0, sigma1, probability, N);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "VPseedK41", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("VPseedK41", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_VPseedK41("", "VPseedK41", _W_5);
|
|
extern "C" {
|
|
void PbRegister_VPseedK41()
|
|
{
|
|
KEEP_UNUSED(_RP_VPseedK41);
|
|
}
|
|
}
|
|
|
|
//! Vortex-in-cell integration
|
|
|
|
void VICintegration(VortexSheetMesh &mesh,
|
|
Real sigma,
|
|
Grid<Vec3> &vel,
|
|
const FlagGrid &flags,
|
|
Grid<Vec3> *vorticity = nullptr,
|
|
Real cgMaxIterFac = 1.5,
|
|
Real cgAccuracy = 1e-3,
|
|
Real scale = 0.01,
|
|
int precondition = 0)
|
|
{
|
|
|
|
MuTime t0;
|
|
const Real fac = 16.0; // experimental factor to balance out regularization
|
|
|
|
// if no vort grid is given, use a temporary one
|
|
Grid<Vec3> vortTemp(mesh.getParent());
|
|
Grid<Vec3> &vort = (vorticity) ? (*vorticity) : (vortTemp);
|
|
vort.clear();
|
|
|
|
// map vorticity to grid using Peskin kernel
|
|
int sgi = ceil(sigma);
|
|
Real pkfac = M_PI / sigma;
|
|
const int numTris = mesh.numTris();
|
|
for (int t = 0; t < numTris; t++) {
|
|
Vec3 pos = mesh.getFaceCenter(t);
|
|
Vec3 v = mesh.sheet(t).vorticity * mesh.getFaceArea(t) * fac;
|
|
|
|
// inner kernel
|
|
// first, summate
|
|
Real sum = 0;
|
|
for (int i = -sgi; i < sgi; i++) {
|
|
if (pos.x + i < 0 || (int)pos.x + i >= vort.getSizeX())
|
|
continue;
|
|
for (int j = -sgi; j < sgi; j++) {
|
|
if (pos.y + j < 0 || (int)pos.y + j >= vort.getSizeY())
|
|
continue;
|
|
for (int k = -sgi; k < sgi; k++) {
|
|
if (pos.z + k < 0 || (int)pos.z + k >= vort.getSizeZ())
|
|
continue;
|
|
Vec3i cell(pos.x + i, pos.y + j, pos.z + k);
|
|
if (!flags.isFluid(cell))
|
|
continue;
|
|
Vec3 d = pos -
|
|
Vec3(i + 0.5 + floor(pos.x), j + 0.5 + floor(pos.y), k + 0.5 + floor(pos.z));
|
|
Real dl = norm(d);
|
|
if (dl > sigma)
|
|
continue;
|
|
// precalc Peskin kernel
|
|
sum += 1.0 + cos(dl * pkfac);
|
|
}
|
|
}
|
|
}
|
|
// then, apply normalized kernel
|
|
Real wnorm = 1.0 / sum;
|
|
for (int i = -sgi; i < sgi; i++) {
|
|
if (pos.x + i < 0 || (int)pos.x + i >= vort.getSizeX())
|
|
continue;
|
|
for (int j = -sgi; j < sgi; j++) {
|
|
if (pos.y + j < 0 || (int)pos.y + j >= vort.getSizeY())
|
|
continue;
|
|
for (int k = -sgi; k < sgi; k++) {
|
|
if (pos.z + k < 0 || (int)pos.z + k >= vort.getSizeZ())
|
|
continue;
|
|
Vec3i cell(pos.x + i, pos.y + j, pos.z + k);
|
|
if (!flags.isFluid(cell))
|
|
continue;
|
|
Vec3 d = pos -
|
|
Vec3(i + 0.5 + floor(pos.x), j + 0.5 + floor(pos.y), k + 0.5 + floor(pos.z));
|
|
Real dl = norm(d);
|
|
if (dl > sigma)
|
|
continue;
|
|
Real w = (1.0 + cos(dl * pkfac)) * wnorm;
|
|
vort(cell) += v * w;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Prepare grids for poisson solve
|
|
Grid<Vec3> vortexCurl(mesh.getParent());
|
|
Grid<Real> rhs(mesh.getParent());
|
|
Grid<Real> solution(mesh.getParent());
|
|
Grid<Real> residual(mesh.getParent());
|
|
Grid<Real> search(mesh.getParent());
|
|
Grid<Real> temp1(mesh.getParent());
|
|
Grid<Real> A0(mesh.getParent());
|
|
Grid<Real> Ai(mesh.getParent());
|
|
Grid<Real> Aj(mesh.getParent());
|
|
Grid<Real> Ak(mesh.getParent());
|
|
Grid<Real> pca0(mesh.getParent());
|
|
Grid<Real> pca1(mesh.getParent());
|
|
Grid<Real> pca2(mesh.getParent());
|
|
Grid<Real> pca3(mesh.getParent());
|
|
|
|
MakeLaplaceMatrix(flags, A0, Ai, Aj, Ak);
|
|
CurlOp(vort, vortexCurl);
|
|
|
|
// Solve vector poisson equation
|
|
for (int c = 0; c < 3; c++) {
|
|
// construct rhs
|
|
if (vel.getType() & GridBase::TypeMAC)
|
|
GetShiftedComponent(vortexCurl, rhs, c);
|
|
else
|
|
GetComponent(vortexCurl, rhs, c);
|
|
|
|
// prepare CG solver
|
|
const int maxIter = (int)(cgMaxIterFac * vel.getSize().max());
|
|
vector<Grid<Real> *> matA{&A0, &Ai, &Aj, &Ak};
|
|
|
|
GridCgInterface *gcg = new GridCg<ApplyMatrix>(
|
|
solution, rhs, residual, search, flags, temp1, matA);
|
|
gcg->setAccuracy(cgAccuracy);
|
|
gcg->setUseL2Norm(true);
|
|
gcg->setICPreconditioner(
|
|
(GridCgInterface::PreconditionType)precondition, &pca0, &pca1, &pca2, &pca3);
|
|
|
|
// iterations
|
|
for (int iter = 0; iter < maxIter; iter++) {
|
|
if (!gcg->iterate())
|
|
iter = maxIter;
|
|
}
|
|
debMsg("VICintegration CG iterations:" << gcg->getIterations() << ", res:" << gcg->getSigma(),
|
|
1);
|
|
delete gcg;
|
|
|
|
// copy back
|
|
solution *= scale;
|
|
SetComponent(vel, solution, c);
|
|
}
|
|
}
|
|
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, "VICintegration", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
VortexSheetMesh &mesh = *_args.getPtr<VortexSheetMesh>("mesh", 0, &_lock);
|
|
Real sigma = _args.get<Real>("sigma", 1, &_lock);
|
|
Grid<Vec3> &vel = *_args.getPtr<Grid<Vec3>>("vel", 2, &_lock);
|
|
const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock);
|
|
Grid<Vec3> *vorticity = _args.getPtrOpt<Grid<Vec3>>("vorticity", 4, nullptr, &_lock);
|
|
Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 5, 1.5, &_lock);
|
|
Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 6, 1e-3, &_lock);
|
|
Real scale = _args.getOpt<Real>("scale", 7, 0.01, &_lock);
|
|
int precondition = _args.getOpt<int>("precondition", 8, 0, &_lock);
|
|
_retval = getPyNone();
|
|
VICintegration(
|
|
mesh, sigma, vel, flags, vorticity, cgMaxIterFac, cgAccuracy, scale, precondition);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "VICintegration", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("VICintegration", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_VICintegration("", "VICintegration", _W_6);
|
|
extern "C" {
|
|
void PbRegister_VICintegration()
|
|
{
|
|
KEEP_UNUSED(_RP_VICintegration);
|
|
}
|
|
}
|
|
|
|
//! Obtain density field from levelset with linear gradient of size sigma over the interface
|
|
void densityFromLevelset(const LevelsetGrid &phi,
|
|
Grid<Real> &density,
|
|
Real value = 1.0,
|
|
Real sigma = 1.0)
|
|
{
|
|
FOR_IJK(phi)
|
|
{
|
|
// remove boundary
|
|
if (i < 2 || j < 2 || k < 2 || i >= phi.getSizeX() - 2 || j >= phi.getSizeY() - 2 ||
|
|
k >= phi.getSizeZ() - 2)
|
|
density(i, j, k) = 0;
|
|
else if (phi(i, j, k) < -sigma)
|
|
density(i, j, k) = value;
|
|
else if (phi(i, j, k) > sigma)
|
|
density(i, j, k) = 0;
|
|
else
|
|
density(i, j, k) = clamp(
|
|
(Real)(0.5 * value / sigma * (1.0 - phi(i, j, k))), (Real)0.0, value);
|
|
}
|
|
}
|
|
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, "densityFromLevelset", !noTiming);
|
|
PyObject *_retval = nullptr;
|
|
{
|
|
ArgLocker _lock;
|
|
const LevelsetGrid &phi = *_args.getPtr<LevelsetGrid>("phi", 0, &_lock);
|
|
Grid<Real> &density = *_args.getPtr<Grid<Real>>("density", 1, &_lock);
|
|
Real value = _args.getOpt<Real>("value", 2, 1.0, &_lock);
|
|
Real sigma = _args.getOpt<Real>("sigma", 3, 1.0, &_lock);
|
|
_retval = getPyNone();
|
|
densityFromLevelset(phi, density, value, sigma);
|
|
_args.check();
|
|
}
|
|
pbFinalizePlugin(parent, "densityFromLevelset", !noTiming);
|
|
return _retval;
|
|
}
|
|
catch (std::exception &e) {
|
|
pbSetError("densityFromLevelset", e.what());
|
|
return 0;
|
|
}
|
|
}
|
|
static const Pb::Register _RP_densityFromLevelset("", "densityFromLevelset", _W_7);
|
|
extern "C" {
|
|
void PbRegister_densityFromLevelset()
|
|
{
|
|
KEEP_UNUSED(_RP_densityFromLevelset);
|
|
}
|
|
}
|
|
|
|
} // namespace Manta
|