Page MenuHome
Differential D3850 Diff 17619 intern/mantaflow/intern/manta_develop/preprocessed/omp/plugin/flip.cpp

Changeset View

Standalone View

View Options

intern/mantaflow/intern/manta_develop/preprocessed/omp/plugin/flip.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
*
* FLIP (fluid implicit particles)
* for use with particle data fields
*
******************************************************************************/
#include "particle.h"
#include "grid.h"
#include "commonkernels.h"
#include "randomstream.h"
#include "levelset.h"
#include "shapes.h"
#include "matrixbase.h"
using namespace std;
namespace Manta {
// init
//! note - this is a simplified version , sampleLevelsetWithParticles has more functionality
void sampleFlagsWithParticles(const FlagGrid& flags, BasicParticleSystem& parts, const int discretization, const Real randomness) {
const bool is3D = flags.is3D();
const Real jlen = randomness / discretization;
const Vec3 disp (1.0 / discretization, 1.0 / discretization, 1.0/discretization);
RandomStream mRand(9832);
FOR_IJK_BND(flags, 0) {
if ( flags.isObstacle(i,j,k) ) continue;
if ( flags.isFluid(i,j,k) ) {
const Vec3 pos (i,j,k);
for (int dk=0; dk<(is3D ? discretization : 1); dk++)
for (int dj=0; dj<discretization; dj++)
for (int di=0; di<discretization; di++) {
Vec3 subpos = pos + disp * Vec3(0.5+di, 0.5+dj, 0.5+dk);
subpos += jlen * (Vec3(1,1,1) - 2.0 * mRand.getVec3());
if(!is3D) subpos[2] = 0.5;
parts.addBuffered(subpos);
}
}
}
parts.insertBufferedParticles();
} 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, "sampleFlagsWithParticles" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",1,&_lock); const int discretization = _args.get<int >("discretization",2,&_lock); const Real randomness = _args.get<Real >("randomness",3,&_lock); _retval = getPyNone(); sampleFlagsWithParticles(flags,parts,discretization,randomness); _args.check(); }
pbFinalizePlugin(parent,"sampleFlagsWithParticles", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("sampleFlagsWithParticles",e.what()); return 0; }
}
static const Pb::Register _RP_sampleFlagsWithParticles ("","sampleFlagsWithParticles",_W_0);
extern "C" {
void PbRegister_sampleFlagsWithParticles() {
KEEP_UNUSED(_RP_sampleFlagsWithParticles); }
}
//! sample a level set with particles, use reset to clear the particle buffer,
//! and skipEmpty for a continuous inflow (in the latter case, only empty cells will
//! be re-filled once they empty when calling sampleLevelsetWithParticles during
//! the main loop).
void sampleLevelsetWithParticles(const LevelsetGrid& phi, const FlagGrid& flags, BasicParticleSystem& parts, const int discretization, const Real randomness, const bool reset=false, const bool refillEmpty=false, const int particleFlag=-1) {
const bool is3D = phi.is3D();
const Real jlen = randomness / discretization;
const Vec3 disp (1.0 / discretization, 1.0 / discretization, 1.0/discretization);
RandomStream mRand(9832);
if(reset) {
parts.clear();
parts.doCompress();
}
FOR_IJK_BND(phi, 0) {
if ( flags.isObstacle(i,j,k) ) continue;
if ( refillEmpty && flags.isFluid(i,j,k) ) continue;
if ( phi(i,j,k) < 1.733 ) {
const Vec3 pos (i,j,k);
for (int dk=0; dk<(is3D ? discretization : 1); dk++)
for (int dj=0; dj<discretization; dj++)
for (int di=0; di<discretization; di++) {
Vec3 subpos = pos + disp * Vec3(0.5+di, 0.5+dj, 0.5+dk);
subpos += jlen * (Vec3(1,1,1) - 2.0 * mRand.getVec3());
if(!is3D) subpos[2] = 0.5;
if( phi.getInterpolated(subpos) > 0. ) continue;
if(particleFlag < 0){
parts.addBuffered(subpos);
}
else{
parts.addBuffered(subpos, particleFlag);
}
}
}
}
parts.insertBufferedParticles();
} 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, "sampleLevelsetWithParticles" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const LevelsetGrid& phi = *_args.getPtr<LevelsetGrid >("phi",0,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",1,&_lock); BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",2,&_lock); const int discretization = _args.get<int >("discretization",3,&_lock); const Real randomness = _args.get<Real >("randomness",4,&_lock); const bool reset = _args.getOpt<bool >("reset",5,false,&_lock); const bool refillEmpty = _args.getOpt<bool >("refillEmpty",6,false,&_lock); const int particleFlag = _args.getOpt<int >("particleFlag",7,-1,&_lock); _retval = getPyNone(); sampleLevelsetWithParticles(phi,flags,parts,discretization,randomness,reset,refillEmpty,particleFlag); _args.check(); }
pbFinalizePlugin(parent,"sampleLevelsetWithParticles", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("sampleLevelsetWithParticles",e.what()); return 0; }
}
static const Pb::Register _RP_sampleLevelsetWithParticles ("","sampleLevelsetWithParticles",_W_1);
extern "C" {
void PbRegister_sampleLevelsetWithParticles() {
KEEP_UNUSED(_RP_sampleLevelsetWithParticles); }
}
//! sample a shape with particles, use reset to clear the particle buffer,
//! and skipEmpty for a continuous inflow (in the latter case, only empty cells will
//! be re-filled once they empty when calling sampleShapeWithParticles during
//! the main loop).
void sampleShapeWithParticles(const Shape& shape, const FlagGrid& flags, BasicParticleSystem& parts, const int discretization, const Real randomness, const bool reset=false, const bool refillEmpty=false, const LevelsetGrid *exclude=NULL) {
const bool is3D = flags.is3D();
const Real jlen = randomness / discretization;
const Vec3 disp (1.0 / discretization, 1.0 / discretization, 1.0/discretization);
RandomStream mRand(9832);
if(reset) {
parts.clear();
parts.doCompress();
}
FOR_IJK_BND(flags, 0) {
if ( flags.isObstacle(i,j,k) ) continue;
if ( refillEmpty && flags.isFluid(i,j,k) ) continue;
const Vec3 pos (i,j,k);
for (int dk=0; dk<(is3D ? discretization : 1); dk++)
for (int dj=0; dj<discretization; dj++)
for (int di=0; di<discretization; di++) {
Vec3 subpos = pos + disp * Vec3(0.5+di, 0.5+dj, 0.5+dk);
subpos += jlen * (Vec3(1,1,1) - 2.0 * mRand.getVec3());
if(!is3D) subpos[2] = 0.5;
if(exclude && exclude->getInterpolated(subpos) <= 0.) continue;
if(!shape.isInside(subpos)) continue;
parts.addBuffered(subpos);
}
}
parts.insertBufferedParticles();
} 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, "sampleShapeWithParticles" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const Shape& shape = *_args.getPtr<Shape >("shape",0,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",1,&_lock); BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",2,&_lock); const int discretization = _args.get<int >("discretization",3,&_lock); const Real randomness = _args.get<Real >("randomness",4,&_lock); const bool reset = _args.getOpt<bool >("reset",5,false,&_lock); const bool refillEmpty = _args.getOpt<bool >("refillEmpty",6,false,&_lock); const LevelsetGrid* exclude = _args.getPtrOpt<LevelsetGrid >("exclude",7,NULL,&_lock); _retval = getPyNone(); sampleShapeWithParticles(shape,flags,parts,discretization,randomness,reset,refillEmpty,exclude); _args.check(); }
pbFinalizePlugin(parent,"sampleShapeWithParticles", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("sampleShapeWithParticles",e.what()); return 0; }
}
static const Pb::Register _RP_sampleShapeWithParticles ("","sampleShapeWithParticles",_W_2);
extern "C" {
void PbRegister_sampleShapeWithParticles() {
KEEP_UNUSED(_RP_sampleShapeWithParticles); }
}
//! mark fluid cells and helpers
struct knClearFluidFlags : public KernelBase {
knClearFluidFlags(FlagGrid& flags, int dummy=0) : KernelBase(&flags,0) ,flags(flags),dummy(dummy) {
runMessage(); run(); }
inline void op(int i, int j, int k, FlagGrid& flags, int dummy=0 ) {
if (flags.isFluid(i,j,k)) {
flags(i,j,k) = (flags(i,j,k) | FlagGrid::TypeEmpty) & ~FlagGrid::TypeFluid;
}
} inline FlagGrid& getArg0() {
return flags; }
typedef FlagGrid type0;inline int& getArg1() {
return dummy; }
typedef int type1; void runMessage() { debMsg("Executing kernel knClearFluidFlags ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,flags,dummy); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,flags,dummy); }
}
}
FlagGrid& flags; int dummy; }
;
struct knSetNbObstacle : public KernelBase {
knSetNbObstacle(FlagGrid& nflags, const FlagGrid& flags, const Grid<Real>& phiObs) : KernelBase(&nflags,1) ,nflags(nflags),flags(flags),phiObs(phiObs) {
runMessage(); run(); }
inline void op(int i, int j, int k, FlagGrid& nflags, const FlagGrid& flags, const Grid<Real>& phiObs ) {
if ( phiObs(i,j,k)>0. ) return;
if (flags.isEmpty(i,j,k)) {
bool set=false;
if( (flags.isFluid(i-1,j,k)) && (phiObs(i+1,j,k)<=0.) ) set=true;
if( (flags.isFluid(i+1,j,k)) && (phiObs(i-1,j,k)<=0.) ) set=true;
if( (flags.isFluid(i,j-1,k)) && (phiObs(i,j+1,k)<=0.) ) set=true;
if( (flags.isFluid(i,j+1,k)) && (phiObs(i,j-1,k)<=0.) ) set=true;
if(flags.is3D()) {
if( (flags.isFluid(i,j,k-1)) && (phiObs(i,j,k+1)<=0.) ) set=true;
if( (flags.isFluid(i,j,k+1)) && (phiObs(i,j,k-1)<=0.) ) set=true;
}
if(set) nflags(i,j,k) = (flags(i,j,k) | FlagGrid::TypeFluid) & ~FlagGrid::TypeEmpty;
}
} inline FlagGrid& getArg0() {
return nflags; }
typedef FlagGrid type0;inline const FlagGrid& getArg1() {
return flags; }
typedef FlagGrid type1;inline const Grid<Real>& getArg2() {
return phiObs; }
typedef Grid<Real> type2; void runMessage() { debMsg("Executing kernel knSetNbObstacle ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,nflags,flags,phiObs); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,nflags,flags,phiObs); }
}
}
FlagGrid& nflags; const FlagGrid& flags; const Grid<Real>& phiObs; }
;
void markFluidCells(const BasicParticleSystem& parts, FlagGrid& flags, const Grid<Real>* phiObs=NULL, const ParticleDataImpl<int>* ptype=NULL, const int exclude=0) {
// remove all fluid cells
knClearFluidFlags(flags, 0);
// mark all particles in flaggrid as fluid
for(IndexInt idx=0; idx<parts.size(); idx++) {
if (!parts.isActive(idx) || (ptype && ((*ptype)[idx] & exclude))) continue;
Vec3i p = toVec3i( parts.getPos(idx) );
if (flags.isInBounds(p) && flags.isEmpty(p))
flags(p) = (flags(p) | FlagGrid::TypeFluid) & ~FlagGrid::TypeEmpty;
}
// special for second order obstacle BCs, check empty cells in boundary region
if(phiObs) {
FlagGrid tmp(flags);
knSetNbObstacle(tmp, flags, *phiObs);
flags.swap(tmp);
}
} 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, "markFluidCells" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",0,&_lock); FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",1,&_lock); const Grid<Real>* phiObs = _args.getPtrOpt<Grid<Real> >("phiObs",2,NULL,&_lock); const ParticleDataImpl<int>* ptype = _args.getPtrOpt<ParticleDataImpl<int> >("ptype",3,NULL,&_lock); const int exclude = _args.getOpt<int >("exclude",4,0,&_lock); _retval = getPyNone(); markFluidCells(parts,flags,phiObs,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"markFluidCells", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("markFluidCells",e.what()); return 0; }
}
static const Pb::Register _RP_markFluidCells ("","markFluidCells",_W_3);
extern "C" {
void PbRegister_markFluidCells() {
KEEP_UNUSED(_RP_markFluidCells); }
}
// for testing purposes only...
void testInitGridWithPos(Grid<Real>& grid) {
FOR_IJK(grid) { grid(i,j,k) = norm( Vec3(i,j,k) ); }
} 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, "testInitGridWithPos" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; Grid<Real>& grid = *_args.getPtr<Grid<Real> >("grid",0,&_lock); _retval = getPyNone(); testInitGridWithPos(grid); _args.check(); }
pbFinalizePlugin(parent,"testInitGridWithPos", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("testInitGridWithPos",e.what()); return 0; }
}
static const Pb::Register _RP_testInitGridWithPos ("","testInitGridWithPos",_W_4);
extern "C" {
void PbRegister_testInitGridWithPos() {
KEEP_UNUSED(_RP_testInitGridWithPos); }
}
//! helper to calculate particle radius factor to cover the diagonal of a cell in 2d/3d
inline Real calculateRadiusFactor(const Grid<Real>& grid, Real factor) {
return (grid.is3D() ? sqrt(3.) : sqrt(2.) ) * (factor+.01); // note, a 1% safety factor is added here
}
//! re-sample particles based on an input levelset
// optionally skip seeding new particles in "exclude" SDF
void adjustNumber(BasicParticleSystem& parts, const MACGrid& vel, const FlagGrid& flags, int minParticles, int maxParticles, const LevelsetGrid& phi, Real radiusFactor=1. , Real narrowBand=-1. , const Grid<Real>* exclude=NULL ) {
// which levelset to use as threshold
const Real SURFACE_LS = -1.0 * calculateRadiusFactor(phi, radiusFactor);
Grid<int> tmp( vel.getParent() );
std::ostringstream out;
// count particles in cells, and delete excess particles
for (IndexInt idx=0; idx<(int)parts.size(); idx++) {
if (parts.isActive(idx)) {
Vec3i p = toVec3i( parts.getPos(idx) );
if (!tmp.isInBounds(p) ) {
parts.kill(idx); // out of domain, remove
continue;
}
Real phiv = phi.getInterpolated( parts.getPos(idx) );
if( phiv > 0 ) { parts.kill(idx); continue; }
if( narrowBand>0. && phiv < -narrowBand) { parts.kill(idx); continue; }
bool atSurface = false;
if (phiv > SURFACE_LS) atSurface = true;
int num = tmp(p);
// dont delete particles in non fluid cells here, the particles are "always right"
if ( num > maxParticles && (!atSurface) ) {
parts.kill(idx);
} else {
tmp(p) = num+1;
}
}
}
// seed new particles
RandomStream mRand(9832);
FOR_IJK(tmp) {
int cnt = tmp(i,j,k);
// skip cells near surface
if (phi(i,j,k) > SURFACE_LS) continue;
if( narrowBand>0. && phi(i,j,k) < -narrowBand ) { continue; }
if( exclude && ( (*exclude)(i,j,k) < 0.) ) { continue; }
if (flags.isFluid(i,j,k) && cnt < minParticles) {
for (int m=cnt; m < minParticles; m++) {
Vec3 pos = Vec3(i,j,k) + mRand.getVec3();
//Vec3 pos (i + 0.5, j + 0.5, k + 0.5); // cell center
parts.addBuffered( pos );
}
}
}
parts.doCompress();
parts.insertBufferedParticles();
} 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, "adjustNumber" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",0,&_lock); const MACGrid& vel = *_args.getPtr<MACGrid >("vel",1,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",2,&_lock); int minParticles = _args.get<int >("minParticles",3,&_lock); int maxParticles = _args.get<int >("maxParticles",4,&_lock); const LevelsetGrid& phi = *_args.getPtr<LevelsetGrid >("phi",5,&_lock); Real radiusFactor = _args.getOpt<Real >("radiusFactor",6,1. ,&_lock); Real narrowBand = _args.getOpt<Real >("narrowBand",7,-1. ,&_lock); const Grid<Real>* exclude = _args.getPtrOpt<Grid<Real> >("exclude",8,NULL ,&_lock); _retval = getPyNone(); adjustNumber(parts,vel,flags,minParticles,maxParticles,phi,radiusFactor,narrowBand,exclude); _args.check(); }
pbFinalizePlugin(parent,"adjustNumber", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("adjustNumber",e.what()); return 0; }
}
static const Pb::Register _RP_adjustNumber ("","adjustNumber",_W_5);
extern "C" {
void PbRegister_adjustNumber() {
KEEP_UNUSED(_RP_adjustNumber); }
}
// simple and slow helper conversion to show contents of int grids like a real grid in the ui
// (use eg to quickly display contents of the particle-index grid)
void debugIntToReal(const Grid<int>& source, Grid<Real>& dest, Real factor=1. ) {
FOR_IJK( source ) { dest(i,j,k) = (Real)source(i,j,k) * factor; }
} 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, "debugIntToReal" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const Grid<int>& source = *_args.getPtr<Grid<int> >("source",0,&_lock); Grid<Real>& dest = *_args.getPtr<Grid<Real> >("dest",1,&_lock); Real factor = _args.getOpt<Real >("factor",2,1. ,&_lock); _retval = getPyNone(); debugIntToReal(source,dest,factor); _args.check(); }
pbFinalizePlugin(parent,"debugIntToReal", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("debugIntToReal",e.what()); return 0; }
}
static const Pb::Register _RP_debugIntToReal ("","debugIntToReal",_W_6);
extern "C" {
void PbRegister_debugIntToReal() {
KEEP_UNUSED(_RP_debugIntToReal); }
}
// build a grid that contains indices for a particle system
// the particles in a cell i,j,k are particles[index(i,j,k)] to particles[index(i+1,j,k)-1]
// (ie, particles[index(i+1,j,k)] already belongs to cell i+1,j,k)
void gridParticleIndex(const BasicParticleSystem& parts, ParticleIndexSystem& indexSys, const FlagGrid& flags, Grid<int>& index, Grid<int>* counter=NULL ) {
bool delCounter = false;
if(!counter) { counter = new Grid<int>( flags.getParent() ); delCounter=true; }
else { counter->clear(); }
// count particles in cells, and delete excess particles
index.clear();
int inactive = 0;
for (IndexInt idx=0; idx<(IndexInt)parts.size(); idx++) {
if (parts.isActive(idx)) {
// check index for validity...
Vec3i p = toVec3i( parts.getPos(idx) );
if (! index.isInBounds(p)) { inactive++; continue; }
index(p)++;
} else {
inactive++;
}
}
// note - this one might be smaller...
indexSys.resize( parts.size()-inactive );
// convert per cell number to continuous index
IndexInt idx=0;
FOR_IJK( index ) {
int num = index(i,j,k);
index(i,j,k) = idx;
idx += num;
}
// add particles to indexed array, we still need a per cell particle counter
for (IndexInt idx=0; idx<(IndexInt)parts.size(); idx++) {
if (!parts.isActive(idx)) continue;
Vec3i p = toVec3i( parts.getPos(idx) );
if (! index.isInBounds(p)) { continue; }
// initialize position and index into original array
//indexSys[ index(p)+(*counter)(p) ].pos = parts[idx].pos;
indexSys[ index(p)+(*counter)(p) ].sourceIndex = idx;
(*counter)(p)++;
}
if(delCounter) delete counter;
} 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, "gridParticleIndex" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",0,&_lock); ParticleIndexSystem& indexSys = *_args.getPtr<ParticleIndexSystem >("indexSys",1,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",2,&_lock); Grid<int>& index = *_args.getPtr<Grid<int> >("index",3,&_lock); Grid<int>* counter = _args.getPtrOpt<Grid<int> >("counter",4,NULL ,&_lock); _retval = getPyNone(); gridParticleIndex(parts,indexSys,flags,index,counter); _args.check(); }
pbFinalizePlugin(parent,"gridParticleIndex", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("gridParticleIndex",e.what()); return 0; }
}
static const Pb::Register _RP_gridParticleIndex ("","gridParticleIndex",_W_7);
extern "C" {
void PbRegister_gridParticleIndex() {
KEEP_UNUSED(_RP_gridParticleIndex); }
}
struct ComputeUnionLevelsetPindex : public KernelBase {
ComputeUnionLevelsetPindex(const Grid<int>& index, const BasicParticleSystem& parts, const ParticleIndexSystem& indexSys, LevelsetGrid& phi, const Real radius, const ParticleDataImpl<int> *ptype, const int exclude) : KernelBase(&index,0) ,index(index),parts(parts),indexSys(indexSys),phi(phi),radius(radius),ptype(ptype),exclude(exclude) {
runMessage(); run(); }
inline void op(int i, int j, int k, const Grid<int>& index, const BasicParticleSystem& parts, const ParticleIndexSystem& indexSys, LevelsetGrid& phi, const Real radius, const ParticleDataImpl<int> *ptype, const int exclude ) {
const Vec3 gridPos = Vec3(i,j,k) + Vec3(0.5); // shifted by half cell
Real phiv = radius * 1.0; // outside
int r = int(radius) + 1;
int rZ = phi.is3D() ? r : 0;
for(int zj=k-rZ; zj<=k+rZ; zj++)
for(int yj=j-r ; yj<=j+r ; yj++)
for(int xj=i-r ; xj<=i+r ; xj++) {
if (!phi.isInBounds(Vec3i(xj,yj,zj))) continue;
// note, for the particle indices in indexSys the access is periodic (ie, dont skip for eg inBounds(sx,10,10)
IndexInt isysIdxS = index.index(xj,yj,zj);
IndexInt pStart = index(isysIdxS), pEnd=0;
if(phi.isInBounds(isysIdxS+1)) pEnd = index(isysIdxS+1);
else pEnd = indexSys.size();
// now loop over particles in cell
for(IndexInt p=pStart; p<pEnd; ++p) {
const int psrc = indexSys[p].sourceIndex;
if(ptype && ((*ptype)[psrc] & exclude)) continue;
const Vec3 pos = parts[psrc].pos;
phiv = std::min( phiv , fabs( norm(gridPos-pos) )-radius );
}
}
phi(i,j,k) = phiv;
} inline const Grid<int>& getArg0() {
return index; }
typedef Grid<int> type0;inline const BasicParticleSystem& getArg1() {
return parts; }
typedef BasicParticleSystem type1;inline const ParticleIndexSystem& getArg2() {
return indexSys; }
typedef ParticleIndexSystem type2;inline LevelsetGrid& getArg3() {
return phi; }
typedef LevelsetGrid type3;inline const Real& getArg4() {
return radius; }
typedef Real type4;inline const ParticleDataImpl<int> * getArg5() {
return ptype; }
typedef ParticleDataImpl<int> type5;inline const int& getArg6() {
return exclude; }
typedef int type6; void runMessage() { debMsg("Executing kernel ComputeUnionLevelsetPindex ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,index,parts,indexSys,phi,radius,ptype,exclude); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,index,parts,indexSys,phi,radius,ptype,exclude); }
}
}
const Grid<int>& index; const BasicParticleSystem& parts; const ParticleIndexSystem& indexSys; LevelsetGrid& phi; const Real radius; const ParticleDataImpl<int> * ptype; const int exclude; }
;
void unionParticleLevelset(const BasicParticleSystem& parts, const ParticleIndexSystem& indexSys, const FlagGrid& flags, const Grid<int>& index, LevelsetGrid& phi, const Real radiusFactor=1., const ParticleDataImpl<int> *ptype=NULL, const int exclude=0) {
// use half a cell diagonal as base radius
const Real radius = 0.5 * calculateRadiusFactor(phi, radiusFactor);
// no reset of phi necessary here
ComputeUnionLevelsetPindex(index, parts, indexSys, phi, radius, ptype, exclude);
phi.setBound(0.5, 0);
} static PyObject* _W_8 (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, "unionParticleLevelset" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",0,&_lock); const ParticleIndexSystem& indexSys = *_args.getPtr<ParticleIndexSystem >("indexSys",1,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",2,&_lock); const Grid<int>& index = *_args.getPtr<Grid<int> >("index",3,&_lock); LevelsetGrid& phi = *_args.getPtr<LevelsetGrid >("phi",4,&_lock); const Real radiusFactor = _args.getOpt<Real >("radiusFactor",5,1.,&_lock); const ParticleDataImpl<int> * ptype = _args.getPtrOpt<ParticleDataImpl<int> >("ptype",6,NULL,&_lock); const int exclude = _args.getOpt<int >("exclude",7,0,&_lock); _retval = getPyNone(); unionParticleLevelset(parts,indexSys,flags,index,phi,radiusFactor,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"unionParticleLevelset", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("unionParticleLevelset",e.what()); return 0; }
}
static const Pb::Register _RP_unionParticleLevelset ("","unionParticleLevelset",_W_8);
extern "C" {
void PbRegister_unionParticleLevelset() {
KEEP_UNUSED(_RP_unionParticleLevelset); }
}
//! kernel for computing averaged particle level set weights
struct ComputeAveragedLevelsetWeight : public KernelBase {
ComputeAveragedLevelsetWeight(const BasicParticleSystem& parts, const Grid<int>& index, const ParticleIndexSystem& indexSys, LevelsetGrid& phi, const Real radius, const ParticleDataImpl<int>* ptype, const int exclude, Grid<Vec3>* save_pAcc = NULL, Grid<Real>* save_rAcc = NULL) : KernelBase(&index,0) ,parts(parts),index(index),indexSys(indexSys),phi(phi),radius(radius),ptype(ptype),exclude(exclude),save_pAcc(save_pAcc),save_rAcc(save_rAcc) {
runMessage(); run(); }
inline void op(int i, int j, int k, const BasicParticleSystem& parts, const Grid<int>& index, const ParticleIndexSystem& indexSys, LevelsetGrid& phi, const Real radius, const ParticleDataImpl<int>* ptype, const int exclude, Grid<Vec3>* save_pAcc = NULL, Grid<Real>* save_rAcc = NULL ) {
const Vec3 gridPos = Vec3(i,j,k) + Vec3(0.5); // shifted by half cell
Real phiv = radius * 1.0; // outside
// loop over neighborhood, similar to ComputeUnionLevelsetPindex
const Real sradiusInv = 1. / (4. * radius * radius) ;
int r = int(1. * radius) + 1;
int rZ = phi.is3D() ? r : 0;
// accumulators
Real wacc = 0.;
Vec3 pacc = Vec3(0.);
Real racc = 0.;
for(int zj=k-rZ; zj<=k+rZ; zj++)
for(int yj=j-r ; yj<=j+r ; yj++)
for(int xj=i-r ; xj<=i+r ; xj++) {
if (! phi.isInBounds(Vec3i(xj,yj,zj)) ) continue;
IndexInt isysIdxS = index.index(xj,yj,zj);
IndexInt pStart = index(isysIdxS), pEnd=0;
if(phi.isInBounds(isysIdxS+1)) pEnd = index(isysIdxS+1);
else pEnd = indexSys.size();
for(IndexInt p=pStart; p<pEnd; ++p) {
IndexInt psrc = indexSys[p].sourceIndex;
if(ptype && ((*ptype)[psrc] & exclude)) continue;
Vec3 pos = parts[psrc].pos;
Real s = normSquare(gridPos-pos) * sradiusInv;
//Real w = std::max(0., cubed(1.-s) );
Real w = std::max(0., (1.-s)); // a bit smoother
wacc += w;
racc += radius * w;
pacc += pos * w;
}
}
if(wacc > VECTOR_EPSILON) {
racc /= wacc;
pacc /= wacc;
phiv = fabs( norm(gridPos-pacc) )-racc;
if (save_pAcc) (*save_pAcc)(i, j, k) = pacc;
if (save_rAcc) (*save_rAcc)(i, j, k) = racc;
}
phi(i,j,k) = phiv;
} inline const BasicParticleSystem& getArg0() {
return parts; }
typedef BasicParticleSystem type0;inline const Grid<int>& getArg1() {
return index; }
typedef Grid<int> type1;inline const ParticleIndexSystem& getArg2() {
return indexSys; }
typedef ParticleIndexSystem type2;inline LevelsetGrid& getArg3() {
return phi; }
typedef LevelsetGrid type3;inline const Real& getArg4() {
return radius; }
typedef Real type4;inline const ParticleDataImpl<int>* getArg5() {
return ptype; }
typedef ParticleDataImpl<int> type5;inline const int& getArg6() {
return exclude; }
typedef int type6;inline Grid<Vec3>* getArg7() {
return save_pAcc; }
typedef Grid<Vec3> type7;inline Grid<Real>* getArg8() {
return save_rAcc; }
typedef Grid<Real> type8; void runMessage() { debMsg("Executing kernel ComputeAveragedLevelsetWeight ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,parts,index,indexSys,phi,radius,ptype,exclude,save_pAcc,save_rAcc); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,parts,index,indexSys,phi,radius,ptype,exclude,save_pAcc,save_rAcc); }
}
}
const BasicParticleSystem& parts; const Grid<int>& index; const ParticleIndexSystem& indexSys; LevelsetGrid& phi; const Real radius; const ParticleDataImpl<int>* ptype; const int exclude; Grid<Vec3>* save_pAcc; Grid<Real>* save_rAcc; }
;
template<class T> T smoothingValue(const Grid<T> val, int i, int j, int k, T center) {
return val(i,j,k);
}
// smoothing, and
template <class T> struct knSmoothGrid : public KernelBase {
knSmoothGrid(const Grid<T>& me, Grid<T>& tmp, Real factor) : KernelBase(&me,1) ,me(me),tmp(tmp),factor(factor) {
runMessage(); run(); }
inline void op(int i, int j, int k, const Grid<T>& me, Grid<T>& tmp, Real factor ) {
T val = me(i,j,k) +
me(i+1,j,k) + me(i-1,j,k) +
me(i,j+1,k) + me(i,j-1,k) ;
if(me.is3D()) {
val += me(i,j,k+1) + me(i,j,k-1);
}
tmp(i,j,k) = val * factor;
} inline const Grid<T>& getArg0() {
return me; }
typedef Grid<T> type0;inline Grid<T>& getArg1() {
return tmp; }
typedef Grid<T> type1;inline Real& getArg2() {
return factor; }
typedef Real type2; void runMessage() { debMsg("Executing kernel knSmoothGrid ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,me,tmp,factor); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,me,tmp,factor); }
}
}
const Grid<T>& me; Grid<T>& tmp; Real factor; }
;
template <class T> struct knSmoothGridNeg : public KernelBase {
knSmoothGridNeg(const Grid<T>& me, Grid<T>& tmp, Real factor) : KernelBase(&me,1) ,me(me),tmp(tmp),factor(factor) {
runMessage(); run(); }
inline void op(int i, int j, int k, const Grid<T>& me, Grid<T>& tmp, Real factor ) {
T val = me(i,j,k) +
me(i+1,j,k) + me(i-1,j,k) +
me(i,j+1,k) + me(i,j-1,k) ;
if(me.is3D()) {
val += me(i,j,k+1) + me(i,j,k-1);
}
val *= factor;
if(val<tmp(i,j,k)) tmp(i,j,k) = val;
else tmp(i,j,k) = me(i,j,k);
} inline const Grid<T>& getArg0() {
return me; }
typedef Grid<T> type0;inline Grid<T>& getArg1() {
return tmp; }
typedef Grid<T> type1;inline Real& getArg2() {
return factor; }
typedef Real type2; void runMessage() { debMsg("Executing kernel knSmoothGridNeg ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,me,tmp,factor); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,me,tmp,factor); }
}
}
const Grid<T>& me; Grid<T>& tmp; Real factor; }
;
//! Zhu & Bridson particle level set creation
void averagedParticleLevelset(const BasicParticleSystem& parts, const ParticleIndexSystem& indexSys, const FlagGrid& flags, const Grid<int>& index, LevelsetGrid& phi, const Real radiusFactor=1., const int smoothen=1, const int smoothenNeg=1, const ParticleDataImpl<int>* ptype=NULL, const int exclude=0) {
// use half a cell diagonal as base radius
const Real radius = 0.5 * calculateRadiusFactor(phi, radiusFactor);
ComputeAveragedLevelsetWeight(parts, index, indexSys, phi, radius, ptype, exclude);
// post-process level-set
for(int i=0; i<std::max(smoothen,smoothenNeg); ++i) {
LevelsetGrid tmp(flags.getParent());
if(i<smoothen) {
knSmoothGrid <Real> (phi,tmp, 1./(phi.is3D() ? 7. : 5.) );
phi.swap(tmp);
}
if(i<smoothenNeg) {
knSmoothGridNeg <Real> (phi,tmp, 1./(phi.is3D() ? 7. : 5.) );
phi.swap(tmp);
}
}
phi.setBound(0.5, 0);
} static PyObject* _W_9 (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, "averagedParticleLevelset" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",0,&_lock); const ParticleIndexSystem& indexSys = *_args.getPtr<ParticleIndexSystem >("indexSys",1,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",2,&_lock); const Grid<int>& index = *_args.getPtr<Grid<int> >("index",3,&_lock); LevelsetGrid& phi = *_args.getPtr<LevelsetGrid >("phi",4,&_lock); const Real radiusFactor = _args.getOpt<Real >("radiusFactor",5,1.,&_lock); const int smoothen = _args.getOpt<int >("smoothen",6,1,&_lock); const int smoothenNeg = _args.getOpt<int >("smoothenNeg",7,1,&_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(); averagedParticleLevelset(parts,indexSys,flags,index,phi,radiusFactor,smoothen,smoothenNeg,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"averagedParticleLevelset", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("averagedParticleLevelset",e.what()); return 0; }
}
static const Pb::Register _RP_averagedParticleLevelset ("","averagedParticleLevelset",_W_9);
extern "C" {
void PbRegister_averagedParticleLevelset() {
KEEP_UNUSED(_RP_averagedParticleLevelset); }
}
//! kernel for improvedParticleLevelset
struct correctLevelset : public KernelBase {
correctLevelset(LevelsetGrid& phi, const Grid<Vec3>& pAcc, const Grid<Real>& rAcc, const Real radius, const Real t_low, const Real t_high) : KernelBase(&phi,1) ,phi(phi),pAcc(pAcc),rAcc(rAcc),radius(radius),t_low(t_low),t_high(t_high) {
runMessage(); run(); }
inline void op(int i, int j, int k, LevelsetGrid& phi, const Grid<Vec3>& pAcc, const Grid<Real>& rAcc, const Real radius, const Real t_low, const Real t_high ) {
if (rAcc(i, j, k) <= VECTOR_EPSILON) return; //outside nothing happens
Real x = pAcc(i, j, k).x;
// create jacobian of pAcc via central differences
Matrix3x3f jacobian = Matrix3x3f(
0.5 * (pAcc(i+1, j, k ).x - pAcc(i-1, j, k ).x),
0.5 * (pAcc(i, j+1, k ).x - pAcc(i, j-1, k ).x),
0.5 * (pAcc(i, j , k+1).x - pAcc(i, j, k-1).x),
0.5 * (pAcc(i+1, j, k ).y - pAcc(i-1, j, k ).y),
0.5 * (pAcc(i, j+1, k ).y - pAcc(i, j-1, k ).y),
0.5 * (pAcc(i, j, k+1).y - pAcc(i, j, k-1).y),
0.5 * (pAcc(i+1, j, k ).z - pAcc(i-1, j, k ).z),
0.5 * (pAcc(i, j+1, k ).z - pAcc(i, j-1, k ).z),
0.5 * (pAcc(i, j, k+1).z - pAcc(i, j, k-1).z)
);
// compute largest eigenvalue of jacobian
Vec3 EV = jacobian.eigenvalues();
Real maxEV = std::max(std::max(EV.x, EV.y), EV.z);
// calculate correction factor
Real correction = 1;
if (maxEV >= t_low) {
Real t = (t_high - maxEV) / (t_high - t_low);
correction = t*t*t - 3 * t*t + 3 * t;
}
correction = (correction < 0) ? 0 : correction; // enforce correction factor to [0,1] (not explicitly in paper)
const Vec3 gridPos = Vec3(i, j, k) + Vec3(0.5); // shifted by half cell
const Real correctedPhi = fabs(norm(gridPos - pAcc(i, j, k))) - rAcc(i, j, k) * correction;
phi(i, j, k) = (correctedPhi > radius) ? radius : correctedPhi; // adjust too high outside values when too few particles are
// nearby to make smoothing possible (not in paper)
} inline LevelsetGrid& getArg0() {
return phi; }
typedef LevelsetGrid type0;inline const Grid<Vec3>& getArg1() {
return pAcc; }
typedef Grid<Vec3> type1;inline const Grid<Real>& getArg2() {
return rAcc; }
typedef Grid<Real> type2;inline const Real& getArg3() {
return radius; }
typedef Real type3;inline const Real& getArg4() {
return t_low; }
typedef Real type4;inline const Real& getArg5() {
return t_high; }
typedef Real type5; void runMessage() { debMsg("Executing kernel correctLevelset ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,phi,pAcc,rAcc,radius,t_low,t_high); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,phi,pAcc,rAcc,radius,t_low,t_high); }
}
}
LevelsetGrid& phi; const Grid<Vec3>& pAcc; const Grid<Real>& rAcc; const Real radius; const Real t_low; const Real t_high; }
;
//! Approach from "A unified particle model for fluid-solid interactions" by Solenthaler et al. in 2007
void improvedParticleLevelset(const BasicParticleSystem& parts, const ParticleIndexSystem& indexSys, const FlagGrid& flags, const Grid<int>& index, LevelsetGrid& phi, const Real radiusFactor = 1., const int smoothen = 1,const int smoothenNeg = 1, const Real t_low = 0.4, const Real t_high = 3.5, const ParticleDataImpl<int>* ptype = NULL, const int exclude = 0) {
// create temporary grids to store values from levelset weight computation
Grid<Vec3> save_pAcc(flags.getParent());
Grid<Real> save_rAcc(flags.getParent());
const Real radius = 0.5 * calculateRadiusFactor(phi, radiusFactor); // use half a cell diagonal as base radius
ComputeAveragedLevelsetWeight(parts, index, indexSys, phi, radius, ptype, exclude, &save_pAcc, &save_rAcc);
correctLevelset(phi, save_pAcc, save_rAcc, radius, t_low, t_high);
// post-process level-set
for (int i = 0; i<std::max(smoothen, smoothenNeg); ++i) {
LevelsetGrid tmp(flags.getParent());
if (i<smoothen) {
knSmoothGrid <Real>(phi, tmp, 1. / (phi.is3D() ? 7. : 5.));
phi.swap(tmp);
}
if (i<smoothenNeg) {
knSmoothGridNeg <Real>(phi, tmp, 1. / (phi.is3D() ? 7. : 5.));
phi.swap(tmp);
}
}
phi.setBound(0.5, 0);
} static PyObject* _W_10 (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, "improvedParticleLevelset" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",0,&_lock); const ParticleIndexSystem& indexSys = *_args.getPtr<ParticleIndexSystem >("indexSys",1,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",2,&_lock); const Grid<int>& index = *_args.getPtr<Grid<int> >("index",3,&_lock); LevelsetGrid& phi = *_args.getPtr<LevelsetGrid >("phi",4,&_lock); const Real radiusFactor = _args.getOpt<Real >("radiusFactor",5,1.,&_lock); const int smoothen = _args.getOpt<int >("smoothen",6,1,&_lock); const int smoothenNeg = _args.getOpt<int >("smoothenNeg",7,1,&_lock); const Real t_low = _args.getOpt<Real >("t_low",8,0.4,&_lock); const Real t_high = _args.getOpt<Real >("t_high",9,3.5,&_lock); const ParticleDataImpl<int>* ptype = _args.getPtrOpt<ParticleDataImpl<int> >("ptype",10,NULL,&_lock); const int exclude = _args.getOpt<int >("exclude",11,0,&_lock); _retval = getPyNone(); improvedParticleLevelset(parts,indexSys,flags,index,phi,radiusFactor,smoothen,smoothenNeg,t_low,t_high,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"improvedParticleLevelset", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("improvedParticleLevelset",e.what()); return 0; }
}
static const Pb::Register _RP_improvedParticleLevelset ("","improvedParticleLevelset",_W_10);
extern "C" {
void PbRegister_improvedParticleLevelset() {
KEEP_UNUSED(_RP_improvedParticleLevelset); }
}
struct knPushOutofObs : public KernelBase {
knPushOutofObs(BasicParticleSystem& parts, const FlagGrid& flags, const Grid<Real>& phiObs, const Real shift, const Real thresh, const ParticleDataImpl<int>* ptype, const int exclude) : KernelBase(parts.size()) ,parts(parts),flags(flags),phiObs(phiObs),shift(shift),thresh(thresh),ptype(ptype),exclude(exclude) {
runMessage(); run(); }
inline void op(IndexInt idx, BasicParticleSystem& parts, const FlagGrid& flags, const Grid<Real>& phiObs, const Real shift, const Real thresh, const ParticleDataImpl<int>* ptype, const int exclude ) {
if (!parts.isActive(idx) || (ptype && ((*ptype)[idx] & exclude))) return;
Vec3i p = toVec3i( parts.getPos(idx) );
if (!flags.isInBounds(p)) return;
Real v = phiObs.getInterpolated(parts.getPos(idx));
if(v < thresh) {
Vec3 grad = getGradient( phiObs, p.x,p.y,p.z );
if( normalize(grad) < VECTOR_EPSILON ) return;
parts.setPos(idx, parts.getPos(idx) + grad*(thresh - v + shift));
}
} inline BasicParticleSystem& getArg0() {
return parts; }
typedef BasicParticleSystem type0;inline const FlagGrid& getArg1() {
return flags; }
typedef FlagGrid type1;inline const Grid<Real>& getArg2() {
return phiObs; }
typedef Grid<Real> type2;inline const Real& getArg3() {
return shift; }
typedef Real type3;inline const Real& getArg4() {
return thresh; }
typedef Real type4;inline const ParticleDataImpl<int>* getArg5() {
return ptype; }
typedef ParticleDataImpl<int> type5;inline const int& getArg6() {
return exclude; }
typedef int type6; void runMessage() { debMsg("Executing kernel knPushOutofObs ", 3); debMsg("Kernel range" << " size "<< size << " " , 4); }; void run() {
const IndexInt _sz = size;
#pragma omp parallel
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,parts,flags,phiObs,shift,thresh,ptype,exclude); }
}
BasicParticleSystem& parts; const FlagGrid& flags; const Grid<Real>& phiObs; const Real shift; const Real thresh; const ParticleDataImpl<int>* ptype; const int exclude; }
;
//! push particles out of obstacle levelset
void pushOutofObs(BasicParticleSystem& parts, const FlagGrid& flags, const Grid<Real>& phiObs, const Real shift=0, const Real thresh=0, const ParticleDataImpl<int>* ptype=NULL, const int exclude=0) {
knPushOutofObs(parts, flags, phiObs, shift, thresh, ptype, exclude);
} static PyObject* _W_11 (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, "pushOutofObs" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",0,&_lock); const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",1,&_lock); const Grid<Real>& phiObs = *_args.getPtr<Grid<Real> >("phiObs",2,&_lock); const Real shift = _args.getOpt<Real >("shift",3,0,&_lock); const Real thresh = _args.getOpt<Real >("thresh",4,0,&_lock); const ParticleDataImpl<int>* ptype = _args.getPtrOpt<ParticleDataImpl<int> >("ptype",5,NULL,&_lock); const int exclude = _args.getOpt<int >("exclude",6,0,&_lock); _retval = getPyNone(); pushOutofObs(parts,flags,phiObs,shift,thresh,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"pushOutofObs", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("pushOutofObs",e.what()); return 0; }
}
static const Pb::Register _RP_pushOutofObs ("","pushOutofObs",_W_11);
extern "C" {
void PbRegister_pushOutofObs() {
KEEP_UNUSED(_RP_pushOutofObs); }
}
//******************************************************************************
// grid interpolation functions
template <class T> struct knSafeDivReal : public KernelBase {
knSafeDivReal(Grid<T>& me, const Grid<Real>& other, Real cutoff=VECTOR_EPSILON) : KernelBase(&me,0) ,me(me),other(other),cutoff(cutoff) {
runMessage(); run(); }
inline void op(IndexInt idx, Grid<T>& me, const Grid<Real>& other, Real cutoff=VECTOR_EPSILON ) {
if(other[idx]<cutoff) {
me[idx] = 0.;
} else {
T div( other[idx] );
me[idx] = safeDivide(me[idx], div );
}
} inline Grid<T>& getArg0() {
return me; }
typedef Grid<T> type0;inline const Grid<Real>& getArg1() {
return other; }
typedef Grid<Real> type1;inline Real& getArg2() {
return cutoff; }
typedef Real type2; void runMessage() { debMsg("Executing kernel knSafeDivReal ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const IndexInt _sz = size;
#pragma omp parallel
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,me,other,cutoff); }
}
Grid<T>& me; const Grid<Real>& other; Real cutoff; }
;
// Set velocities on the grid from the particle system
struct knMapLinearVec3ToMACGrid : public KernelBase {
knMapLinearVec3ToMACGrid(const BasicParticleSystem& p, const FlagGrid& flags, const MACGrid& vel, Grid<Vec3>& tmp, const ParticleDataImpl<Vec3>& pvel, const ParticleDataImpl<int>* ptype, const int exclude) : KernelBase(p.size()) ,p(p),flags(flags),vel(vel),tmp(tmp),pvel(pvel),ptype(ptype),exclude(exclude) {
runMessage(); run(); }
inline void op(IndexInt idx, const BasicParticleSystem& p, const FlagGrid& flags, const MACGrid& vel, Grid<Vec3>& tmp, const ParticleDataImpl<Vec3>& pvel, const ParticleDataImpl<int>* ptype, const int exclude ) {
unusedParameter(flags);
if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude))) return;
vel.setInterpolated( p[idx].pos, pvel[idx], &tmp[0] );
} inline const BasicParticleSystem& getArg0() {
return p; }
typedef BasicParticleSystem type0;inline const FlagGrid& getArg1() {
return flags; }
typedef FlagGrid type1;inline const MACGrid& getArg2() {
return vel; }
typedef MACGrid type2;inline Grid<Vec3>& getArg3() {
return tmp; }
typedef Grid<Vec3> type3;inline const ParticleDataImpl<Vec3>& getArg4() {
return pvel; }
typedef ParticleDataImpl<Vec3> type4;inline const ParticleDataImpl<int>* getArg5() {
return ptype; }
typedef ParticleDataImpl<int> type5;inline const int& getArg6() {
return exclude; }
typedef int type6; void runMessage() { debMsg("Executing kernel knMapLinearVec3ToMACGrid ", 3); debMsg("Kernel range" << " size "<< size << " " , 4); }; void run() {
const IndexInt _sz = size; for (IndexInt i = 0; i < _sz; i++) op(i, p,flags,vel,tmp,pvel,ptype,exclude); }
const BasicParticleSystem& p; const FlagGrid& flags; const MACGrid& vel; Grid<Vec3>& tmp; const ParticleDataImpl<Vec3>& pvel; const ParticleDataImpl<int>* ptype; const int exclude; }
;
// optionally , this function can use an existing vec3 grid to store the weights
// this is useful in combination with the simple extrapolation function
void mapPartsToMAC(const FlagGrid& flags, MACGrid& vel, MACGrid& velOld, const BasicParticleSystem& parts, const ParticleDataImpl<Vec3>& partVel, Grid<Vec3>* weight=NULL, const ParticleDataImpl<int>* ptype=NULL, const int exclude=0) {
// interpol -> grid. tmpgrid for particle contribution weights
bool freeTmp = false;
if(!weight) {
weight = new Grid<Vec3>(flags.getParent());
freeTmp = true;
} else {
weight->clear(); // make sure we start with a zero grid!
}
vel.clear();
knMapLinearVec3ToMACGrid( parts, flags, vel, *weight, partVel, ptype, exclude );
// stomp small values in weight to zero to prevent roundoff errors
weight->stomp(Vec3(VECTOR_EPSILON));
vel.safeDivide(*weight);
// store original state
velOld.copyFrom( vel );
if(freeTmp) delete weight;
} static PyObject* _W_12 (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, "mapPartsToMAC" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); MACGrid& vel = *_args.getPtr<MACGrid >("vel",1,&_lock); MACGrid& velOld = *_args.getPtr<MACGrid >("velOld",2,&_lock); const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",3,&_lock); const ParticleDataImpl<Vec3>& partVel = *_args.getPtr<ParticleDataImpl<Vec3> >("partVel",4,&_lock); Grid<Vec3>* weight = _args.getPtrOpt<Grid<Vec3> >("weight",5,NULL,&_lock); const ParticleDataImpl<int>* ptype = _args.getPtrOpt<ParticleDataImpl<int> >("ptype",6,NULL,&_lock); const int exclude = _args.getOpt<int >("exclude",7,0,&_lock); _retval = getPyNone(); mapPartsToMAC(flags,vel,velOld,parts,partVel,weight,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"mapPartsToMAC", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("mapPartsToMAC",e.what()); return 0; }
}
static const Pb::Register _RP_mapPartsToMAC ("","mapPartsToMAC",_W_12);
extern "C" {
void PbRegister_mapPartsToMAC() {
KEEP_UNUSED(_RP_mapPartsToMAC); }
}
template <class T> struct knMapLinear : public KernelBase {
knMapLinear(const BasicParticleSystem& p, const FlagGrid& flags, const Grid<T>& target, Grid<Real>& gtmp, const ParticleDataImpl<T>& psource ) : KernelBase(p.size()) ,p(p),flags(flags),target(target),gtmp(gtmp),psource(psource) {
runMessage(); run(); }
inline void op(IndexInt idx, const BasicParticleSystem& p, const FlagGrid& flags, const Grid<T>& target, Grid<Real>& gtmp, const ParticleDataImpl<T>& psource ) {
unusedParameter(flags);
if (!p.isActive(idx)) return;
target.setInterpolated( p[idx].pos, psource[idx], gtmp );
} inline const BasicParticleSystem& getArg0() {
return p; }
typedef BasicParticleSystem type0;inline const FlagGrid& getArg1() {
return flags; }
typedef FlagGrid type1;inline const Grid<T>& getArg2() {
return target; }
typedef Grid<T> type2;inline Grid<Real>& getArg3() {
return gtmp; }
typedef Grid<Real> type3;inline const ParticleDataImpl<T>& getArg4() {
return psource; }
typedef ParticleDataImpl<T> type4; void runMessage() { debMsg("Executing kernel knMapLinear ", 3); debMsg("Kernel range" << " size "<< size << " " , 4); }; void run() {
const IndexInt _sz = size; for (IndexInt i = 0; i < _sz; i++) op(i, p,flags,target,gtmp,psource); }
const BasicParticleSystem& p; const FlagGrid& flags; const Grid<T>& target; Grid<Real>& gtmp; const ParticleDataImpl<T>& psource; }
;
template<class T>
void mapLinearRealHelper(const FlagGrid& flags, Grid<T>& target,
const BasicParticleSystem& parts, const ParticleDataImpl<T>& source )
{
Grid<Real> tmp(flags.getParent());
target.clear();
knMapLinear<T>( parts, flags, target, tmp, source );
knSafeDivReal<T>( target, tmp );
}
void mapPartsToGrid(const FlagGrid& flags, Grid<Real>& target , const BasicParticleSystem& parts , const ParticleDataImpl<Real>& source ) {
mapLinearRealHelper<Real>(flags,target,parts,source);
} static PyObject* _W_13 (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, "mapPartsToGrid" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); Grid<Real>& target = *_args.getPtr<Grid<Real> >("target",1,&_lock); const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",2,&_lock); const ParticleDataImpl<Real>& source = *_args.getPtr<ParticleDataImpl<Real> >("source",3,&_lock); _retval = getPyNone(); mapPartsToGrid(flags,target,parts,source); _args.check(); }
pbFinalizePlugin(parent,"mapPartsToGrid", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("mapPartsToGrid",e.what()); return 0; }
}
static const Pb::Register _RP_mapPartsToGrid ("","mapPartsToGrid",_W_13);
extern "C" {
void PbRegister_mapPartsToGrid() {
KEEP_UNUSED(_RP_mapPartsToGrid); }
}
void mapPartsToGridVec3(const FlagGrid& flags, Grid<Vec3>& target , const BasicParticleSystem& parts , const ParticleDataImpl<Vec3>& source ) {
mapLinearRealHelper<Vec3>(flags,target,parts,source);
} static PyObject* _W_14 (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, "mapPartsToGridVec3" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); Grid<Vec3>& target = *_args.getPtr<Grid<Vec3> >("target",1,&_lock); const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",2,&_lock); const ParticleDataImpl<Vec3>& source = *_args.getPtr<ParticleDataImpl<Vec3> >("source",3,&_lock); _retval = getPyNone(); mapPartsToGridVec3(flags,target,parts,source); _args.check(); }
pbFinalizePlugin(parent,"mapPartsToGridVec3", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("mapPartsToGridVec3",e.what()); return 0; }
}
static const Pb::Register _RP_mapPartsToGridVec3 ("","mapPartsToGridVec3",_W_14);
extern "C" {
void PbRegister_mapPartsToGridVec3() {
KEEP_UNUSED(_RP_mapPartsToGridVec3); }
}
// integers need "max" mode, not yet implemented
//PYTHON() void mapPartsToGridInt ( FlagGrid& flags, Grid<int >& target , BasicParticleSystem& parts , ParticleDataImpl<int >& source ) {
// mapLinearRealHelper<int >(flags,target,parts,source);
//}
template <class T> struct knMapFromGrid : public KernelBase {
knMapFromGrid(const BasicParticleSystem& p, const Grid<T>& gsrc, ParticleDataImpl<T>& target ) : KernelBase(p.size()) ,p(p),gsrc(gsrc),target(target) {
runMessage(); run(); }
inline void op(IndexInt idx, const BasicParticleSystem& p, const Grid<T>& gsrc, ParticleDataImpl<T>& target ) {
if (!p.isActive(idx)) return;
target[idx] = gsrc.getInterpolated( p[idx].pos );
} inline const BasicParticleSystem& getArg0() {
return p; }
typedef BasicParticleSystem type0;inline const Grid<T>& getArg1() {
return gsrc; }
typedef Grid<T> type1;inline ParticleDataImpl<T>& getArg2() {
return target; }
typedef ParticleDataImpl<T> type2; void runMessage() { debMsg("Executing kernel knMapFromGrid ", 3); debMsg("Kernel range" << " size "<< size << " " , 4); }; void run() {
const IndexInt _sz = size;
#pragma omp parallel
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,p,gsrc,target); }
}
const BasicParticleSystem& p; const Grid<T>& gsrc; ParticleDataImpl<T>& target; }
;
void mapGridToParts(const Grid<Real>& source , const BasicParticleSystem& parts , ParticleDataImpl<Real>& target ) {
knMapFromGrid<Real>(parts, source, target);
} static PyObject* _W_15 (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, "mapGridToParts" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const Grid<Real>& source = *_args.getPtr<Grid<Real> >("source",0,&_lock); const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",1,&_lock); ParticleDataImpl<Real>& target = *_args.getPtr<ParticleDataImpl<Real> >("target",2,&_lock); _retval = getPyNone(); mapGridToParts(source,parts,target); _args.check(); }
pbFinalizePlugin(parent,"mapGridToParts", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("mapGridToParts",e.what()); return 0; }
}
static const Pb::Register _RP_mapGridToParts ("","mapGridToParts",_W_15);
extern "C" {
void PbRegister_mapGridToParts() {
KEEP_UNUSED(_RP_mapGridToParts); }
}
void mapGridToPartsVec3(const Grid<Vec3>& source , const BasicParticleSystem& parts , ParticleDataImpl<Vec3>& target ) {
knMapFromGrid<Vec3>(parts, source, target);
} static PyObject* _W_16 (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, "mapGridToPartsVec3" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const Grid<Vec3>& source = *_args.getPtr<Grid<Vec3> >("source",0,&_lock); const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",1,&_lock); ParticleDataImpl<Vec3>& target = *_args.getPtr<ParticleDataImpl<Vec3> >("target",2,&_lock); _retval = getPyNone(); mapGridToPartsVec3(source,parts,target); _args.check(); }
pbFinalizePlugin(parent,"mapGridToPartsVec3", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("mapGridToPartsVec3",e.what()); return 0; }
}
static const Pb::Register _RP_mapGridToPartsVec3 ("","mapGridToPartsVec3",_W_16);
extern "C" {
void PbRegister_mapGridToPartsVec3() {
KEEP_UNUSED(_RP_mapGridToPartsVec3); }
}
// Get velocities from grid
struct knMapLinearMACGridToVec3_PIC : public KernelBase {
knMapLinearMACGridToVec3_PIC(const BasicParticleSystem& p, const FlagGrid& flags, const MACGrid& vel, ParticleDataImpl<Vec3>& pvel, const ParticleDataImpl<int>* ptype, const int exclude) : KernelBase(p.size()) ,p(p),flags(flags),vel(vel),pvel(pvel),ptype(ptype),exclude(exclude) {
runMessage(); run(); }
inline void op(IndexInt idx, const BasicParticleSystem& p, const FlagGrid& flags, const MACGrid& vel, ParticleDataImpl<Vec3>& pvel, const ParticleDataImpl<int>* ptype, const int exclude ) {
if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude))) return;
// pure PIC
pvel[idx] = vel.getInterpolated( p[idx].pos );
} inline const BasicParticleSystem& getArg0() {
return p; }
typedef BasicParticleSystem type0;inline const FlagGrid& getArg1() {
return flags; }
typedef FlagGrid type1;inline const MACGrid& getArg2() {
return vel; }
typedef MACGrid type2;inline ParticleDataImpl<Vec3>& getArg3() {
return pvel; }
typedef ParticleDataImpl<Vec3> type3;inline const ParticleDataImpl<int>* getArg4() {
return ptype; }
typedef ParticleDataImpl<int> type4;inline const int& getArg5() {
return exclude; }
typedef int type5; void runMessage() { debMsg("Executing kernel knMapLinearMACGridToVec3_PIC ", 3); debMsg("Kernel range" << " size "<< size << " " , 4); }; void run() {
const IndexInt _sz = size;
#pragma omp parallel
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,p,flags,vel,pvel,ptype,exclude); }
}
const BasicParticleSystem& p; const FlagGrid& flags; const MACGrid& vel; ParticleDataImpl<Vec3>& pvel; const ParticleDataImpl<int>* ptype; const int exclude; }
;
void mapMACToParts(const FlagGrid& flags, const MACGrid& vel , const BasicParticleSystem& parts , ParticleDataImpl<Vec3>& partVel, const ParticleDataImpl<int>* ptype=NULL, const int exclude=0) {
knMapLinearMACGridToVec3_PIC( parts, flags, vel, partVel, ptype, exclude );
} static PyObject* _W_17 (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, "mapMACToParts" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); const MACGrid& vel = *_args.getPtr<MACGrid >("vel",1,&_lock); const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",2,&_lock); ParticleDataImpl<Vec3>& partVel = *_args.getPtr<ParticleDataImpl<Vec3> >("partVel",3,&_lock); const ParticleDataImpl<int>* ptype = _args.getPtrOpt<ParticleDataImpl<int> >("ptype",4,NULL,&_lock); const int exclude = _args.getOpt<int >("exclude",5,0,&_lock); _retval = getPyNone(); mapMACToParts(flags,vel,parts,partVel,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"mapMACToParts", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("mapMACToParts",e.what()); return 0; }
}
static const Pb::Register _RP_mapMACToParts ("","mapMACToParts",_W_17);
extern "C" {
void PbRegister_mapMACToParts() {
KEEP_UNUSED(_RP_mapMACToParts); }
}
// with flip delta interpolation
struct knMapLinearMACGridToVec3_FLIP : public KernelBase {
knMapLinearMACGridToVec3_FLIP(const BasicParticleSystem& p, const FlagGrid& flags, const MACGrid& vel, const MACGrid& oldVel, ParticleDataImpl<Vec3>& pvel, const Real flipRatio, const ParticleDataImpl<int>* ptype, const int exclude) : KernelBase(p.size()) ,p(p),flags(flags),vel(vel),oldVel(oldVel),pvel(pvel),flipRatio(flipRatio),ptype(ptype),exclude(exclude) {
runMessage(); run(); }
inline void op(IndexInt idx, const BasicParticleSystem& p, const FlagGrid& flags, const MACGrid& vel, const MACGrid& oldVel, ParticleDataImpl<Vec3>& pvel, const Real flipRatio, const ParticleDataImpl<int>* ptype, const int exclude ) {
if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude))) return;
Vec3 v = vel.getInterpolated(p[idx].pos);
Vec3 delta = v - oldVel.getInterpolated(p[idx].pos);
pvel[idx] = flipRatio * (pvel[idx] + delta) + (1.0 - flipRatio) * v;
} inline const BasicParticleSystem& getArg0() {
return p; }
typedef BasicParticleSystem type0;inline const FlagGrid& getArg1() {
return flags; }
typedef FlagGrid type1;inline const MACGrid& getArg2() {
return vel; }
typedef MACGrid type2;inline const MACGrid& getArg3() {
return oldVel; }
typedef MACGrid type3;inline ParticleDataImpl<Vec3>& getArg4() {
return pvel; }
typedef ParticleDataImpl<Vec3> type4;inline const Real& getArg5() {
return flipRatio; }
typedef Real type5;inline const ParticleDataImpl<int>* getArg6() {
return ptype; }
typedef ParticleDataImpl<int> type6;inline const int& getArg7() {
return exclude; }
typedef int type7; void runMessage() { debMsg("Executing kernel knMapLinearMACGridToVec3_FLIP ", 3); debMsg("Kernel range" << " size "<< size << " " , 4); }; void run() {
const IndexInt _sz = size;
#pragma omp parallel
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,p,flags,vel,oldVel,pvel,flipRatio,ptype,exclude); }
}
const BasicParticleSystem& p; const FlagGrid& flags; const MACGrid& vel; const MACGrid& oldVel; ParticleDataImpl<Vec3>& pvel; const Real flipRatio; const ParticleDataImpl<int>* ptype; const int exclude; }
;
void flipVelocityUpdate(const FlagGrid& flags, const MACGrid& vel, const MACGrid& velOld, const BasicParticleSystem& parts, ParticleDataImpl<Vec3>& partVel, const Real flipRatio, const ParticleDataImpl<int>* ptype=NULL, const int exclude=0) {
knMapLinearMACGridToVec3_FLIP( parts, flags, vel, velOld, partVel, flipRatio, ptype, exclude );
} static PyObject* _W_18 (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, "flipVelocityUpdate" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); const MACGrid& vel = *_args.getPtr<MACGrid >("vel",1,&_lock); const MACGrid& velOld = *_args.getPtr<MACGrid >("velOld",2,&_lock); const BasicParticleSystem& parts = *_args.getPtr<BasicParticleSystem >("parts",3,&_lock); ParticleDataImpl<Vec3>& partVel = *_args.getPtr<ParticleDataImpl<Vec3> >("partVel",4,&_lock); const Real flipRatio = _args.get<Real >("flipRatio",5,&_lock); const ParticleDataImpl<int>* ptype = _args.getPtrOpt<ParticleDataImpl<int> >("ptype",6,NULL,&_lock); const int exclude = _args.getOpt<int >("exclude",7,0,&_lock); _retval = getPyNone(); flipVelocityUpdate(flags,vel,velOld,parts,partVel,flipRatio,ptype,exclude); _args.check(); }
pbFinalizePlugin(parent,"flipVelocityUpdate", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("flipVelocityUpdate",e.what()); return 0; }
}
static const Pb::Register _RP_flipVelocityUpdate ("","flipVelocityUpdate",_W_18);
extern "C" {
void PbRegister_flipVelocityUpdate() {
KEEP_UNUSED(_RP_flipVelocityUpdate); }
}
//******************************************************************************
// narrow band
struct knCombineVels : public KernelBase {
knCombineVels(MACGrid& vel, const Grid<Vec3>& w, MACGrid& combineVel, const LevelsetGrid* phi, Real narrowBand, Real thresh ) : KernelBase(&vel,0) ,vel(vel),w(w),combineVel(combineVel),phi(phi),narrowBand(narrowBand),thresh(thresh) {
runMessage(); run(); }
inline void op(int i, int j, int k, MACGrid& vel, const Grid<Vec3>& w, MACGrid& combineVel, const LevelsetGrid* phi, Real narrowBand, Real thresh ) {
int idx = vel.index(i,j,k);
for(int c=0; c<3; ++c)
{
// Correct narrow-band FLIP
if(phi) {
Vec3 pos(i,j,k);
pos[(c+1)%3] += Real(0.5);
pos[(c+2)%3] += Real(0.5);
Real p = phi->getInterpolated(pos);
if (p < -narrowBand) { vel[idx][c] = 0; continue; }
}
if (w[idx][c] > thresh) {
combineVel[idx][c] = vel[idx][c];
vel[idx][c] = -1;
} else {
vel[idx][c] = 0;
}
}
} inline MACGrid& getArg0() {
return vel; }
typedef MACGrid type0;inline const Grid<Vec3>& getArg1() {
return w; }
typedef Grid<Vec3> type1;inline MACGrid& getArg2() {
return combineVel; }
typedef MACGrid type2;inline const LevelsetGrid* getArg3() {
return phi; }
typedef LevelsetGrid type3;inline Real& getArg4() {
return narrowBand; }
typedef Real type4;inline Real& getArg5() {
return thresh; }
typedef Real type5; void runMessage() { debMsg("Executing kernel knCombineVels ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; k++) for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,vel,w,combineVel,phi,narrowBand,thresh); }
}
else {
const int k=0;
#pragma omp parallel
{
#pragma omp for
for (int j=0; j < _maxY; j++) for (int i=0; i < _maxX; i++) op(i,j,k,vel,w,combineVel,phi,narrowBand,thresh); }
}
}
MACGrid& vel; const Grid<Vec3>& w; MACGrid& combineVel; const LevelsetGrid* phi; Real narrowBand; Real thresh; }
;
//! narrow band velocity combination
void combineGridVel( MACGrid& vel, const Grid<Vec3>& weight, MACGrid& combineVel, const LevelsetGrid* phi=NULL, Real narrowBand=0.0, Real thresh=0.0) {
knCombineVels(vel, weight, combineVel, phi, narrowBand, thresh);
} static PyObject* _W_19 (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, "combineGridVel" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; MACGrid& vel = *_args.getPtr<MACGrid >("vel",0,&_lock); const Grid<Vec3>& weight = *_args.getPtr<Grid<Vec3> >("weight",1,&_lock); MACGrid& combineVel = *_args.getPtr<MACGrid >("combineVel",2,&_lock); const LevelsetGrid* phi = _args.getPtrOpt<LevelsetGrid >("phi",3,NULL,&_lock); Real narrowBand = _args.getOpt<Real >("narrowBand",4,0.0,&_lock); Real thresh = _args.getOpt<Real >("thresh",5,0.0,&_lock); _retval = getPyNone(); combineGridVel(vel,weight,combineVel,phi,narrowBand,thresh); _args.check(); }
pbFinalizePlugin(parent,"combineGridVel", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("combineGridVel",e.what()); return 0; }
}
static const Pb::Register _RP_combineGridVel ("","combineGridVel",_W_19);
extern "C" {
void PbRegister_combineGridVel() {
KEEP_UNUSED(_RP_combineGridVel); }
}
//! surface tension helper
void getLaplacian(Grid<Real> &laplacian, const Grid<Real> &grid) {
LaplaceOp(laplacian, grid);
} static PyObject* _W_20 (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, "getLaplacian" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; Grid<Real> & laplacian = *_args.getPtr<Grid<Real> >("laplacian",0,&_lock); const Grid<Real> & grid = *_args.getPtr<Grid<Real> >("grid",1,&_lock); _retval = getPyNone(); getLaplacian(laplacian,grid); _args.check(); }
pbFinalizePlugin(parent,"getLaplacian", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("getLaplacian",e.what()); return 0; }
}
static const Pb::Register _RP_getLaplacian ("","getLaplacian",_W_20);
extern "C" {
void PbRegister_getLaplacian() {
KEEP_UNUSED(_RP_getLaplacian); }
}
void getCurvature(Grid<Real> &curv, const Grid<Real> &grid, const Real h=1.0) {
CurvatureOp(curv, grid, h);
} static PyObject* _W_21 (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, "getCurvature" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; Grid<Real> & curv = *_args.getPtr<Grid<Real> >("curv",0,&_lock); const Grid<Real> & grid = *_args.getPtr<Grid<Real> >("grid",1,&_lock); const Real h = _args.getOpt<Real >("h",2,1.0,&_lock); _retval = getPyNone(); getCurvature(curv,grid,h); _args.check(); }
pbFinalizePlugin(parent,"getCurvature", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("getCurvature",e.what()); return 0; }
}
static const Pb::Register _RP_getCurvature ("","getCurvature",_W_21);
extern "C" {
void PbRegister_getCurvature() {
KEEP_UNUSED(_RP_getCurvature); }
}
} // namespace