Page MenuHome
Differential D3850 Diff 17619 intern/mantaflow/intern/manta_develop/preprocessed/tbb/fastmarch.cpp

Changeset View

Standalone View

View Options

intern/mantaflow/intern/manta_develop/preprocessed/tbb/fastmarch.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
*
* Fast marching and extrapolation
*
******************************************************************************/
#include "fastmarch.h"
#include "levelset.h"
#include "kernel.h"
#include <algorithm>
using namespace std;
namespace Manta {
template<class COMP, int TDIR>
FastMarch<COMP,TDIR>::FastMarch(const FlagGrid& flags, Grid<int>& fmFlags, Grid<Real>& levelset, Real maxTime, MACGrid* velTransport )
: mLevelset(levelset), mFlags(flags), mFmFlags(fmFlags)
{
if (velTransport)
mVelTransport.initMarching(velTransport, &flags);
mMaxTime = maxTime * TDIR;
}
// helper for individual components to calculateDistance
template<class COMP, int TDIR> template<int C>
Real FastMarch<COMP,TDIR>::calcWeights(int& okcnt, int& invcnt, Real* v, const Vec3i& idx) {
Real val = 0.;
Vec3i idxPlus(idx), idxMinus(idx);
idxPlus[C]++;
idxMinus[C]--;
mWeights[C*2] = mWeights[C*2+1] = 0.;
if (mFmFlags(idxPlus)==FlagInited) {
// somewhat arbitrary - choose +1 value over -1 ...
val = mLevelset(idxPlus);
v[okcnt] = val; okcnt++;
mWeights[C*2] = 1.;
} else if (mFmFlags(idxMinus)==FlagInited) {
val = mLevelset(idxMinus);
v[okcnt] = val; okcnt++;
mWeights[C*2+1] = 1.;
}
else {
invcnt++;
}
return val;
}
template<class COMP, int TDIR>
inline Real FastMarch<COMP,TDIR>::calculateDistance(const Vec3i& idx) {
//int invflag = 0;
int invcnt = 0;
Real v[3];
int okcnt = 0;
Real aVal = calcWeights<0>(okcnt, invcnt, v, idx);
Real bVal = calcWeights<1>(okcnt, invcnt, v, idx);
Real cVal = 0.;
if (mLevelset.is3D()) cVal = calcWeights<2>(okcnt, invcnt, v, idx);
else { invcnt++; mWeights[4] = mWeights[5] = 0.; }
Real ret = InvalidTime();
switch(invcnt) {
case 0: {
// take all values
const Real ca=v[0], cb=v[1], cc=v[2];
const Real csqrt = max(0. ,
-2.*(ca*ca+cb*cb- cb*cc + cc*cc - ca*(cb+cc)) + 3 );
// clamp to make sure the sqrt is valid
ret = 0.333333*( ca+cb+cc+ TDIR*sqrt(csqrt) );
// weights needed for transport (transpTouch)
mWeights[0] *= fabs(ret-ca);
mWeights[1] *= fabs(ret-ca);
mWeights[2] *= fabs(ret-cb);
mWeights[3] *= fabs(ret-cb);
mWeights[4] *= fabs(ret-cc);
mWeights[5] *= fabs(ret-cc);
Real norm = 0.0; // try to force normalization
for(int i=0;i<6;i++) {
norm += mWeights[i];
}
norm = 1.0/norm;
for(int i=0;i<6;i++) { mWeights[i] *= norm; }
} break;
case 1: {
// take just the 2 ok values
// t=0.5*( a+b+ (2*g*g-(b-a)*(b-a))^0.5)
const Real csqrt = max(0. , 2.-(v[1]-v[0])*(v[1]-v[0]) );
// clamp to make sure the sqrt is valid
ret = 0.5*( v[0]+v[1]+ TDIR*sqrt(csqrt) );
// weights needed for transport (transpTouch)
mWeights[0] *= fabs(ret-aVal);
mWeights[1] *= fabs(ret-aVal);
mWeights[2] *= fabs(ret-bVal);
mWeights[3] *= fabs(ret-bVal);
mWeights[4] *= fabs(ret-cVal);
mWeights[5] *= fabs(ret-cVal);
Real norm = 0.0; // try to force normalization
for(int i=0;i<6;i++) {
norm += mWeights[i];
}
norm = 1.0/norm;
for(int i=0;i<6;i++) { mWeights[i] *= norm; }
// */
} break;
case 2: {
// just use the one remaining value
ret = v[0]+ (Real)(TDIR) ; // direction = +- 1
} break;
default:
errMsg("FastMarch :: Invalid invcnt");
break;
}
return ret;
}
template<class COMP, int TDIR>
void FastMarch<COMP,TDIR>::addToList(const Vec3i& p, const Vec3i& src) {
if (!mLevelset.isInBounds(p,1)) return;
const IndexInt idx = mLevelset.index(p);
// already known value, value alreay set to valid value? skip cell...
if(mFmFlags[idx] == FlagInited) return;
// discard by source time now , TODO do instead before calling all addtolists?
Real srct = mLevelset(src);
if(COMP::compare(srct, mMaxTime)) return;
Real ttime = calculateDistance(p);
// remove old entry if larger
bool found=false;
Real oldt = mLevelset[idx];
if (mFmFlags[idx] == FlagIsOnHeap) {
found = true;
// is old time better?
if(COMP::compare(ttime,oldt)) return;
}
// update field
mFmFlags[idx] = FlagIsOnHeap;
mLevelset[idx] = ttime;
// debug info std::cout<<"set "<< idx <<","<< ttime <<"\n";
if (mVelTransport.isInitialized())
mVelTransport.transpTouch(p.x, p.y, p.z, mWeights, ttime);
// the following adds entries to the heap of active cells
// current: (!found) , previous: always add, might lead to duplicate
// entries, but the earlier will be handled earlier, the second one will skip to the FlagInited check above
if(!found)
{
// add list entry with source value
COMP entry;
entry.p = p;
entry.time = mLevelset[idx];
mHeap.push( entry );
// debug info std::cout<<"push "<< entry.p <<","<< entry.time <<"\n";
}
}
//! Enforce delta_phi = 0 on boundaries
struct SetLevelsetBoundaries : public KernelBase {
SetLevelsetBoundaries(Grid<Real>& phi) : KernelBase(&phi,0) ,phi(phi) {
runMessage(); run(); }
inline void op(int i, int j, int k, Grid<Real>& phi ) {
if (i==0) phi(i,j,k) = phi(1,j,k);
if (i==maxX-1) phi(i,j,k) = phi(i-1,j,k);
if (j==0) phi(i,j,k) = phi(i,1,k);
if (j==maxY-1) phi(i,j,k) = phi(i,j-1,k);
if(phi.is3D()) {
if (k==0) phi(i,j,k) = phi(i,j,1);
if (k==maxZ-1) phi(i,j,k) = phi(i,j,k-1);
}
} inline Grid<Real>& getArg0() {
return phi; }
typedef Grid<Real> type0; void runMessage() { debMsg("Executing kernel SetLevelsetBoundaries ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; 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, phi); }
Grid<Real>& phi; }
;
/*****************************************************************************/
//! Walk...
template<class COMP, int TDIR>
void FastMarch<COMP,TDIR>::performMarching() {
mReheapVal = 0.0;
while(mHeap.size() > 0) {
const COMP& ce = mHeap.top();
Vec3i p = ce.p;
mFmFlags(p) = FlagInited;
mHeap.pop();
// debug info std::cout<<"pop "<< ce.p <<","<< ce.time <<"\n";
addToList(Vec3i(p.x-1,p.y,p.z), p);
addToList(Vec3i(p.x+1,p.y,p.z), p);
addToList(Vec3i(p.x,p.y-1,p.z), p);
addToList(Vec3i(p.x,p.y+1,p.z), p);
if(mLevelset.is3D()) {
addToList(Vec3i(p.x,p.y,p.z-1), p);
addToList(Vec3i(p.x,p.y,p.z+1), p);
}
}
// set boundary for plain array
SetLevelsetBoundaries setls(mLevelset);
setls.getArg0(); // get rid of compiler warning...
}
// explicit instantiation
template class FastMarch<FmHeapEntryIn, -1>;
template class FastMarch<FmHeapEntryOut, +1>;
/*****************************************************************************/
// simpler extrapolation functions (primarily for FLIP)
struct knExtrapolateMACSimple : public KernelBase {
knExtrapolateMACSimple(MACGrid& vel, int distance , Grid<int>& tmp , const int d , const int c ) : KernelBase(&vel,1) ,vel(vel),distance(distance),tmp(tmp),d(d),c(c) {
runMessage(); run(); }
inline void op(int i, int j, int k, MACGrid& vel, int distance , Grid<int>& tmp , const int d , const int c ) const {
static const Vec3i nb[6] = {
Vec3i(1 ,0,0), Vec3i(-1,0,0),
Vec3i(0,1 ,0), Vec3i(0,-1,0),
Vec3i(0,0,1 ), Vec3i(0,0,-1) };
const int dim = (vel.is3D() ? 3:2);
if (tmp(i,j,k) != 0) return;
// copy from initialized neighbors
Vec3i p(i,j,k);
int nbs = 0;
Real avgVel = 0.;
for (int n=0; n<2*dim; ++n) {
if (tmp(p+nb[n]) == d) {
//vel(p)[c] = (c+1.)*0.1;
avgVel += vel(p+nb[n])[c];
nbs++;
}
}
if(nbs>0) {
tmp(p) = d+1;
vel(p)[c] = avgVel / nbs;
}
} inline MACGrid& getArg0() {
return vel; }
typedef MACGrid type0;inline int& getArg1() {
return distance; }
typedef int type1;inline Grid<int>& getArg2() {
return tmp; }
typedef Grid<int> type2;inline const int& getArg3() {
return d; }
typedef int type3;inline const int& getArg4() {
return c; }
typedef int type4; void runMessage() { debMsg("Executing kernel knExtrapolateMACSimple ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void operator() (const tbb::blocked_range<IndexInt>& __r) const {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ>1) {
for (int k=__r.begin(); k!=(int)__r.end(); k++) for (int j=1; j<_maxY; j++) for (int i=1; i<_maxX; i++) op(i,j,k,vel,distance,tmp,d,c); }
else {
const int k=0; for (int j=__r.begin(); j!=(int)__r.end(); j++) for (int i=1; i<_maxX; i++) op(i,j,k,vel,distance,tmp,d,c); }
}
void run() {
if (maxZ>1) tbb::parallel_for (tbb::blocked_range<IndexInt>(minZ, maxZ), *this); else tbb::parallel_for (tbb::blocked_range<IndexInt>(1, maxY), *this); }
MACGrid& vel; int distance; Grid<int>& tmp; const int d; const int c; }
;
//! copy velocity into domain side, note - don't read & write same grid, hence velTmp copy
struct knExtrapolateIntoBnd : public KernelBase {
knExtrapolateIntoBnd(FlagGrid& flags, MACGrid& vel, const MACGrid& velTmp) : KernelBase(&flags,0) ,flags(flags),vel(vel),velTmp(velTmp) {
runMessage(); run(); }
inline void op(int i, int j, int k, FlagGrid& flags, MACGrid& vel, const MACGrid& velTmp ) const {
int c=0;
Vec3 v(0,0,0);
const bool isObs = flags.isObstacle(i,j,k);
if( i==0 ) {
v = velTmp(i+1,j,k);
if(isObs && v[0] < 0.) v[0] = 0.;
c++;
}
else if( i==(flags.getSizeX()-1) ) {
v = velTmp(i-1,j,k);
if(isObs && v[0] > 0.) v[0] = 0.;
c++;
}
if( j==0 ) {
v = velTmp(i,j+1,k);
if(isObs && v[1] < 0.) v[1] = 0.;
c++;
}
else if( j==(flags.getSizeY()-1) ) {
v = velTmp(i,j-1,k);
if(isObs && v[1] > 0.) v[1] = 0.;
c++;
}
if(flags.is3D()) {
if( k==0 ) {
v = velTmp(i,j,k+1);
if(isObs && v[2] < 0.) v[2] = 0.;
c++;
}
else if( k==(flags.getSizeZ()-1) ) {
v = velTmp(i,j,k-1);
if(isObs && v[2] > 0.) v[2] = 0.;
c++;
} }
if(c>0) {
vel(i,j,k) = v/(Real)c;
}
} inline FlagGrid& getArg0() {
return flags; }
typedef FlagGrid type0;inline MACGrid& getArg1() {
return vel; }
typedef MACGrid type1;inline const MACGrid& getArg2() {
return velTmp; }
typedef MACGrid type2; void runMessage() { debMsg("Executing kernel knExtrapolateIntoBnd ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void operator() (const tbb::blocked_range<IndexInt>& __r) const {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ>1) {
for (int k=__r.begin(); k!=(int)__r.end(); k++) for (int j=0; j<_maxY; j++) for (int i=0; i<_maxX; i++) op(i,j,k,flags,vel,velTmp); }
else {
const int k=0; for (int j=__r.begin(); j!=(int)__r.end(); j++) for (int i=0; i<_maxX; i++) op(i,j,k,flags,vel,velTmp); }
}
void run() {
if (maxZ>1) tbb::parallel_for (tbb::blocked_range<IndexInt>(minZ, maxZ), *this); else tbb::parallel_for (tbb::blocked_range<IndexInt>(0, maxY), *this); }
FlagGrid& flags; MACGrid& vel; const MACGrid& velTmp; }
;
// todo - use getGradient instead?
inline Vec3 getNormal(const Grid<Real>& data, int i, int j, int k) {
if (i > data.getSizeX()-2) i= data.getSizeX()-2;
if (i < 1) i = 1;
if (j > data.getSizeY()-2) j= data.getSizeY()-2;
if (j < 1) j = 1;
int kd = 1;
if(data.is3D()) {
if (k > data.getSizeZ()-2) k= data.getSizeZ()-2;
if (k < 1) k = 1;
} else { kd=0; }
return Vec3( data(i+1,j ,k ) - data(i-1,j ,k ) ,
data(i ,j+1,k ) - data(i ,j-1,k ) ,
data(i ,j ,k+kd) - data(i ,j ,k-kd) );
}
struct knUnprojectNormalComp : public KernelBase {
knUnprojectNormalComp(FlagGrid& flags, MACGrid& vel, Grid<Real>& phi, Real maxDist) : KernelBase(&flags,1) ,flags(flags),vel(vel),phi(phi),maxDist(maxDist) {
runMessage(); run(); }
inline void op(int i, int j, int k, FlagGrid& flags, MACGrid& vel, Grid<Real>& phi, Real maxDist ) const {
// apply inside, within range near obstacle surface
if(phi(i,j,k)>0. || phi(i,j,k)<-maxDist) return;
Vec3 n = getNormal(phi, i,j,k);
Vec3 v = vel(i,j,k);
if(dot(n,v) < 0.) {
normalize(n);
Real l = dot(n,v);
vel(i,j,k) -= n*l;
}
} inline FlagGrid& getArg0() {
return flags; }
typedef FlagGrid type0;inline MACGrid& getArg1() {
return vel; }
typedef MACGrid type1;inline Grid<Real>& getArg2() {
return phi; }
typedef Grid<Real> type2;inline Real& getArg3() {
return maxDist; }
typedef Real type3; void runMessage() { debMsg("Executing kernel knUnprojectNormalComp ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void operator() (const tbb::blocked_range<IndexInt>& __r) const {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ>1) {
for (int k=__r.begin(); k!=(int)__r.end(); k++) for (int j=1; j<_maxY; j++) for (int i=1; i<_maxX; i++) op(i,j,k,flags,vel,phi,maxDist); }
else {
const int k=0; for (int j=__r.begin(); j!=(int)__r.end(); j++) for (int i=1; i<_maxX; i++) op(i,j,k,flags,vel,phi,maxDist); }
}
void run() {
if (maxZ>1) tbb::parallel_for (tbb::blocked_range<IndexInt>(minZ, maxZ), *this); else tbb::parallel_for (tbb::blocked_range<IndexInt>(1, maxY), *this); }
FlagGrid& flags; MACGrid& vel; Grid<Real>& phi; Real maxDist; }
;
// a simple extrapolation step , used for cases where there's no levelset
// (note, less accurate than fast marching extrapolation.)
// into obstacle is a special mode for second order obstable boundaries (extrapolating
// only fluid velocities, not those at obstacles)
void extrapolateMACSimple(FlagGrid& flags, MACGrid& vel, int distance = 4, LevelsetGrid* phiObs=NULL , bool intoObs = false ) {
Grid<int> tmp( flags.getParent() );
int dim = (flags.is3D() ? 3:2);
for(int c=0; c<dim; ++c) {
Vec3i dir = 0;
dir[c] = 1;
tmp.clear();
// remove all fluid cells (not touching obstacles)
FOR_IJK_BND(flags,1) {
Vec3i p(i,j,k);
bool mark = false;
if(!intoObs) {
if( flags.isFluid(p) || flags.isFluid(p-dir) ) mark = true;
} else {
if( (flags.isFluid(p) || flags.isFluid(p-dir) ) &&
(!flags.isObstacle(p)) && (!flags.isObstacle(p-dir)) ) mark = true;
}
if(mark) tmp(p) = 1;
}
// extrapolate for distance
for(int d=1; d<1+distance; ++d) {
knExtrapolateMACSimple(vel, distance, tmp, d, c);
} // d
}
if(phiObs) {
knUnprojectNormalComp( flags, vel, *phiObs, distance );
}
// copy tangential values into sides of domain
MACGrid velTmp( flags.getParent() ); velTmp.copyFrom(vel);
knExtrapolateIntoBnd(flags, vel, velTmp);
} 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, "extrapolateMACSimple" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); MACGrid& vel = *_args.getPtr<MACGrid >("vel",1,&_lock); int distance = _args.getOpt<int >("distance",2,4,&_lock); LevelsetGrid* phiObs = _args.getPtrOpt<LevelsetGrid >("phiObs",3,NULL ,&_lock); bool intoObs = _args.getOpt<bool >("intoObs",4,false ,&_lock); _retval = getPyNone(); extrapolateMACSimple(flags,vel,distance,phiObs,intoObs); _args.check(); }
pbFinalizePlugin(parent,"extrapolateMACSimple", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("extrapolateMACSimple",e.what()); return 0; }
}
static const Pb::Register _RP_extrapolateMACSimple ("","extrapolateMACSimple",_W_0);
extern "C" {
void PbRegister_extrapolateMACSimple() {
KEEP_UNUSED(_RP_extrapolateMACSimple); }
}
struct knExtrapolateMACFromWeight : public KernelBase {
knExtrapolateMACFromWeight( MACGrid& vel, Grid<Vec3>& weight, int distance , const int d, const int c ) : KernelBase(&vel,1) ,vel(vel),weight(weight),distance(distance),d(d),c(c) {
runMessage(); run(); }
inline void op(int i, int j, int k, MACGrid& vel, Grid<Vec3>& weight, int distance , const int d, const int c ) const {
static const Vec3i nb[6] = {
Vec3i(1 ,0,0), Vec3i(-1,0,0),
Vec3i(0,1 ,0), Vec3i(0,-1,0),
Vec3i(0,0,1 ), Vec3i(0,0,-1) };
const int dim = (vel.is3D() ? 3:2);
if (weight(i,j,k)[c] != 0) return;
// copy from initialized neighbors
Vec3i p(i,j,k);
int nbs = 0;
Real avgVel = 0.;
for (int n=0; n<2*dim; ++n) {
if (weight(p+nb[n])[c] == d) {
avgVel += vel(p+nb[n])[c];
nbs++;
}
}
if(nbs>0) {
weight(p)[c] = d+1;
vel(p)[c] = avgVel / nbs;
}
} inline MACGrid& getArg0() {
return vel; }
typedef MACGrid type0;inline Grid<Vec3>& getArg1() {
return weight; }
typedef Grid<Vec3> type1;inline int& getArg2() {
return distance; }
typedef int type2;inline const int& getArg3() {
return d; }
typedef int type3;inline const int& getArg4() {
return c; }
typedef int type4; void runMessage() { debMsg("Executing kernel knExtrapolateMACFromWeight ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void operator() (const tbb::blocked_range<IndexInt>& __r) const {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ>1) {
for (int k=__r.begin(); k!=(int)__r.end(); k++) for (int j=1; j<_maxY; j++) for (int i=1; i<_maxX; i++) op(i,j,k,vel,weight,distance,d,c); }
else {
const int k=0; for (int j=__r.begin(); j!=(int)__r.end(); j++) for (int i=1; i<_maxX; i++) op(i,j,k,vel,weight,distance,d,c); }
}
void run() {
if (maxZ>1) tbb::parallel_for (tbb::blocked_range<IndexInt>(minZ, maxZ), *this); else tbb::parallel_for (tbb::blocked_range<IndexInt>(1, maxY), *this); }
MACGrid& vel; Grid<Vec3>& weight; int distance; const int d; const int c; }
;
// same as extrapolateMACSimple, but uses weight vec3 grid instead of flags to check
// for valid values (to be used in combination with mapPartsToMAC)
// note - the weight grid values are destroyed! the function is necessary due to discrepancies
// between velocity mapping on surface-levelset / fluid-flag creation. With this
// extrapolation we make sure the fluid region is covered by initial velocities
void extrapolateMACFromWeight( MACGrid& vel, Grid<Vec3>& weight, int distance = 2) {
const int dim = (vel.is3D() ? 3:2);
for(int c=0; c<dim; ++c) {
Vec3i dir = 0;
dir[c] = 1;
// reset weight values to 0 (uninitialized), and 1 (initialized inner values)
FOR_IJK_BND(vel,1) {
Vec3i p(i,j,k);
if(weight(p)[c]>0.) weight(p)[c] = 1.0;
}
// extrapolate for distance
for(int d=1; d<1+distance; ++d) {
knExtrapolateMACFromWeight(vel, weight, distance, d, c);
} // d
}
} 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, "extrapolateMACFromWeight" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; MACGrid& vel = *_args.getPtr<MACGrid >("vel",0,&_lock); Grid<Vec3>& weight = *_args.getPtr<Grid<Vec3> >("weight",1,&_lock); int distance = _args.getOpt<int >("distance",2,2,&_lock); _retval = getPyNone(); extrapolateMACFromWeight(vel,weight,distance); _args.check(); }
pbFinalizePlugin(parent,"extrapolateMACFromWeight", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("extrapolateMACFromWeight",e.what()); return 0; }
}
static const Pb::Register _RP_extrapolateMACFromWeight ("","extrapolateMACFromWeight",_W_1);
extern "C" {
void PbRegister_extrapolateMACFromWeight() {
KEEP_UNUSED(_RP_extrapolateMACFromWeight); }
}
// simple extrapolation functions for levelsets
static const Vec3i nb[6] = {
Vec3i(1 ,0,0), Vec3i(-1,0,0),
Vec3i(0,1 ,0), Vec3i(0,-1,0),
Vec3i(0,0,1 ), Vec3i(0,0,-1) };
template <class S> struct knExtrapolateLsSimple : public KernelBase {
knExtrapolateLsSimple(Grid<S>& val, int distance , Grid<int>& tmp , const int d , S direction ) : KernelBase(&val,1) ,val(val),distance(distance),tmp(tmp),d(d),direction(direction) {
runMessage(); run(); }
inline void op(int i, int j, int k, Grid<S>& val, int distance , Grid<int>& tmp , const int d , S direction ) const {
const int dim = (val.is3D() ? 3:2);
if (tmp(i,j,k) != 0) return;
// copy from initialized neighbors
Vec3i p(i,j,k);
int nbs = 0;
S avg(0.);
for (int n=0; n<2*dim; ++n) {
if (tmp(p+nb[n]) == d) {
avg += val(p+nb[n]);
nbs++;
}
}
if(nbs>0) {
tmp(p) = d+1;
val(p) = avg / nbs + direction;
}
} inline Grid<S>& getArg0() {
return val; }
typedef Grid<S> type0;inline int& getArg1() {
return distance; }
typedef int type1;inline Grid<int>& getArg2() {
return tmp; }
typedef Grid<int> type2;inline const int& getArg3() {
return d; }
typedef int type3;inline S& getArg4() {
return direction; }
typedef S type4; void runMessage() { debMsg("Executing kernel knExtrapolateLsSimple ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void operator() (const tbb::blocked_range<IndexInt>& __r) const {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ>1) {
for (int k=__r.begin(); k!=(int)__r.end(); k++) for (int j=1; j<_maxY; j++) for (int i=1; i<_maxX; i++) op(i,j,k,val,distance,tmp,d,direction); }
else {
const int k=0; for (int j=__r.begin(); j!=(int)__r.end(); j++) for (int i=1; i<_maxX; i++) op(i,j,k,val,distance,tmp,d,direction); }
}
void run() {
if (maxZ>1) tbb::parallel_for (tbb::blocked_range<IndexInt>(minZ, maxZ), *this); else tbb::parallel_for (tbb::blocked_range<IndexInt>(1, maxY), *this); }
Grid<S>& val; int distance; Grid<int>& tmp; const int d; S direction; }
;
template <class S> struct knSetRemaining : public KernelBase {
knSetRemaining(Grid<S>& phi, Grid<int>& tmp, S distance ) : KernelBase(&phi,1) ,phi(phi),tmp(tmp),distance(distance) {
runMessage(); run(); }
inline void op(int i, int j, int k, Grid<S>& phi, Grid<int>& tmp, S distance ) const {
if (tmp(i,j,k) != 0) return;
phi(i,j,k) = distance;
} inline Grid<S>& getArg0() {
return phi; }
typedef Grid<S> type0;inline Grid<int>& getArg1() {
return tmp; }
typedef Grid<int> type1;inline S& getArg2() {
return distance; }
typedef S type2; void runMessage() { debMsg("Executing kernel knSetRemaining ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void operator() (const tbb::blocked_range<IndexInt>& __r) const {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ>1) {
for (int k=__r.begin(); k!=(int)__r.end(); k++) for (int j=1; j<_maxY; j++) for (int i=1; i<_maxX; i++) op(i,j,k,phi,tmp,distance); }
else {
const int k=0; for (int j=__r.begin(); j!=(int)__r.end(); j++) for (int i=1; i<_maxX; i++) op(i,j,k,phi,tmp,distance); }
}
void run() {
if (maxZ>1) tbb::parallel_for (tbb::blocked_range<IndexInt>(minZ, maxZ), *this); else tbb::parallel_for (tbb::blocked_range<IndexInt>(1, maxY), *this); }
Grid<S>& phi; Grid<int>& tmp; S distance; }
;
void extrapolateLsSimple(Grid<Real>& phi, int distance = 4, bool inside=false ) {
Grid<int> tmp( phi.getParent() );
tmp.clear();
const int dim = (phi.is3D() ? 3:2);
// by default, march outside
Real direction = 1.;
if(!inside) {
// mark all inside
FOR_IJK_BND(phi,1) {
if ( phi(i,j,k) < 0. ) { tmp(i,j,k) = 1; }
}
} else {
direction = -1.;
FOR_IJK_BND(phi,1) {
if ( phi(i,j,k) > 0. ) { tmp(i,j,k) = 1; }
}
}
// + first layer around
FOR_IJK_BND(phi,1) {
Vec3i p(i,j,k);
if ( tmp(p) ) continue;
for (int n=0; n<2*dim; ++n) {
if (tmp(p+nb[n])==1) {
tmp(i,j,k) = 2; n=2*dim;
}
}
}
// extrapolate for distance
for(int d=2; d<1+distance; ++d) {
knExtrapolateLsSimple<Real>(phi, distance, tmp, d, direction );
}
// set all remaining cells to max
knSetRemaining<Real>(phi, tmp, Real(direction * (distance+2)) );
} 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, "extrapolateLsSimple" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; Grid<Real>& phi = *_args.getPtr<Grid<Real> >("phi",0,&_lock); int distance = _args.getOpt<int >("distance",1,4,&_lock); bool inside = _args.getOpt<bool >("inside",2,false ,&_lock); _retval = getPyNone(); extrapolateLsSimple(phi,distance,inside); _args.check(); }
pbFinalizePlugin(parent,"extrapolateLsSimple", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("extrapolateLsSimple",e.what()); return 0; }
}
static const Pb::Register _RP_extrapolateLsSimple ("","extrapolateLsSimple",_W_2);
extern "C" {
void PbRegister_extrapolateLsSimple() {
KEEP_UNUSED(_RP_extrapolateLsSimple); }
}
// extrapolate centered vec3 values from marked fluid cells
void extrapolateVec3Simple(Grid<Vec3>& vel, Grid<Real>& phi, int distance = 4, bool inside=false) {
Grid<int> tmp( vel.getParent() );
tmp.clear();
const int dim = (vel.is3D() ? 3:2);
// mark initial cells, by default, march outside
if(!inside) {
// mark all inside
FOR_IJK_BND(phi,1) {
if ( phi(i,j,k) < 0. ) { tmp(i,j,k) = 1; }
}
} else {
FOR_IJK_BND(phi,1) {
if ( phi(i,j,k) > 0. ) { tmp(i,j,k) = 1; }
}
}
// + first layer next to initial cells
FOR_IJK_BND(vel,1) {
Vec3i p(i,j,k);
if ( tmp(p) ) continue;
for (int n=0; n<2*dim; ++n) {
if (tmp(p+nb[n])==1) {
tmp(i,j,k) = 2; n=2*dim;
}
}
}
for(int d=2; d<1+distance; ++d) {
knExtrapolateLsSimple<Vec3>(vel, distance, tmp, d, Vec3(0.) );
}
knSetRemaining<Vec3>(vel, tmp, Vec3(0.) );
} 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, "extrapolateVec3Simple" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; Grid<Vec3>& vel = *_args.getPtr<Grid<Vec3> >("vel",0,&_lock); Grid<Real>& phi = *_args.getPtr<Grid<Real> >("phi",1,&_lock); int distance = _args.getOpt<int >("distance",2,4,&_lock); bool inside = _args.getOpt<bool >("inside",3,false,&_lock); _retval = getPyNone(); extrapolateVec3Simple(vel,phi,distance,inside); _args.check(); }
pbFinalizePlugin(parent,"extrapolateVec3Simple", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("extrapolateVec3Simple",e.what()); return 0; }
}
static const Pb::Register _RP_extrapolateVec3Simple ("","extrapolateVec3Simple",_W_3);
extern "C" {
void PbRegister_extrapolateVec3Simple() {
KEEP_UNUSED(_RP_extrapolateVec3Simple); }
}
} // namespace